# Computational Design of Microarchitected Flow-Through
Electrodes for Energy Storage

###### Abstract

Porous flow-through electrodes are used as the core reactive component across
electrochemical technologies. Controlling the fluid flow, species transport,
and reactive environment is critical to attaining high performance. However,
conventional electrode materials like felts and papers provide few
opportunities for precise engineering of the electrode and its microstructure.
To address these limitations, architected electrodes composed of unit cells
with spatially varying geometry determined via computational optimization are
proposed. Resolved simulation is employed to develop a homogenized description
of the constituent unit cells. These effective properties serve as inputs to a
continuum model for the electrode when used in the negative half cell of a
vanadium redox flow battery. Porosity distributions minimizing power loss are
then determined via computational design optimization to generate architected
porosity electrodes. The architected electrodes are compared to bulk, uniform
porosity electrodes and found to lead to increased power efficiency across
operating flow rates and currents. The design methodology is further used to
generate a scaled-up electrode with comparable power efficiency to the bench-
scale systems. The variable porosity architecture and computational design
methodology presented here thus offers a novel pathway for automatically
generating spatially engineered electrode structures with improved power
performance.

###### keywords:

Flow-through electrodes, Redox flow batteries, Optimization, Energy storage

Victor A. Beck*, Jonathan J. Wong, Charles F. Jekel, Daniel A. Tortorelli,
Sarah E. Baker, Eric B. Duoss, _and_ Marcus A. Worsley

V. A. Beck, J. J. Wong, C. F. Jekel, S. E. Baker, E. B. Duoss, M. A. Worsley
Lawrence Livermore National Laboratory, Livermore, CA 94550
*Email Address<EMAIL_ADDRESS>

D. A. Tortorelli
University of Illinois, Urbana, IL 61801
Lawrence Livermore National Laboratory, Livermore, CA 94550

##  1 Introduction

Electricity generated from renewable sources is a continually growing
component of global energy production and a key driver for a sustainable
energy future [1, 2, 3]. Further expansion requires efficient and cost-
effective integration into existing power distribution systems but
intermittency and curtailment remain a challenge [4, 5, 6]. A number of
strategies have emerged to address these issues, including electrochemical
energy storage [7] and repurposing otherwise wasted electricity to electrify
chemical manufacturing [8, 9, 10, 11]. The direct electrochemical conversion
of CO2 is an especially powerful avenue as it simultaneously combines storage,
chemical synthesis, and carbon-removal [12, 13, 14]. As these advances are
translated from the laboratory to industrial scale, energy efficient operation
will become increasingly important to ensure economic viability [2, 7, 15, 6].

Porous flow-through electrodes are routinely used as the core reactor
component across these applications and are ubiquitously employed for
electrochemical energy storage using redox flow batteries [4, 2]. Flow battery
performance is closely tied to the porous electrode properties. The electrode
is often a disordered, homogeneous collection of micron-scale, electroactive
particles like carbon fibers, felts or spherical substrates, any of which may
be further coated with catalyst [4, 16, 17]. These materials seek to maximize
the surface reaction while minimizing overpotential and hydraulic losses.
However, the material properties to meet these requirements are inherently
adversarial and present a major challenge in attaining high performance. Open
structures are necessary to allow fluid penetration, promote mass transfer,
reduce pumping losses and supply reactants to the surface, but permeable
geometries will reduce the solid fraction and require low hydrodynamicaly
accessible surface area. In turn, the increasing porosity, decreasing
intrinsic surface area, and lower overall conductance lead to greater kinetic
and Ohmic losses [4, 17].

Previous demonstrations of high power flow batteries have circumvented these
issues by engineering the assembly to enable the use of very thin electrodes
[18, 19]. The improved performance is attributed in part to the significantly
decreased area specific resistances relative to thicker electrodes, like
uncompressed carbon felts [20, 21]. Controlling electrode thickness and
compression becomes an effective, bulk parameter to control the gross
electrode microstructure, impacting average conductance, permeability, and
mass transfer [20, 21, 22, 23].

To further drive performance, these architectures use sophisticated flowfields
to appropriately distribute reactants across the electrode surface [24, 25,
26]. This important approach partially externalizes the challenges of
balancing mass transport and electrochemical losses from the electrode to the
fluid distribution system, providing further design freedom but at the cost of
increased complexity. Previous studies have thus employed a combination of
numerical [27, 28, 29] and combined numerical and experimental [30] approaches
to develop engineering guidelines for flowfield channel dimensions and layouts
to maximize peak power and efficiency. More recent work has employed X-ray
computed tomography to simultaneously assess the impact of non-uniform
compression and flowfield arrangement, thus connecting bulk engineering
parameters to the electrode microstructure and its effective properties [31].
Indeed, a growing body of research has focused on further establishing the
connection between microstructure and hydraulic, mass transport, and
electrochemical properties [23, 31, 32, 33].

As a complement to developing new assembly architectures, engineering the
electrode structure directly is emerging as a viable route for improving
performance [15, 17]. Holes made using laser perforation were used to create
mass transport channels in carbon paper electrodes and increase peak power
[34]. Similarly, slots milled into a large scale carbon felt electrode
improved fluid distribution and decreased pumping losses without employing a
costly flowfield [35]. Dual-scale electrodes created by etching carbon papers
[36] or combining electrospun fiber mats with a backing layer [37] have
enabled even more granular engineering of the structure. Similar dual-scale
concepts have been introduced for lithium-ion electrodes [38, 39] and extended
to create continuously variable porosity electrodes [40] which have been
recently demonstrated to lead to improved rate capability while maximizing
energy density [41]. However, to date, this novel idea has not been applied to
make flow-through electrodes.

Additive and advanced manufacturing techniques can be employed to further
extend and control the structural complexity of electrode materials [42, 43,
44, 45]. Porous electrodes with superior mass transport have been created from
carbon and graphene aerogels using direct ink writing for use in
supercapacitors [42]. Flow through electrodes made from metals [43, 45],
including nickel and stainless steel [44], have been produced at varying
scales with complex, flow-controlling architectures. The resolution of the 3D
printed, flow-through electrodes leads to feature sizes that exceed those of
conventional electrodes [17] by 1-2 orders of magnitude. However, a number of
advanced manufacturing technologies exist with resolution approaching
micrometers [46, 47] and tenths of micrometers [48] using materials that are,
or can be readily transformed to be, suitable for use as electrodes.

The near arbitrary controlled provided by these techniques cannot be fully
exploited without advanced analysis and design tools to guide the electrode
architecture. Simulation has been used extensively and effectively to develop
a more detailed understanding of the transport and electrochemical processes
in flow batteries [49, 50]. The computational efforts have additionally
provided design guidance, identified key control variables, and provided
useful heuristics highlighting the importance of flow uniformity when
engineering the electrode assembly [51, 52, 53, 54].

The methodology necessarily involves beginning with an initial system
architecture, analyzing the system, and then improving it through human-driven
iteration. This process can be laborious and, crucially, explores only a
limited portion of the design phase-space. A novel, alternative technique is
to use automatic design algorithms, like topology optimization[55], to invert
the design process and aide in the phase-space exploration. Instead of
evaluating the performance of a proposed architecture, a performance target is
specified and the algorithm iterates over permissible architectures to meet
the target. This can lead to intriguing, non-intuitive designs which are
nevertheless high performance, as have been recently demonstrated for
flowfield design in flow batteries [56] and fuel cells [57]. Similar inverse
design concepts have also been used to optimize the porosity of lithium ion
batteries [40] but have never been applied to design the structure of flow
through electrodes.

We introduce the concept of algorithmically designed, microarchitected
variable porosity 3D flow through electrochemical reactors. As a specific
application, we focus on energy storage by designing flow through electrodes
for vanadium redox flow batteries. We begin by describing our modeling,
simulation and optimization methodology, including using high resolution
continuum simulation to develop a homogenized description of the constituent
microstructure unit cell. High performance computing is then employed to
determine optimal distributions of the spatially varying unit cell porosities
to maximize power efficiency across operating conditions. The resultant
architectures are evaluated for their power performance and compared to bulk,
porous flow through electrodes. The mechanisms leading to improved power
efficiency are identified and connected to the underlying electrode structure.
Finally, we demonstrate how this computational design methodology can be used
to scale-up electrodes while minimizing power efficiency losses. The design
methodology provides a framework for automatically generating high
performance, architected 3D electrodes which can fully exploit the design
freedom from advanced and additive materials manufacturing techniques.

##  2 Computational Design of Electrode Architecture

Figure 1: (a) An example of a variable porosity, microarchitected electrode.
The electrode is split into portions with unit cell with porosity above and
below 0.5 for visualization. The porosity is controlled by fixing the unit
cell length but allowing the rod diameter to vary. The unit cell size is
arbitrary and is chosen to be large in this example to aide with
visualization. (b) The microarchitected electrode can be equivalently
described with a spatially varying continuous porosity field (same legend as
(a)). (c) The design domain is the negative half-cell compartment of a
vanadium flow battery. The current collector and membrane have dimension of 2
cm x 2 cm. The inlet and outlet are square with area 0.02 cm x 0.02 cm. The
compartment overall dimension is 2.4 cm x 2 cm x 0.5 cm.

A combination of modeling, simulation and computational design optimization is
used to generate the 3D, architected porosity electrodes presented in this
work. The techniques described here are generally applicable to dilute,
single-phase, flow-through porous electrochemical reactors. Here, this
methodology is applied to determine optimal electrode architectures for
electrochemical energy storage applications. Specifically, we model the
electrode when used in the negative half-cell of a vanadium flow-through
battery and seek to optimize the electrode architecture to minimize the power
loss at fixed flow rates and discharge currents.

As shown in Figure 1, the electrode is divided into regular, isotruss unit
cells composed of rods/ fibers with equal radius. The fiber radius of these
component cells is allowed to vary spatially. Previous advanced manufacturing
work suggests scales and resolutions that would imminently permit the
manufacture of porous carbon electrodes at a scale of 100 cm2 with square unit
cells of length $L=50$ $\mu$m [42, 47, 48, 46]. Assuming equal print
resolution for both the void and the solid, the fiber radius, $r$, is
arbitrarily permitted to vary from $r_{min}=1.25$ $\mu$m to $r_{max}=11.25$
$\mu$m. The unit cell length is fixed, and thus the fiber radius can be used
to directly control the porosity and surface area per volume of the unit cell,
directly impacting the electrochemical, mass transport, and hydrodynamic
response of the electrode as described below. The resultant designed lattice
architecture can be readily transformed into a surface file and used as input
for advanced manufacturing techniques. Note that though this work employ an
isotruss geometry, these techniques are readily applicable to any unit cell
structure.

### Continuum simulation

Following previous 3D simulations [58, 27, 59, 54], the negative half-cell of
an all vanadium flow-through battery is modeled using porous electrode theory
[60, 61, 62]. Every point in the continuum represents both liquid and solid
and is characterized by the local porosity as determined by the local rod
radius of the unit cell,
$\epsilon\equiv\epsilon\left(r\left(\vec{x}\right)\right)$. Note that because
the unit cell rod radius changes with position, $r=r\left(\vec{x}\right)$, the
transport properties of the system will also be position dependent.

The electrolyte is a solution of $V^{2+}/V^{3+}$ in 1 M sulfuric acid. The
mass balance expression for the reactive species
$i\in\left\\{V^{2+},V^{3+}\right\\}$ is,

$\vec{\nabla}\cdot\left(\vec{v}c_{i}-D_{i}\vec{\nabla}c_{i}\right)=aj_{n,i},$
(1)

where $c_{i}$ is the species concentration, $D_{i}$ is the effective
diffusivity of the species in the liquid, $a$ is the specific area per volume,
and $\vec{v}$ is the superficial velocity. Additionally, a high, constant
conductivity solution is assumed and electromigration is ignored. The mass
transfer flux from the solid is the product of the mass transfer coefficient
and the difference between the surface and bulk concentrations:
$j_{n,i}=k_{m}\left(c^{s}_{i}-c_{i}\right).$ The current density from the
surface reaction ($V^{3+}+\>\>e^{-}\rightarrow\>\>V^{2+},U_{0}$) is modeled
using the Butler-Volmer expression,

$i_{n}=k_{0}\left(c^{s}_{\text{V(II)}}e^{\beta\Delta\Phi}-c^{s}_{\text{V(III)}}e^{-\beta\Delta\Phi}\right),$
(2)

where $k_{0}$ is the rate constant, $\beta\equiv F/2RT$, and we define
$\Delta\Phi\equiv\Phi_{1}-\Phi_{2}-U_{0}$ using the solid, $\Phi_{1}$, and
liquid, $\Phi_{2}$, potentials. The current density is related to the surface
species flux through Faraday’s Law, $i_{n}=Fj_{n,V^{3+}}=-Fj_{n,V^{2+}}$. The
liquid and solid potentials are modeled using Ohm’s law,

$-\vec{\nabla}\cdot\left(-\sigma\vec{\nabla}\Phi_{1}\right)=\vec{\nabla}\cdot\left(-\kappa\vec{\nabla}\Phi_{2}\right)=ai_{n},$
(3)

where $\sigma$ and $\kappa$ are the effective conductivities of the solid and
liquid, respectively. Finally, the flow field in the porous medium is
determined from the Navier-Stokes equation augmented with a Darcy drag,

$\rho\vec{v}\cdot\vec{\nabla}\vec{v}+\frac{\mu}{\alpha}\vec{v}=-\vec{\nabla}p+\mu\nabla^{2}\vec{v},$
(4)

where $\rho$ is the electrolyle density, $\mu$ is the electrolyle viscosity,
$p$ is the pressure, and $\alpha$ is the position dependent permeability.
Though it is included in the solution algorithm, the non-linear term makes a
negligible contribution since the Reynolds number is small throughout the
porous electrode, $Re\ll 1$.

The design domain in this work is the negative half-cell electrode compartment
of a discharging vanadium flow through battery as show in Figure 1c. For a
specified flow rate, $Q$, and inlet area $A_{in}$, the velocity at the inlet
is set to a uniform velocity normal to the boundary,

$\vec{v}=-\frac{Q}{A_{in}}\vec{n},$ (5)

The inlet species concentration is fixed at $c_{i}=1$ M. The top portion of
the domain is adjacent to the membrane, and a fixed current density is
specified as the boundary condition,

$-\kappa\frac{\partial\Phi_{2}}{\partial\vec{n}}=\frac{I}{A},$ (6)

where $I$ is the applied current and $A$ is the membrane area. The bottom
boundary opposite the membrane is the surface of the current collector with
boundary condition, $\Phi_{1}=0$. The domain exits to zero pressure, and all
other flow boundary conditions are wall-type boundaries with no-slip boundary
conditions for the velocity and no-flux boundary conditions for the pressure.
All other scalar boundary conditions are no-flux (i.e., homogeneous Neumann)
boundary conditions.

Table 1:  Physical Parameters Parameter | Value | Units | Ref
---|---|---|---
$D_{V^{2+},0}$ | 2.4 x $10^{-10}$ | m2/s | [63]
$D_{V^{3+},0}$ | 2.4 x $10^{-10}$ | m2/s | [63]
$\kappa_{0}$ | 40 | S/m | Estimate
$\sigma_{0}$ | $10^{4}$ | S/m | [64]
$U_{0}$ | -0.25 | V | [5]
$k_{0}$ | 1.75 x $10^{-7}$ | m/s | [62]
$T$ | 300 | K | Assumed
$\mu$ | 8.9 x $10^{-3}$ | Pa-s | Assumed

### Homogenization

Figure 2: (a) The isotruss unit cell is pictured with an illustrative planar
slice of the resolved CFD computation in the interior of the cell. (b) The
dependence of the porosity, area per volume, and permeability on $r/L$. The
blue regions are outside of the range of permitted fiber radii given by
$r_{min}$ and $r_{max}$ for the architectures designed in this work. The
permeability plot includes a comparison, in black, to the correlation from
Davies et al. for fibrous media [65]. The symbols are the results from CFD and
the blue lines are spline fits (c) The unit cell Sherwood number, $Sh$, is
presented as a function of the unit cell Péclet number, $Pe$, and the unit
cell porosity, $\epsilon$. The symbols are the results from CFD and the curves
correspond to the correlation equation, as shown.

The electrochemical, transport, and hydrodynamic responses of the unit cells
as a function of the rod radius must be determined to apply the porous
electrode model. Briefly, a 3D resolved, microscopic model of the isotruss
unit cell in Figure 1a is developed using the commercial CFD software package
Starccm+ (Siemens), a more detailed description of the procedure can be found
in the Supporting Information and from our previous work [66, 67]. The
intrinsic area per volume and porosity depend on the ratio of the internal rod
radius to the fixed unit cell edge length, $L$. The CAD package in Starccm+ is
used to calculate this relationship for several values of $r/L$. These points
are then fit using a spline to produce the curves in Figure 2b, and we
emphasize that these properties are position dependent:

$\displaystyle a$ $\displaystyle=a\left(r\left(\vec{x}\right)\right),$ (7)
$\displaystyle\epsilon$
$\displaystyle=\epsilon\left(r\left(\vec{x}\right)\right).$ (8)

For the diffusive and conductive properties, the Bruggeman relation is
employed[68, 61, 62, 60], and these properties are also position dependent,

$\displaystyle D_{i}$
$\displaystyle=D_{i,0}\left[\epsilon\left(r\left(\vec{x}\right)\right)\right]^{3/2},$
(9) $\displaystyle\kappa$
$\displaystyle=\kappa_{0}\left[\epsilon\left(r\left(\vec{x}\right)\right)\right]^{3/2},$
(10) $\displaystyle\sigma$
$\displaystyle=\sigma_{0}\left[1-\epsilon\left(r\left(\vec{x}\right)\right)\right]^{3/2},$
(11)

with $D_{i,0}$ the molecular diffusivity of species $i$, and $\kappa_{0}$ and
$\sigma_{0}$ are the conductivities of the liquid and solid, respectively. The
physical parameter values used are listed in Table 1.

The permeability is determined by meshing the void domain in the unit cell and
calculating the steady, fully-developed velocity field for an applied pressure
drop (recall $Re\ll 1$). Following Darcy’s Law, the slope of the linear
response of the superficial velocity, $\vec{v}$, against pressure is used to
determine the permeability [69],

$\alpha=\alpha\left(r\left(\vec{x}\right)\right).$ (12)

Note that for the isotruss the permeability tensor is isotropic and
characterized by a single scalar component. These values of permeability are
also fit to a spline and are compared in Figure 2b against correlations for
fibrous porous media [65, 69]. The correlations are expected to be accurate
only for dilute fiber beds (i.e., high porosity), and there is good agreement
with the calculation as the fiber radius to unit cell length ratio decreases.

The mass transport properties of the unit cell are determined by using the
convection-diffusion equation to simulate the transport and surface-
consumption of a dilute species assuming perfectly adsorbing fiber surfaces.
Because creeping flow is assumed, there is no $Re$ dependence, and the only
parameters are the species $Pe=r\left|\vec{v}\right|/D_{i,0}$ and the
porosity. Simulations are performed across $Pe$ and $\epsilon$ to calculate
the effective mass transfer coefficient,

$k_{m}\equiv\frac{\left|\vec{v}\right|}{aL}\left(1-\frac{\left<c_{out}\right>}{\left<c_{in}\right>}\right),$
(13)

where $\left<\cdot\right>$ is the axial velocity average (i.e., mixing-cup
average). The simulated values of the effective non-dimensional mass transfer
coefficient, $Sh\equiv k_{m}r/D_{0}$, is then fit to yield the correlation
shown in Figure 2c. In the dilute limit, the expression reduces to $Sh\approx
0.468Pe^{0.342}/\epsilon$ and is in nearly exact agreement with a previous
experimentally informed and frequently used correlation [70, 61, 24, 30]. Note
that because of the porosity and velocity dependence, the mass transfer
coefficient is also position dependent:

$k_{m}=k_{m}\left[\epsilon\left(r\left(\vec{x}\right)\right),Pe\left(\vec{v}\left(\vec{x}\right)\right)\right].$
(14)

### Optimization

Using the continuum, forward model for the electrochemical and homogenized
response of the system, the total power loss for any porosity distribution can
be calculated. The power loss objective function is defined as the sum of the
electric power losses, $P_{elec}$, and the hydraulic power losses, $P_{flow}$:

$P_{tot}=P_{elec}+P_{flow},$ (15)

with,

$P_{elec}=\int_{mem}\eta\;\vec{i}_{2}\cdot\vec{n}\;d\Gamma,$ (16)

an integral over the electrode-membrane interface, and

$P_{flow}=\int_{in}\frac{p\;\vec{v}\cdot\vec{n}}{\Psi_{pump}}\;d\Gamma,$ (17)

an integral over the inlet [30]. The current density in the liquid is given by
$\vec{i}_{2}$. The local overpotential is defined as,
$\eta=\Phi_{1,cc}-\Phi_{2,mem}-U_{0}$, where $\Phi_{1,cc}\equiv 0$ is the
potential at the current collector, $\Phi_{2,mem}$ is the potential at the
interface between the electrode and the membrane, and the Nernst potential
only contributes the standard potential since the inlet concentrations of
oxidant and reductant are set equal. The power efficiency is thus defined as
$\Xi=1-P_{tot}/IU_{0}$. In this work the pump efficiency is idealized and
assumed to be $\Psi_{pump}=1$. Lower pump efficiency is equivalent to an
increased weighting of the flow contribution to the objective function as seen
in Equation (17).

We seek to determine the distribution of unit cells with rod radius,
$r\left(\vec{x}\right)$, which will minimize the total power loss, $P_{tot}$:

$\displaystyle\min_{r\left(\vec{x}\right)}\>\>P_{tot}$ (18) $\displaystyle
s.t.~{}~{}~{}r_{min}\leq r\left(\vec{x}\right)\leq r_{max}.$

The total derivative, $dP_{tot}/dr\left(\vec{x}\right)$, subject to the
constraints imposed by the physical model in Equation (1)-(4) is calculated
using standard techniques from PDE constrained optimization (see Supporting
Information). Briefly, the continuous adjoint approach is employed to derive
analytical expressions for the adjoint PDEs [71]. This results in one adjoint
PDE per forward model PDE. For a given solution of the forward model, the
adjoint PDEs are numerically solved. The total derivitative (i.e.,
sensitivity) is then computed from the forward solution, adjoint solution, and
partial derivative of the Lagrangian function with respect to the design
variables. The porosity distribution is updated using the Method of Moving
Asymptotes (MMA) [72]. Iteration continues until the average relative change
in the objective function over the last 5 steps is less than $10^{-4}$, to
arrive at a local optimum for the rod distribution. A Helmholtz filter is used
to regularize the optimization problem and control the smoothness length scale
of the porosity variation [73]. The length scale parameter is set to $200$
$\mu$m.

The forward simulation, adjoint calculation, and gradient descent update are
all performed using bespoke code written using OpenFOAM. The domain is meshed
using $\approx 1.2M$ cubic, finite-volume cells.

##  3 Results and Discussion

The ultimate power efficiency of a flow through electrode is engineered by
balancing the losses arising from insufficient mass transport to the reactive
surfaces against the hydraulic power necessary to drive fluid to those
surfaces and provide charge conduction pathways. Below we characterize the
power losses in flow through electrodes composed of isotrusses and apply
design optimization to a three-dimensional model of coupled fluid flow,
species transport, and current distribution, to address this engineering
challenge.

### Characterizing power losses in flow through electrodes

Figure 3: (a) The total power loss is separated into total electric and flow
components for a bulk electrode, characterized by a single porosity, operated
at 40 mL/min and 400 mA/cm2. The total electric power loss is further
separated into the portion due to concentration polarization and the portion
due to internal electric resistance. The location of the minimum from each
contribution is marked. (b) The power efficiency at a fixed current density of
400 mA/cm2 as the porosity and flow rate are varied. The solid lines show the
power efficiency for a bulk, monolithic electrode. The black curve is the loci
of power efficiency maxima. The arrows show the increase in efficiency for the
optimized electrode relative to the best performing bulk electrode. (c) The
total power loss, separated into hydraulic losses, internal electric losses,
and electric transport losses is presented for both the bulk and optimized
electrodes at fixed current density of 400 mA/cm2. The histogram couplets
correspond on the left to the power losses for the bulk electrode at its
optimal porosity (i.e., the maxima in (b)) and the right to the optimized,
architected, variable porosity electrode.

The porosity dependence of the total power loss in a bulk, monolothic
electrode used as the negative half-cell of a vanadium flow-through battery
(i.e, a homogeneous electrode which can be described by a single value for the
porosity) is presented in Figure 3a. At a fixed input flow rate of 40 mL/min
and operating current density of 400 mA/cm2 the minimum power loss occurs when
using an electrode with bulk porosity $\epsilon=0.815$. Equivalently, the
electrode is a uniform set of lattices with rod radius $r=3.7$ $\mu$m. The
total power loss is decomposed via Equation (15)-(17) into contributions from
hydraulic losses and electric losses to reveal that this minimum results from
a trade-off between each component. The flow losses decrease as the porosity
of the electrode increases, with a minimum hydraulic power loss occurring at
the highest attainable porosity. This is as expected, as the more open
structure will have higher permeability. Alternatively, $\epsilon\rightarrow
1$ represents a singular limit for the electric losses since both Ohmic and
kinetic overpotential losses grow unbounded as the solid material disappears.
Instead, the electric losses generally decrease with increasing solid fraction
since this increases reactive area and the effective solid conductivity,
lowering Ohmic and surface overpotential losses.

The electric losses pass through a minimum at $\epsilon=0.711$ and eventually
begin to increase as the porosity further decreases. Increasing the solid
fraction hinders transport to the reactive surface by displacing fluid,
changing the local mass transfer coefficients, and decreasing the liquid
conductivity. To better differentiate these effects, the internal electric
loss, $P_{int}$, is determined by simulating an electrode wherein we prescribe
the surface concentrations of all species to be equal to the inlet
concentrations, reducing the problem to the solution of Equation (3). This
idealization removes all concentration variations and hence all concentration
polarization losses, and it is conceptually equivalent to canonical work on
one-dimensional porous electrode models [74]. The concentration polarization
power loss is thus defined as,

$P_{trans}\equiv P_{elec}-P_{int}.$ (19)

The internal electric loss has contributions only from Ohmic and kinetic
overpotential losses [60, 74], while the concentration losses will include
contributions from variations in both the concentration _and_ the mass
transfer coefficients. We thus equivalently refer to the concentration
polarization losses as electric losses due to insufficient species transport,
or simply electric transport losses. The internal power loss plotted in Figure
3a is the dominant contribution to the total electric power loss and the
dominant loss in general. At this flow rate and current, the transport losses
are generally too small to impact the optimal porosity, but it is noted that
these losses increase as the void volume is reduced.

### Minimizing power losses in flow through electrodes

The total power efficiency, $\Xi$, in bulk electrodes operating at flow rates
ranging from 2-100 mL/min is shown in 3b. All of the power efficiency curves
are similar and have a single optimal porosity that balances the flow and
electric power losses. A spline fit to the loci of the curve maxima yields the
expected optimal porosity as the flow rate is varied. At 400 mA/cm2 the most
efficient bulk electrode operates at the maximum of this spline curve.
Equivalently, the minimum power loss is attained when operating at 15 mL/min
with a bulk electrode porosity of $\epsilon=0.761$.

The key contribution in this work arises from allowing the geometry to vary
spatially throughout the electrode by changing and optimizing the rod radius
in each fixed unit cell. The optimization approach discussed above is used to
determine a distribution of unit cell rod radii, $r\left(\vec{x}\right)$,
which minimizes the total power loss. At each flow rate and current density,
the initial input to the optimization process is the corresponding minimum
bulk porosity. As shown in Figure 3b, in each case the optimization algorithm
is able to find a distribution of radii which leads to further decreases in
the total power loss and hence greater efficiency than the bulk electrode.
Indeed, the highest efficiency observed in Figure 3b is $\Xi=0.614$ for the
optimized electrode operating at 15 mL/min, an efficiency increase over the
bulk electrode of 13.5%. Across flow rates, the efficiency relative to the
best performing bulk electrodes is improved by as much as 310%. This
optimization procedure and analysis is repeated at current densities of 100
mA/cm2 and 200 mA/cm2 as presented in the Supporting Information, and, in
general, lower current operation leads to even higher efficiency.

Figure 4: The architected porosity electrodes resulting from the optimization
procedure at a fixed current of 400 mA/cm2 and a flow rate of (a) 2 mL/min,
(b) 5 mL/min, (c) 15 mL/min, (d) 40 mL/min, and (e) 100 mL/min. To aid in
visualization, the electrodes are split into a “solid” half with $\epsilon$
$\leq$ $0.5$, presented on the top row, and a “void” half where $\epsilon$ $>$
$0.5$ presented on the bottom row. The electrodes are colored by height, with
red adjacent to the membrane and blue adjacent to the current collector.

The resultant locally optimal, variable porosity electrode geometries are
visualized by splitting the electrode into two volumes. The “solid” portion of
the electrode lumps all unit cells with $\epsilon\leq 0.5$ while the “void”
portion of the electrode lumps all unit cell with $\epsilon>0.5$. The former
is more closed to fluid flow and will behave more like a pure solid (i.e.,
preventing flow and enhancing electrical conduction), while the latter is more
open and will behave more like a void (i.e., permitting flow and reducing
electrical conduction). The two portions of the designed electrode are
presented in Figure 4, and it is emphasized that though we are visualizing
solid blocks, all regions in the electrode are composed of unit cells (see
Figure 1a).

Figure 5: Planar slices of the electrode porosity parallel to the current
collector and membrane are presented for the architected electrode designed at
flow rate of (a) 2 mL/min, (b) 5 mL/min, (c) 15 mL/min, (d) 40 mL/min, and (e)
100 mL/min. The planes are created at 1.5 mm, 2.5 mm, 3.5 mm, and 4.5 mm from
the current collector. The velocity vector field is also presented, with
vector length proportional to relative velocity magnitude within each column.

Representative slices of the porosity distribution and the flow vector field
in planes parallel to the current collector are presented in Figure 5. At low
flow rates, the optimal electrode geometry tends towards a large, mostly solid
block with a thin high-flow region near the membrane. While at higher flow
rates several large channels are carved into the solid block, and these grow
larger with increasing flow rate. Across the designed electrodes, the porosity
increases with distance from the current collector and all electrodes show a
thin, high porosity, distributed flow region adjacent to the membrane.

The highest efficiency, optimized electrode is pictured in Figure 4c and
Figure 5c. It is a mostly solid structure near the current collector with
several low porosity channels tunneling through the block that direct fluid
from the current collector to the membrane. Near the membrane the number of
channels increases, eventually leading to the high porosity region and
appearing to spread the fluid along the membrane surface. In contrast, the
electrode operated at 5 mL/min (see Figure 4a and Figure 5a) begins by
distributing the majority of the fluid along the top, upstream edge of the
electrode. The channels in the highest effiency electrode thus appear to
increase the number of injection points for the fluid and distribute them
throughout the volume. In essence, the variable porosity electrode is
performing an analogous function to a flowfield plate [24, 25, 26, 27, 28, 29,
30], but now the flow control is integral to the electrode. Similar flow-
directing structures are seen in Figure 4d and Figure 5d, but, interestingly,
these are most pronounced in the best performing, architected electrode.

### Comparing power losses in bulk electrodes to variable porosity electrodes

The contributions to the total power loss for the bulk electrodes are
presented in the first bar of each pair in Figure 3c. At low flow rates, the
electric transport losses are major contributors to the total loss in the bulk
electrodes. There is insufficient mass transfer at low flow rates.
Alternatively, at high flow rates the electric transport losses are small but
at the cost of increased hydraulic losses. At each of these extrema the
performance of the bulk electrode is maximized by using a more open porosity.
At low flow rates, low porosity is necessary to enable greater mass transport,
while at high flow rates it is instead the reduced flow resistance from higher
permeability which is necessary to improve performance. Unfortunately, the
higher porosity also leads to lower overall conductivity and an increase in
internal electric losses. The loci of maxima in Figure 3b can thus be
understood as an optimization over operational flow rate to find the lowest
porosity value which minimizes losses associated with fluid distribution:
Operation at 15 mL/min _enables_ the use of the more solid, conductive
electrode.

The architected electrodes are compared to the best performing bulk electrodes
(i.e., the electrodes corresponding to the minima of the solid curves in
Figure 3b) operating at the same current and flow rate. In all cases, the
variable porosity, optimized electrodes lead to a reduction in all
contributions to the power loss. At the low flow rates the electric transport
losses due to insufficient species flow are minimized by forcing all of the
flow into a thin region adjacent to the membrane, effectively creating a thin,
porous electrode (i.e., Figure 4a). For high flow rates, an excess of flow
amply supplies the reaction near the membrane but leads to large hydraulic
losses which are minimized by carving by-passes in the thick bulk region
(i.e., Figure 4e). The internal resistance of the electrodes is also lowered
by using the architected porosities. However, unlike for a bulk electrodes
above, architected porosity relaxes the constraint between porosity and
conductivity observed above for the bulk electrode allowing for increased
power efficiency. Interestingly, the internal losses for the bulk electrode
all appear to reach similar limiting values around 0.15 W. These are all lower
than the internal losses observed in the bulk electrode and imply that
variable porosity can improve performance even when mass transport losses are
unimportant. Similar trends are observed at lower operating current densities
(see Supporting Information).

Figure 6: Averages of the (a) porosity, (b) area per volume, (c) planar
component of the velocity, $u_{x,z}$, (d) reaction current made dimensionless
with the applied current divided by the electrode volume, $V$, and (e) the
concentration of the reductant in planes parallel to the current collector and
located at the given height above the current collector. In (f) the one-way
flow rate defined in Equation (21) is presented. The solid lines correspond to
the variable porosity electrodes while the dashed lines correspond to the bulk
electrodes at the optimal porosity (i.e., maxima of curves in Figure 3b. The
current density is fixed at 400 mA/cm2.

### Characterizing the porosity distribution

The porosity of the optimized electrodes is averaged on planes directly above
the current collector as presented in Figure 6a. At low flow rates, the planar
averages show that the region near the current collector is near the lowest
available value of the porosity, confirming that in Figure 4a-b the “solid”
regions are indeed monolithic and as solid as possible. The optimization
algorithm has effectively created a conductive block extending from the
current collector with a thickness of approximately 3 mm stacked below a 2 mm
porous electrode with gradually increasing porosity. This is reminiscent of
commonly recognized heuristics for engineering high power fuel cells and flow
batteries which employ thin geometries with high porosity papers or felts [18,
19, 75]. An important departure, however, is that in previous work the
geometry is a thin flow channel with a constant porosity electrode. The
variable porosity electrodes show that even for the thin, membrane adjacent 2
mm region, a smooth, gradual increase in porosity leads to greater power
efficiency since the optimization procedure would adjust the porosity to be
uniform near the electrode if this led to continued improvement.

The gradual increase in porosity as the membrane is approached is observed
across all of the optimized electrodes, and the transition occurs at
approximately the same position, 3 mm. Further, the ultimate value of the
porosity at the membrane is fairly constant across all architected electrodes,
$\epsilon\approx 0.84$. This value is within the range of optimal porosities
determined for the bulk electrodes but exceeds the porosity of the best
performing bulk electrode (e.g., the location of the maxima in Figure 3b). The
physical parameters of the system suggest that in the absence of concentration
polarization a characteristic length scale for the bulk electrode is the
penetration depth [60],

$L_{p}=\sqrt{\frac{RT}{ai_{0}F\left(\kappa^{-1}+\sigma^{-1}\right)}}$ (20)

which will vary with porosity, but for the range of porosities and intrinsic
areas studied here is $\approx 0.1-0.9$mm. It is not obvious how to assign an
equivalent characteristic length for the variable porosity electrodes, but the
values are of the same magnitude as the characteristic length for the porosity
variation in Figure 6a.

### Elucidating mechanisms for enhanced performance

The architected electrodes lead to higher power efficiencies across applied
current densities and flow rates (see Supporting Information for additional
results). The resultant flow, concentration and current fields for electrodes
operating at 400 mA/cm2 are analyzed and compared to the bulk electrode fields
to understand the factors that lead to these improvements. These quantities
are again averaged on planes parallel and directly above the current collector
and presented in Figure 6.

The planar averages of the reaction current density, the magnitude of the
righthand side of Equation 3, normalized by the volume averaged current,
$I/V$, are shown in Figure 6d, revealing the spatial distribution of the
electrochemical reaction and thus the electrode utilization. Driving this
reaction at minimal power loss is the main goal of the electrochemical
engineering problem posed in this work. Recall that for the bulk electrodes
the highest power efficiency was attained for the bulk electrode operating at
a flow rate of 15 mL/min. When compared to the other bulk electrodes, it is
evident from the current distribution that in this electrode the reaction rate
near the current collector is lowest and, because the operating current is
fixed, the reaction rate near the membrane is greater. All of the architected
electrodes exceed the power efficiency of this bulk electrode and additionally
lead to lower internal losses. These electrodes all show further decreased
reaction rate in the bottom portion (i.e., the region $<$3 mm from current
collector), and consequently greater reaction rate near the membrane, relative
to the bulk electrodes. High performance electrodes drive reaction closer to
the membrane, and all of the optimized electrodes show exponential growth of
the reaction rate as the membrane is approached. In the electrodes studied
here, improved volumetric utilization of the electrode is not a requisite for
improved power efficiency.

To maintain high reaction currents near the membrane, the reductant needs to
be supplied at high molar fluxes, as quantified by the product of the local
velocity and the concentration. The magnitude of the velocity projected onto
the averaging plane is presented in Figure 6b. For the bulk electrodes the
nearly constant average value shows that the flow is approximately evenly
distributed with distance from the current collector. All of the curves are
similar and scale with the applied flow rate since the permeability is
constant and nearly identical across the bulk electrodes and the flow
equations are linear at low $Re$. Because at low $Re$ the flow vector field
orientation does not change if the porosity is uniform, the only avenue for
increasing reductant concentration near the membrane is to increase the flow
rate. Indeed, in Figure 6c it is only at the highest flow rates that a
uniform, high concentration is attained for the bulk electrodes. At lower flow
rates there is a strong depletion of the active species and thus a strong
increase in the electric transport losses (see Figure 3c).

In contrast, the optimized electrodes succeed in attaining higher
concentrations of the reductant through improved flow management without
incurring unacceptably large hydraulic losses. At low flow rates (2-15
mL/min), the planar averages in Figure 6c show that the variable porosity
electrodes drive most of the fluid against the membrane and, with the
exception of the electrode designed at 100 mL/min, exceed the flow rate of the
bulk electrode operating at the same flow rate. The large velocities seen near
the current collector at the highest flow rates are due to the large bypass
channels seen in Figure 4d-e.

The impact of the flow distribution in the optimized electrodes on the
concentration is readily evident in Figure 6e. The concentrations near the
membrane of these architected electrodes are higher than their bulk
counterparts, with the exception of the optimization at 100 mL/min where the
optimized electrodes has led to a decrease in the concentration near the
membrane. Consequently, the electric transport losses seen in Figure Figure 3c
are reduced relative to the bulk electrode at all flow rates except 100
mL/min, where there is oversupply of species.

For planes parallel to the current collector, the vertical transport induced
by the electrode can be determined from the flow rate crossing the plane in
only one-direction:

$\left<Q_{y}|Q_{y}>0\right>_{x,z}\equiv\frac{1}{2}\int_{x,z}\left|\vec{v}\cdot\vec{e}_{y}\right|+\vec{v}\cdot\vec{e}_{y}\;dA,$
(21)

where $\vec{e}_{y}$ is a unit vector normal to the current collector and
pointing in the direction of the membrane. It is equivalent to an integral of
the vertical component of the vectors in Figure 5 and provides a metric to
quantify the vertical exchange of material. Indeed, when plotted in Figure 6f
the average one-way flow rate, $\left<Q_{y}|Q_{y}>0\right>$, shows that the
optimized electrodes all lead to greater flow towards the membrane relative to
the bulk electrode, especially near the membrane. In contrast, for the bulk
electrode the only mechanism to increase vertical transport is to increase the
flow rate, but this leads to the wasteful increase of hydraulic losses
throughout the electrode. Interestingly, as seen from the 40 mL/min and 100
mL/min variable porosity electrode curves, the flow of material to the
membrane is not maximized, instead it is optimized toward an apparent limiting
value that is nearly equivalent to the vertical transport seen in the 15
mL/min variable porosity electrode. Any supply in excess leads to unnecessary
hydraulic losses (see Figure 3c).

Figure 7: The probability distribution function of (a) the magnitude of the
velocity and (b) the perpendicular component of velocity, $v_{y}$, as a
function of the applied flow rate. The velocities are normalized by the area
specific velocity, the ratio of the flow rate to the membrane area. The
probability distribution functions for the optimal, bulk electrodes are
presented in the black curves. Note that due to the linearity of the flow
equations, the normalized flow field of the bulk electrodes is identical
across flow rates and the probability distribution function for the bulk
electrodes can be represented by a single curve labeled, “Bulk.”

The central role of increased vertical transport is further supported by
analyzing the probability distribution function of the vertical component of
the velocity normalized by the area specific velocity, $Q/A$, as shown in
Figure 7a. For all of the designed electrodes the velocity distribution is
much broader, and this is especially pronounced for the electrode optimized at
15 mL/min. The 15mL/min optimized electrodes shows the most induced vertical
flow. As shown in Figure 7b, the variable porosity electrodes lead to a
broader distribution of high velocities relative to the bulk electrodes. In
addition to enhancing convective transport, this serves to increase the local
reaction rates by increasing the flow speed dependent mass transfer
coefficients.

From these observations we can thus conjecture an interpretation of the
resultant electrode architectures in Figure 4. The optimization procedure is
not only sculpting the gross features of the electrode, creating a thin, high
porosity electrode near the membrane, it is also developing an integral flow
distribution system, essentially a porous flowfield, to manage the flow paths.
This is especially evident in the “braided” design seen in Figure 4c, where we
can hypothesize that the architectural features serve to distribute fluid
vertically as well as to allow for maximal flow distribution, with minimal
hydraulic loss, along the electrode. Of course, this is a single
interpretation focused on flow distribution: The optimization procedure is
general and algorithmic, and simultaneously optimizes flow paths, convective
paths, and conduction paths to control reaction rates and automatically arrive
at an optimal architecture irrespective of analogies to other engineering
approaches.

Figure 8: (a) The architected porosity distribution for a scaled-up electrode
designed for a half cell with dimensions 8 cm x 8 cm x 5 mm. (b) The total
power loss, separated into hydraulic losses, internal electric losses, and
electric transport losses for the best performing, scaled bulk electrode
(first bar) and the optimized electrode (second bar).

### Optimization of large scale electrodes

The electrodes analyzed above are of typical experimental bench scale and have
an electrode area of 4 cm2. However, though these dimensions are expedient for
analysis, it is often unclear how to translate performance insights realized
at these smaller scales to the larger scales needed for industrial
application. Industrially useful reactors are often several orders of
magnitude larger [35, 75]. Thus, as a final demonstration of the utility of
the approach described here, the planar dimensions of the half-cell
compartment are scaled by 16x to a size of 8cm x 8cm x 5mm (i.e., electrode
area of 64 cm2) while keeping the current density and area specific velocity
fixed. The mesh resolution is constant requiring the solution of an
optimization problem with $\approx$20M design variables. The best performing
bulk electrode at 4 cm2 is scaled to this new dimension and compared to an
architected electrode designed at this new dimension as shown in Figure 8a.

When the bulk electrode is scaled in the absence of fluid distribution systems
the power efficiency drops from $\Xi=0.541$ to $\Xi=0.384$, a 29.0% reduction.
Alternatively, the architected, optimized electrode experiences a much smaller
reduction in power efficiency dropping from $\Xi=0.614$ to $\Xi=0.539$, only a
12.3% reduction and thus retaining nearly all of the power efficiency upon
scale-up. Equally important, the scaled-up architected electrode exceeds the
power efficiency of the scaled-up bulk electrode by 40.3%.

The enhanced performance is further reflected in the power losses presented in
Figure 8b. The structures in the optimized electrode again lead to lower
contributions to all loss components. A bulk electrode would, in practice,
also use a fluid distribution system, but it is not possible for an external
manifold to arbitrarily control the detailed local flow in the interior of the
electrode. Indeed, even if it were possible the target flow distribution is
generally unknown and would need to be determined via laborious and time
inefficient testing and iteration. In this context, the advantages of the
architected electrodes described here are self-evident, the variable porosity
leads to expanded design flexibility that can be readily harnessed by the
computational design algorithm to automatically generate structures with
significantly enhanced performance.

##  4 Conclusions

This work has introduced the concept of algorithmically designed,
microarchitected 3D flow through electrochemical reactors. The technique was
used to design electrodes for the negative half-cell of a vanadium redox flow
battery. Across flow rates, utilizing architected electrodes led to power
efficiency gains of 13.5% - 310% over uniform porosity, bulk electrodes.
Decomposing the power losses into flow losses, internal electric losses, and
concentration polarization losses revealed that the algorithmically designed
electrodes manipulate the porosity distribution to reduce the dominant loss
contributors. These architectural changes included gross features like
creating thin, high conductivity geometries at low flow rates, and creating
large, low hydraulic resistance bypass channels for operation at higher flow
rates. However, the best performing electrodes also included several more
subtle architectural changes that shape both fluid and conductive pathways,
creating an integrated fluid distribution system without sacrificing
electrochemical performance. These conclusions were supported through analysis
of the local reaction rates, flow, and concentrations, where it was
demonstrated that the designed, variable porosity leads to improved transport
from the current collector toward the membrane, greater reaction adjacent to
the membrane, and more vertical transport between the membrane and current
collector. Algorithmically designed, architected porosity electrodes thus
provide a novel pathway for electrode engineering without necessitating
complex external manifolds like flowfields. As a final demonstration of the
utility of this technique, the framework was employed to scale-up to 64 cm2
from 4 cm2. The scaled electrode incurred only a 12.3% reduction in power
efficiency and exceeded the power efficiency of the best performing bulk
electrodes by 40.3%. This work has thus provided a versatile, novel, variable
porosity electrode architecture with vastly expanded design flexibility and a
computational methodology to automatically generate high performance
structures across length scales. The techniques described here are broadly
applicable across flow-through electrodes and liquid-fed porous reactors.

Supporting Information

Supporting Information is available from the author.

Acknowledgements

This work was performed under the auspices of the U.S. Department of Energy by
Lawrence Livermore National Laboratory under Contract DE-AC52-07-NA27344 and
was supported by the LLNL-LDRD program under project numbers 16-ERD-051,
19-SI-005, and 19-ERD-035. LLNL Release Number LLNL-JRNL-819457.

## References

  * [1] S. Chu, A. Majumdar,  _Nature_ 2012, _488_ , 7411 294.
  * [2] G. L. Soloveichik,  _Chem. Rev._ 2015, _115_ , 20 11533.
  * [3] S. Chu, Y. Cui, N. Liu,  _Nat. Mater._ 2017, 16–33.
  * [4] A. Z. Weber, M. M. Mench, J. P. Meyers, P. N. Ross, J. T. Gostick, Q. Liu,  _J Appl Electrochem_ 2011, _41_ , 10 1137.
  * [5] L. Li, S. Kim, W. Wang, M. Vijayakumar, Z. Nie, B. Chen, J. Zhang, G. Xia, J. Hu, G. Graff, J. Liu, Z. Yang,  _Adv. Energy Mater._ 2011, _1_ , 3 394.
  * [6] L. F. Arenas, C. Ponce de León, F. C. Walsh,  _Current Opinion in Electrochemistry_ 2019, _16_ 117\.
  * [7] T. M. Gür,  _Energy Environ. Sci._ 2018, _11_ , 10 2696.
  * [8] Z. J. Schiffer, K. Manthiram,  _Joule_ 2017, _1_ , 1 10.
  * [9] M. J. Kim, Y. Seo, M. A. Cruz, B. J. Wiley,  _ACS Nano_ 2019, _13_ , 6 6998.
  * [10] A. Shatskiy, H. Lundberg, M. D. Kärkäs,  _ChemElectroChem_ 2019, _6_ , 16 4067.
  * [11] T. Noël, Y. Cao, G. Laudadio,  _Acc. Chem. Res._ 2019, _52_ , 10 2858.
  * [12] J. W. Ager, A. A. Lapkin,  _Science_ 2018, _360_ , 6390 707.
  * [13] C.-T. Dinh, T. Burdyny, M. G. Kibria, A. Seifitokaldani, C. M. Gabardo, F. P. G. de Arquer, A. Kiani, J. P. Edwards, P. De Luna, O. S. Bushuyev, C. Zou, R. Quintero-Bermudez, Y. Pang, D. Sinton, E. H. Sargent,  _Science_ 2018, _360_ , 6390 783.
  * [14] V. Vedharathinam, Z. Qi, C. Horwood, B. Bourcier, M. Stadermann, J. Biener, M. Biener,  _ACS Catal._ 2019, _9_ , 12 10605.
  * [15] M. L. Perry, A. Z. Weber,  _Journal of The Electrochemical Society_ 2016, _163_ , 1 A5064.
  * [16] M. Park, J. Ryu, W. Wang, J. Cho,  _Nat. Rev. Mater._ 2017, _2_ 16080\.
  * [17] A. Forner-Cuenca, F. R. Brushett,  _Current Opinion in Electrochemistry_ 2019, _18_ 113\.
  * [18] D. S. Aaron, Q. Liu, Z. Tang, G. M. Grim, A. B. Papandrew, A. Turhan, T. A. Zawodzinski, M. M. Mench,  _Journal of Power Sources_ 2012, _206_ 450\.
  * [19] Q. H. Liu, G. M. Grim, A. B. Papandrew, A. Turhan, T. A. Zawodzinski, M. M. Mench,  _Journal of The Electrochemical Society_ 2012, _159_ , 8 A1246.
  * [20] L. D. Brown, T. P. Neville, R. Jervis, T. J. Mason, P. R. Shearing, D. J. L. Brett,  _Journal of Energy Storage_ 2016, _8_ 91\.
  * [21] T. J. Davies, J. J. Tummino,  _C_ 2018, _4_ , 1 8.
  * [22] R. Banerjee, N. Bevilacqua, A. Mohseninia, B. Wiedemann, F. Wilhelm, J. Scholta, R. Zeis,  _Journal of Energy Storage_ 2019, _26_ 100997\.
  * [23] A. A. Wong, M. J. Aziz,  _Journal of The Electrochemical Society_ 2020, _167_ , 11 110542.
  * [24] R. M. Darling, M. L. Perry,  _Journal of The Electrochemical Society_ 2014, _161_ , 9 A1381.
  * [25] T. J. Latha, S. Jayanti,  _Journal of Power Sources_ 2014, _248_ 140\.
  * [26] J. Houser, A. Pezeshki, J. T. Clement, D. Aaron, M. M. Mench,  _Journal of Power Sources_ 2017, _351_ 96\.
  * [27] Q. Xu, T. S. Zhao, P. K. Leung,  _Applied Energy_ 2013, _105_ , C 47.
  * [28] E. Knudsen, P. Albertus, K. T. Cho, A. Z. Weber, A. Kojic,  _Journal of Power Sources_ 2015, _299_ 617\.
  * [29] X. Ke, J. M. Prahl, J. I. D. Alexander, R. F. Savinell,  _Journal of Power Sources_ 2018, _384_ 295\.
  * [30] M. R. Gerhardt, A. A. Wong, M. J. Aziz,  _Journal of The Electrochemical Society_ 2018, _165_ , 11 A2625.
  * [31] N. Kulkarni, M. D. R. Kok, R. Jervis, F. Iacoviello, Q. Meyer, P. R. Shearing, D. J. L. Brett,  _Journal of Power Sources_ 2019, _426_ 97\.
  * [32] M. D. R. Kok, R. Jervis, T. G. Trantet, M. A. Sadeghi, D. J. L. Brett, P. R. Shearing, J. T. Gostick,  _Chemical Engineering Science_ 2019, _196_ 104\.
  * [33] A. Forner-Cuenca, E. E. Penn, A. M. Oliveira, F. R. Brushett,  _Journal of The Electrochemical Society_ 2019, _166_ , 10 A2230.
  * [34] I. Mayrhuber, C. R. Dennison, V. Kalra, E. C. Kumbur,  _Journal of Power Sources_ 2014, _260_ , C 251.
  * [35] A. Bhattarai, N. Wai, R. Schweiss, A. Whitehead, G. G. Scherer, P. C. Ghimire, T. M. Lim, H. H. Hng,  _Applied Energy_ 2019, _236_ 437\.
  * [36] X. L. Zhou, Y. K. Zeng, X. B. Zhu, L. Wei, T. S. Zhao,  _Journal of Power Sources_ 2016, _325_ 329\.
  * [37] Q. Wu, Y. Lv, L. Lin, X. Zhang, Y. Liu, X. Zhou,  _Journal of Power Sources_ 2019, _410-411_ 152\.
  * [38] C. L. Cobb, M. Blanco,  _Journal of Power Sources_ 2014, _249_ 357\.
  * [39] V. P. Nemani, S. J. Harris, K. C. Smith,  _Journal of The Electrochemical Society_ 2015, _162_ , 8 A1415.
  * [40] S. Golmon, K. Maute, M. L. Dunn,  _Journal of Power Sources_ 2014, _253_ 239\.
  * [41] X. Lu, A. Bertei, D. P. Finegan, C. Tan, S. R. Daemi, J. S. Weaving, K. B. O’Regan, T. M. M. Heenan, G. Hinds, E. Kendrick, D. J. L. Brett, P. R. Shearing,  _Nat Commun_ 2020, _11_ 2079\.
  * [42] C. Zhu, T. Liu, F. Qian, T. Y.-J. Han, E. B. Duoss, J. D. Kuntz, C. M. Spadaccini, M. A. Worsley, Y. Li,  _Nano Lett._ 2016, _16_ , 6 3448.
  * [43] J. Lölsberg, O. Starck, S. Stiefel, J. Hereijgers, T. Breugelmans, M. Wessling,  _ChemElectroChem_ 2017, _4_ , 12 3309.
  * [44] L. F. Arenas, C. Ponce de León, F. C. Walsh,  _Electrochemistry Communications_ 2017, _77_ 133\.
  * [45] J. Hereijgers, J. Schalck, J. Lölsberg, M. Wessling, T. Breugelmans,  _ChemElectroChem_ 2019, _6_ , 2 378.
  * [46] B. D. Moran,  Large area projection micro stereolithography, 2016,  US Patent App. 14/688,187.
  * [47] R. M. Hensleigh, H. Cui, J. S. Oakdale, J. C. Ye, P. G. Campbell, E. B. Duoss, C. M. Spadaccini, X. Zheng, M. A. Worsley,  _Mater. Horiz._ 2018, _5_ , 6 1035.
  * [48] S. K. Saha, D. Wang, V. H. Nguyen, Y. Chang, J. S. Oakdale, S.-C. Chen,  _Science_ 2019, _366_ , 6461 105.
  * [49] O. C. Esan, X. Shi, Z. Pan, X. Huo, L. An, T. S. Zhao,  _Adv. Energy Mater._ 2020, _22_ 2000758\.
  * [50] B. K. Chakrabarti, E. Kalamaras, A. K. Singh, A. Bertei, J. Rubio-Garcia, V. Yufit, K. M. Tenny, B. Wu, F. Tariq, Y. S. Hajimolana, N. P. Brandon, C. T. John Low, E. P. L. Roberts, Y.-M. Chiang, F. R. Brushett,  _Sustainable Energy Fuels_ 2020, _4_ , 11 5433.
  * [51] J. Escudero-González, P. A. López-Jiménez,  _Journal of Power Sources_ 2014, _251_ , c 243.
  * [52] T. Jyothi Latha, S. Jayanti,  _J Appl Electrochem_ 2014, _44_ , 9 995.
  * [53] Q. Zheng, F. Xing, X. Li, G. Ning, H. Zhang,  _Journal of Power Sources_ 2016, _324_ 402\.
  * [54] C. Yuan, F. Xing, Q. Zheng, H. Zhang, X. Li, X. Ma,  _Energy_ 2020, _194_ 116839\.
  * [55] O. Sigmund, K. Maute,  _Struct Multidisc Optim_ 2013, _48_ , 6 1031.
  * [56] C.-H. Chen, K. Yaji, S. Yamasaki, S. Tsushima, K. Fujita,  _Journal of Energy Storage_ 2019, _26_ 100990\.
  * [57] R. Behrou, A. Pizzolato, A. Forner-Cuenca,  _Int. J. Heat Mass Transfer_ 2019, _135_ 72\.
  * [58] X. Ma, H. Zhang, F. Xing,  _Electrochimica Acta_ 2011, _58_ 238\.
  * [59] Y. Wang, S. C. Cho,  _Journal of The Electrochemical Society_ 2014, _161_ , 9 A1200.
  * [60] J. Newman, W. Tiedemann,  _AIChE Journal_ 1975, _21_ , 1 25.
  * [61] J. A. Trainham, J. Newman,  _Journal of The Electrochemical Society_ 1977, _124_ , 10 1528.
  * [62] A. A. Shah, M. J. Watt-Smith, F. C. Walsh,  _Electrochimica Acta_ 2008, _53_ , 27 8087.
  * [63] T. Yamamura, N. Watanabe, T. Yano, Y. Shiokawa,  _Journal of The Electrochemical Society_ 2005, _152_ , 4 A830.
  * [64] X. Lu, O. Nilsson, J. Fricke, R. W. Pekala,  _Journal of Applied Physics_ 1993, _73_ , 2 581.
  * [65] C. N. Davies,  _Air Filtration_ ,  Academic, London, 1973.
  * [66] C. Zhu, Z. Qi, V. A. Beck, M. Luneau, J. Lattimer, W. Chen, M. A. Worsley, J. Ye, E. B. Duoss, C. M. Spadaccini, C. M. Friend, J. Biener,  _Science Advances_ 2018, _4_ , 8.
  * [67] J. Wicks, M. L. Jue, V. A. Beck, J. S. Oakdale, N. A. Dudukovic, A. L. Clemens, S. Liang, M. E. Ellis, G. Lee, S. E. Baker, E. B. Duoss, E. H. Sargent,  _Adv. Mater._ 2021, 2003855–8.
  * [68] D.-W. Chung, M. Ebner, D. R. Ely, V. Wood, R. Edwin García,  _Modelling Simul. Mater. Sci. Eng._ 2013, _21_ , 7 074009\.
  * [69] M. A. Tahir, H. Vahedi Tafreshi,  _Physics of Fluids_ 2009, _21_ , 8 083604.
  * [70] E. J. Wilson, C. J. Geankoplis,  _Ind. Eng. Chem. Fund._ 1966, _5_ , 1 9.
  * [71] C. Othmer,  _Int. J. Numer. Meth. Fluids_ 2008, _58_ , 8 861.
  * [72] K. Svandberg,  _International Journal for Numerical Methods in Engineering_ 1987, _24_ 359\.
  * [73] B. S. Lazarov, O. Sigmund,  _International Journal for Numerical Methods in Engineering_ 2010, _86_ , 6 765.
  * [74] J. S. Newman, C. W. Tobias,  _Journal of The Electrochemical Society_ 1962, _109_ , 12 1183.
  * [75] X. Ke, J. M. Prahl, J. I. D. Alexander, J. S. Wainright, T. A. Zawodzinski, R. F. Savinell,  _Chem. Soc. Rev._ 2018, _47_ , 23 8721.