# Stepwise ionization of Mo14+ ions in EBIT$:$ The importance of the
metastable level

Cunqiang Wu  Key Laboratory of Atomic and Molecular Physics $\&$ Functional
Materials of Gansu Province, College of Physics and Electronic Engineering,
Northwest Normal University, Lanzhou 730070, P.R. China  Xiaobin Ding
<EMAIL_ADDRESS>Key Laboratory of Atomic and Molecular Physics $\&$
Functional Materials of Gansu Province, College of Physics and Electronic
Engineering, Northwest Normal University, Lanzhou 730070, P.R. China  Gansu
International Scientific and Technological Cooperation Base of Laser Plasma
Spectroscopy, Lanzhou 730070,China  Qi Guo  Shanghai EBIT laboratory, and Key
Laboratory of Nuclear Physics and Ion-Beam Application (MOE), Institute of
Modern Physics, Fudan University, Shanghai 200433, China  Ke Yao  Shanghai
EBIT laboratory, and Key Laboratory of Nuclear Physics and Ion-Beam
Application (MOE), Institute of Modern Physics, Fudan University, Shanghai
200433, China  Jialin Liu  Shanghai EBIT laboratory, and Key Laboratory of
Nuclear Physics and Ion-Beam Application (MOE), Institute of Modern Physics,
Fudan University, Shanghai 200433, China  Yunqing Fu  Shanghai EBIT
laboratory, and Key Laboratory of Nuclear Physics and Ion-Beam Application
(MOE), Institute of Modern Physics, Fudan University, Shanghai 200433, China
Chenzhong Dong  Key Laboratory of Atomic and Molecular Physics $\&$ Functional
Materials of Gansu Province, College of Physics and Electronic Engineering,
Northwest Normal University, Lanzhou 730070, P.R. China  Gansu International
Scientific and Technological Cooperation Base of Laser Plasma Spectroscopy,
Lanzhou 730070,China

###### Abstract

The visible spectrum of Mo15+ ions was measured using a high-temperature
superconducting electron-beam ion trap at the Shanghai EBIT Laboratory, with
an electron beam energy $E_{e}$=400eV, significantly lower than the ionization
potential (IP=544.0 eV) of Mo14+ ions in ground state. To expound on the
experiment, the energy level structure, radiative transition properties,
electron-impact excitation, and electron-impact ionization cross section for
both the ground state and low-lying excited state of the Mo14+ ions were
calculated using Dirac–Fock-Slater method with a local central potential and
distorted wave approximation. The results demonstrated reasonable agreement
with both available experimental and theoretical data. Through an analysis of
the related atomic processes of Mo14+ ion, a scenario involving the stepwise
ionization of the metastable state 3p63d94s was proposed to explain the
presence of the Mo15+ ions with a lower energy of incident electron. Finally,
the significance of the metastable levels in the ionization of Mo14+ ions is
highlighted.

##  I Introduction

The study of the metastable level structures and properties of highly charged
ions (HCIs) hold importance, not only in testing atomic physics theory that
describes the interaction of atoms or ions with multipole radiative fields [1,
2, 3, 4, 5], but also the radiative emissions from these metastable levels
prove crucial for the effective application of density diagnostics in
astrophysical and laboratory plasmas [6, 7, 8]. The metastable levels of atoms
or ions are characterized by an inability to decay to lower energy levels via
electric dipole or Auger processes, constrained by the limitations of
selection rules, thereby extending the lifetime of the metastable level to
millisecond or more [9, 10, 11]. Certain long-lived metastable levels, which
create an unperturbed, well-defined, ideal two-level system, hold potential
applications in quantum information, plasma-assisted laser spectroscopy [12],
penning-trap mass spectrometry [13], and investigations into the variation of
fundamental constants [14, 15, 16]. Furthermore, long-lived metastable states
with high order forbidden transitions serve as valuable candidates for
studying hyperfine-induced transitions and magnetically-induced transitions,
as new decay channels can significantly influence their lifetimes [17, 18].

The long-lived metastable state can be considered analogous to the ground
state, thereby possessing similar properties and behaviors to the ground
state. In hot plasma environments, Novel electron-impact channels for
excitation, deexcitation, ionization, and recombination may open for the
metastable level [19]. In magnetic confinement fusion devices like the
International Thermonuclear Experimental Reactor (ITER), Tungsten (W) and
Molybdenum (Mo) serve as the divertor and plasma facing component materials.
Research has shown that the population of metastable states has a significant
affects measurements of tungsten erosion from plasma facing components and
impurity transport [20, 21]. Studies from the electron beam ion trap (EBIT)
further underscore that the presence of a long-lived metastable state with
reduced ionization energy can instigate a stepwise ionization process. This
process has a significant impact on establishment and progression of plasma’s
charge state equilibrium, as well as the characteristics of its emission
spectrum [22, 23].

The detection and manipulation of the charge state of HCIs are imperative for
a range of plasma applications. An EBIT [24, 25], chiefly comprising quasi-
monochromatic electron beams, serves as a powerful tool for investigating the
level structure, spectroscopic properties, and atomic processes in plasmas
with metastable state of HCIs [26, 27]. Recent experimental and theoretical
[23, 28, 29, 30, 31] research has indicated the undeniable significance of
diverse atomic processes induced by the metastable states of HCIs.
Investigations into several HCIs reveal that ionization processes initiated by
metastable states can be observed by monitoring the emission spectra, which
showcases incident electron energies significantly below the threshold
required for producing HCIs from the ground state of the preceding low charge
state ions [26, 32, 33, 27, 23, 28, 29, 30, 31]. Numerous ions possessing
metastable states have the potential to prompt complex ionization balance by
stepwise ionization pathways. While ionization of metastable states has been
investigated, such instances remain infrequent, thereby indicating a necessity
for a more systematic and quantitative analysis of stepwise ionization in
HCIs.

In this work, the Mo15+/Mo14+ ratio, the charge state distribution of
molybdenum ions, and the visible spectrum of Mo15+ ion are measured at the
Shanghai low-energy EBIT, CUbic-eBIT (CUBIT) [34, 35]. The observed visible
spectrum is identified as the M1 transition occurring between the doublet
levels of the ground configuration of Mo15+ ion. The ionization energy of
Mo14+ ions stands at 544.0 eV [36], and the experimental results demonstrate
that the M1 transition peak emerges where lower electron beam energy is
deployed. Consequently, by utilize the Dirac-Fock-Slaterrmcdhf (DFS) method
with a local central potential in the Flexible Atomic Code (FAC) [37, 38], we
compute the levels structure, radiative transition probability, and cross
sections for electron collisions of Mo15+ and Mo14+ ions. Additionally, the
GRASP2K code [39, 40], underpinned by the multi-configuration Dirac-Hartree-
Fock (MCDHF) theory, is also utilized to compute the transition wavelength.
Ultimately, a viable collisional-radiative model (CRM) was assembled and
resolved, taking the the metastable levels into account. Furthermore, we
postulated that the production of the Mo15+ ion could be attribute to the
ionization of metastable Mo14+ ions.

##  II Experimental setup

A schematic of the experimental setup is depicted in Fig.1. The setup
primarily consists of the CUBIT [35, 41, 42], a spectrograph, a charge-coupled
device (CCD) camera, and a charge state distribution (CSD) analyzer. CUIBT is
a room temperature permanent magnet electron beam ion trap where highly
charged ions are produced through electron impact ionization. Inside CUIBT, an
electron beam emitted from a LaB6 cathode is compressed by a 0.56T magnetic
field provided by a permanent magnet as it passes through the center drift
tube. The combination of the magnetic field and the electron beam space charge
provides radial trapping, while the voltage applied on the drift tubes forms
the axial electrostatic trapping potential. The electron beam energy is
controlled by adjusting the potential difference between the cathode and the
central drift tube, enabling the production of desired charged states of ions.

As a light source, visible spectra emitted from CUBIT can be observed with a
CCD camera mounted on a Czerny-Turner type visible spectrometer. Additionally,
ions in the trap region can be extracted in a pulse mode by lowering potential
on the end drift tube. These extracted ions are then transported to the CSD
analyzer, whrer ions with different charge to mass ratios are separated using
orthogonal electric and magnetic fields. Only ions with specific charge to
mass ratio can be detected by a microchannel plate (MCP). The CSD of ions is
measured by scanning the voltage applied on the analyzer.

The spectra of molybdenum ions measured in experiment are shown in Fig.2(a).
Simultaneously, the corresponding charge state distributions (CSDs) of
molybdenum ions were measured in a pulsed mode. In this mode, the operation of
the power system was regulated by the clock signal. When the end drift tube
was at high potential for about 1.6s, HCIs were created and trapped in the
middle drift tube, and the light emitted by the trapped ions was observed
during this period. Subsequently, the voltage on the end drift tube was
rapidly lowered to a low potential to open the trap. Therefore, the molybdenum
ions were promptly dumped and extracted, and a pulsed HCIs bunch was detected.
Fig.2(b) shows the measured CSDs at the same electron energy.

Figure 1:  Scheme of the experimental setup.

##  III Theory and computation

All pertinent atomic data, including energy levels, cross sections for
electron-impact excitation/ionization and radiative recombination, alongside
radiative transition probabilities, have been calculated utilizing FAC [37].
The frequency-dependent Breit interaction in the zero-energy limit for the
exchanged photon and higher-order QED effects such as self-energy and vacuum
polarization corrections (treated in the screened hydrogenic approximation)
were added in a subsequent relativistic configuration interaction calculation
[43]. The collisional-radiative model implemented in FAC has been used to
analyze the spectra and to predict the charge-state distributions [37].

The current CRM encompasses 450 and 212 excited levels, originating from
3p63d10, 3p63d${}^{9}4l$ ($l$=0–3), 3p63d${}^{8}nl^{2}$, 3p43d${}^{10}nl^{2}$
($n$=4–6, $l$=0–3), 3p53d${}^{9}4l5l$ ($l$=0–3), 3p53d${}^{9}4l^{2}$
($l$=0–3), 3s13p53d${}^{10}4l^{2}$ ($l$=0–3), and 3p63d9, 3p53d10,
3p63d${}^{8}nl$, 3p63d${}^{7}nl^{2}$ ($n$=4–6, $l$=0–3), 3p53d${}^{9}4l$
($l$=0–3) 3p53d${}^{9}nl$ ($n$=4,5, $l$=0–1), 3p53d${}^{8}4l^{2}$ ($l$=0–3),
3p63d${}^{7}4l5l^{\prime}$ ($l$=0–1, $l^{\prime}$=1–4 ) configurations for
Mo14+ and Mo15+, correspondingly. The populations of the excited levels are
derived by solving the rate equations for all excited states. In the current
model, we postulate the various levels of the ions are interconnected through
only collisional and radiative process within the plasma. Within the context
of low-density EBIT plasma, our primary considerations have been electron-
impact excitation and deexcitation, electron-impact ionization, radiative
transition, and radiative recombination. On the other hand, the other
processes such as dielectronic recombination, three-body recombination, charge
exchange, and ion escape can be disregarded. The following rate equation (1)
can be employed to describe the differential rate of the population of each
excited level:

$\displaystyle\frac{dN_{i}}{dt}=\sum_{\begin{subarray}{c}i\neq
j\end{subarray}}Q_{ij}N_{i}N_{e}+\sum_{i>j}A_{ij}N_{i}+N_{+}N_{e}R_{+j}-\sum_{\begin{subarray}{c}i\neq
j\end{subarray}}Q_{ji}N_{j}N_{e}\quad-\sum_{i<j}A_{ji}N_{j}-N_{j}N_{e}S_{j+}.$
(1)

where, Ne and Ni represent the electron density and the populations of the
i-th level. $A$ is the radiative transition probabilities (including E1, E2,
E3, M1, M2, and M3). $Q$, $S$, and $R$ denote the electron-impact excitation
or deexcitation, electron-impact ionization, and radiative recombination rate
coefficients, respectively. These rate coefficients can be derived by by
convoluting the corresponding cross sections on the free electron energy
distribution function, described by the $\delta$ function for the
monoenergetic electron beam of EBIT. The cross section are computed using the
distorted wave approximation, while the inverse process are inferred in
accordance with the detailed balance principle. Given the quasi-steady-state
approximation $\frac{dN_{i}}{dt}=0$ and normalization condition
$\sum_{i}N_{i}=1$, the population of each level can be ascertained within the
specific electron density and temperature.

##  IV Results and discussion

Figure 2: (a) Spectra of molybdenum ions around 370 nm recorded at electron-
beam energies of 400, 460, 500, and 580 eV. The line at 370.81 nm is the
3p63d9 (2D3/2$\rightarrow$2D5/2) M1 transition of the Mo15+ ion. (b) CSDs of
Mo ions for the same electron energies. Black lines are the results of
measurements, the solid red and blue lines are the Gaussian fits used to infer
the CSDs, and the dashed red lines show the baselines of the fits. Dashed
purple lines indicates the centroid of the specific charged state.

Fig.2(a) depicts the experimental spectra of molybdenum ions observed at
electron-beam energies of 400, 460, 500and 580 eV, respectively. It is evident
from Figure 2(a) that an isolated transition peak of the Mo15+ ion is observed
at a beam energy of 400eV with the wavelength of 370.81nm. Furthermore,
Fig.2(b) illustrates that the abundance of Mo15+ predominates during this
period. Utilizing the MCDHF and DFS methods, we estimate the transition
wavelength and compare it with extant experimental [36] and theoretical [44,
45] results, as presented in Table 1. As illustrated in the table, the
wavelength of the observed transition peak principally corresponds to the M1
transition wavelength of the ground state 3p63d9 (2D5/2$\rightarrow$2D3/2) of
the Mo15+ ion. Notwithstanding, while the ionization energy of Mo14+ ions is
approximately 544.0 eV [36], experimental observation have depicted that Mo15+
ions exist at significantly lower ionization energy, for instance when
$E_{e}$= 400eV. In order to elucidate this phenomenon, a hypothetical
framework of stepwise ionization of the metastable state 3p6$3d^{9}4s$ has
been proposed to account for the presence of Mo15+ ions with lower incident
electron energy.

Table 1:  Wavelength of the 3p63d9 2D5/2$\rightarrow$2D3/2 transition in Co-like Mo. This work (nm) |  | Other work (nm)
---|---|---
MCDHF | DFS | NIST[36] | RMBPT[44] | RCC[45]
371.85 | 372.87 | 370.81 | 371.29 | 371.30
Figure 3: Energy level structure of Mo14+ ion 3p63d10, 3p63d94s, 3p63d94p and Mo15+ ion 3p63d9. The dash line represents electron-impact excitation, and the solid line represents radiation decay, the shot dot dash line represents collision ionization. Table 2:  The excitation energy (in:eV) of the metastable state 3p63d94s of Mo14+ ion and the ionization energy of ground state 3p63d9 of Mo15+ ion. Error=$\frac{(DFS-NIST)}{NIST}\times 100\%$ |  |  | Present |  | Reference |  |
---|---|---|---|---|---|---|---
Charge state | Config | Level | DFS | NIST[36] | MCDF[46] | RPTMP[47] | Error
Mo14+ | 3p63d10 |  1S0 | 0 | 0 | 0 | 0 |
| 3p63d94s |  3D3 | 210.16 | 210.14 | 210.14 | 210.15 | 0.01%
|  |  3D2 | 210.78 | 210.76 | 210.75 | 210.76 | 0.01%
|  |  3D1 | 213.46 | 213.47 | 213.47 | 213.47 | 0.01%
|  |  1D2 | 214.04 | 214.05 | 214.04 | 214.05 | 0.01%
Mo15+ | 3p63d9 |  2D5/2 | 543.88 | 544.00 |  |  | 0.02%
|  |  2D3/2 | 547.20 | 547.34 |  |  | 0.03%

Fig.3 illustrates the possible pathway for the production of Mo15+ ions from
Mo14+ ions. In Fig.3, the dash-dot line symbolizes the electron-impact
excitation process, while the short dash-dot line corresponds to the electron-
impact ionization process, and the solid line stands for the 3p63d${}^{9}4p$
radiative decay process. As observed in Fig.3, the the ionization energy
required to produce Mo15+ from the metastable state 3p63d${}^{9}4s$ is
approximately 338 eV, which is less than the incident electron beam energy
corresponding to the Mo15+ ion spectral line used in the experiment, rendering
this process of energetically feasible. Concurrently, Table 2 displays the
excitation energy for the low-lying excited configuration 3p63d94s of the
Mo14+ ion and the ionization energy of ground state configuration 3p63d9 of
Mo15+ ion. The current results, obtained using the DFS methods, align
remarkably well with NIST data [36] and other theoretical results [46, 47].
The maximum relative error of 0.03% indicates the inclusion of significant
configuration interactions in these calculations.

Fig 4 shows the electron-impact excitation cross section of the Mo14+ ion
3p63d10–3p63d94s (a) and 3p63d10–3p63d94p (b) as a function of incident
electron energy. Both the electron-impact excitation cross section from the
ground state 3p63d10 to the metastable state 3p63d94s and low-lying excited
state 3p63d94p decrease monotonically with the increase in electron energy.
However, Fig.4(b) reveals that the cross section from the ground state 3p63d10
(1S${}_{0})$ to the low-lying excited state 3p63d94p (1P1) is approximately
three times larger than the cross section to 3p63d94s, while the rest of the
cross sections are comparatively small. This observation suggests that the
likelihood of 1S0 being excited to 1P1 by the electron-impact exceeds that of
other states.

Figure 4: The electron-impact excitation cross section from the ground state
3p63d10 to the low-lying excited state 3p63d94s (a) and the metastable state
3p63d94p (b) of Mo14+ ions varies with electron energy.

However, the low-lying excited state 3p63d94p will decay via E1 transition
both to the ground state 3p63d10 (1S0) and the metastable state 3p63d94s
(3D3,2,1,1D2). The decay channels of the Mo14+ ion’s low-lying excited state
and metastable state are portrayed in Fig.5. In addition to the E1 transition
to the ground state 3p63d10, the low-lying excited state 3p63d94p also
exhibits E1 decay to the metastable state 3p63d94s. The metastable states,
including ${}^{3}D_{2}$ and ${}^{1}D_{2}$ exhibits E2 decay to the ground
state, while ${}^{3}D_{1}$ demonstrate M1 decay to both ground state and
internal energy levels 3$D_{2}$. However, 3D3 shows only a weak M3 transition,
which results in a substantial population of ${}^{3}D_{3}$.

Figure 5:  Decay channels from low-lying excited state 3p63d94p of the Mo14+
ions to the ground state 3p63d10 and the metastable state 3p63d94s.

Fig.6 represents the electron-impact ionization cross section of the
metastable state $3p^{6}3d^{9}4s$ of Mo14+ ion. As the electron beam energy
increases, the electron-impact ionization cross section gradually increases.
Moreover, at the same electron energy, the 3D3/3D2–2D5/2 and 3D1/1D2–2D3/2
electron impact ionization cross section are an order of magnitude larger than
other cross sections, and 3D3–2D5/2 electron-impact ionization cross section
is the largest. This demonstrates the increased likelihood of Mo15+ ions
progressing through this channel.

Figure 6: The electron-impact ionization cross section of the Mo14+ ion metastable state 3p63d94s varies with incident electron energy. Table 3: The lifetime and population of the Mo14+ ion metastable 3p63d94s. The computation is conducted at electron-beam energy 400 eV and density $1.0\times 10^{12}$/cm3. Upper | Lower |  $\lambda$ (nm) | Referce | Type | A${}_{ij}(s^{-1})$ | Population$(\%)$ | Lifetime(ms)
---|---|---|---|---|---|---|---
${}^{1}S_{0}$ | $-$ | $-$ | $-$ | $-$ | $-$ | 96.24 | $-$
${}^{3}D_{2}$ | ${}^{1}S_{0}$ | 5.85 | 5.88[36] | E2 | 9.60E+06 | 0.005 | 1.04E$-$04
${}^{3}D_{1}$ | ${}^{1}S_{0}$ | 5.78 | 5.83[48] | M1 | 4.88E$-$01 | 0.45 | 4.13E+00
| ${}^{3}D_{2}$ |  |  | M1 | 2.42E+02 |  |
${}^{1}D_{2}$ | ${}^{3}D_{3}$ |  |  | M1 | 3.22E+02 | 0.001 | 7.30E$-$05
| ${}^{3}D_{2}$ |  |  | M1 | 2.00E+01 |  |
| ${}^{1}S_{0}$ | 5.77 | 5.79[36] | E2 | 1.37E+07 |  |
${}^{3}D_{3}$ | ${}^{1}S_{0}$ | 5.87 | 5.92[48] | M3 | 3.83E$-$01 | 3.30 | 2.85E+03

Table 3 presents the transition wavelength, transition probability, lifetime,
and population of the Mo14+ ion’s metastable state 3p63d94s. The populations
of each metastable levels are estimated by the collisional-radiative model,
using an electron-beam energy 400 eV and density $1.0\times 10^{12}$/cm3. As
observed in Table 3, 3D3 possesses both a high population and a longer
lifetime, and also exhibits the greatest electron-impact ionization cross
section. The large cross section and high population increase the likelihood
that ions will reach to the Mo15+ ion’s ground state through this route.

Based on the foregoing analysis, the presence of Mo15+ ions, which were
experimentally observed at a lower incident electron energy, can be attributed
to the process where the ground state of Mo14+ ion is collisionally excited to
the metastable state 3p63d94s and low-lying excited state 3p63d94p by an
electron. Subsequently, the low-lying excited state undergoes E1 decay to the
metastable state. Furthermore, the metastable state is then collisionally
ionized by an electron, resulting in the production of Mo15+ ions.

Figure 7:  The electron-impact ionization rate and total ionization rate of
the Mo14+ ion in both the ground 3p63d10 (1S0) state and the metastable state
3p63d94s (3D3).

Fig.7 illustrates the electron-impact ionization rate for both the metastable
state and ground state, as well as the total ionization rate of Mo14+ ions.
Below the ionization threshold of Mo14+ ions (IP=554.0 eV), it is solely the
stepwise ionization process of the metastable state that contributes to the
total ionization rate. When over the ionization threshold, both the stepwise
ionization of the metastable state and the direct ionization of the ground
state simultaneously contribute to the total rate. The total ionization rate
is determined by weighted summation of the electron-impact ionization rate for
the metastable $(3.30\%)$ state and the ground state $(96.24\%)$.

The experiment also observed the charge state ratio of Mo15+ and Mo14+ ions.
When the incident electron energy is significantly lower than the ionization
energy of Mo14+ ion, the ratio between the two ions exceeds 1, a phenomenon
primarily attributed to metastable stepwise ionization. Upon considering the
various atomic processes associated with the aforementioned excitation and
ionization channels, a suitable CRM was assembled to analyze the dependency of
the ratio of these two ions on the incident electron energy. Good agreement
was obtained between the experimental and the theoretical results, as
portrayed in Fig.8. Moreover, it is observed that when the incident electron
energy surpasses the ionization energy, the ratio between the two ions
experiences a sudden increase. This can be primarily attributed to the
simultaneous ionization of the metastable and ground state, leading to an
increased total ionization rate.

Figure 8:  Comparison of Mo15+/Mo14+ ion ratio between theoretical calculation
and EBIT experimental measurement.

##  V Summary

In this study, both the Mo15+/Mo14+ ratio, the charge state distribution of
Mo14+ and Mo15+ ions and visible spectrum of Mo15+ ion were measured at the
Shanghai low-energy EBIT, CUbic-eBIT (CUBIT). The energy level, radiative
transitions probability, electron-impact excitation and electron-impact
ionization cross sections of Mo14+ and Mo15+ ions were calculated using the
Dirac-Fock-Slater method with a local central potential and distorted-wave
approximation. All these calculation were based on the FAC package. Through a
detailed analysis of various atomic processes, a theoretically most likely to
occur channel was identified. A reasonable explanation was provided for the
experimental observation of the Mo15+ ion spectrum below the ionization energy
of Mo14+ ion (544.0 eV). The Mo15+/Mo14+ ions ratio was calculated using the
collisional-radiative model, which considered the metastable state related
process, with the resulting calculations showing strong agreement with the
experiment.

##  VI ACKNOWLEDGMENTS

This work was supported by National Nature Science Foundation of China (Grant
No. 12274352), the National Key Research and Development Program of China
(Grant No. 2022YFA1602500).

## References

  * Ralchenko  _et al._ [2011] Y. Ralchenko, I. N. Draganić, D. Osin, J. D. Gillaspy, and J. Reader, Phys. Rev. A 83, 032517 (2011).
  * Guise  _et al._ [2014] N. D. Guise, J. N. Tan, S. M. Brewer, C. F. Fischer, and P. Jönsson, Phys. Rev. A 89, 040502 (2014).
  * Saito  _et al._ [2015] M. Saito, N. Ojima, S. Itoi, and Y. Haruyama, Phys. Rev. A 91, 012508 (2015).
  * Müller  _et al._ [2018] A. Müller, E. Lindroth, S. Bari, A. Borovik, P.-M. Hillenbrand, K. Holste, P. Indelicato, A. L. D. Kilcoyne, S. Klumpp, M. Martins, J. Viefhaus, P. Wilhelm, and S. Schippers, Phys. Rev. A 98, 033416 (2018).
  * Sakaue  _et al._ [2019] H. A. Sakaue, D. Kato, I. Murakami, H. Ohashi, and N. Nakamura, Phys. Rev. A 100, 052515 (2019).
  * Moehs  _et al._ [2000] D. P. Moehs, D. A. Church, M. I. Bhatti, and W. F. Perger, Phys. Rev. Lett. 85, 38 (2000).
  * Steinbrügge  _et al._ [2022] R. Steinbrügge, S. Kühn, F. Nicastro, M. F. Gu, M. Togawa, M. Hoesch, J. Seltmann, I. Sergeev, F. Trinter, S. Bernitt, C. Shah, M. A. Leutenegger, and J. R. C. López-Urrutia, Astrophys. J. 941, 188 (2022).
  * Chen  _et al._ [2023] Y. Chen, W. Li, H. Tian, X. Bai, R. Hutton, and T. Brage, Res. Astron. Astrophys. 23, 022001 (2023).
  * Kim and Desclaux [2002] Y.-K. Kim and J.-P. Desclaux, Phys. Rev. A 66, 012708 (2002).
  * Brenner  _et al._ [2007] G. Brenner, J. R. C. López-Urrutia, Z. Harman, P. H. Mokler, and J. Ullrich, Phys. Rev. A 75, 032504 (2007).
  * Brandt  _et al._ [2022] A. D. Brandt, S. F. Cooper, C. Rasor, Z. Burkley, A. Matveev, and D. C. Yost, Phys. Rev. Lett. 128, 023001 (2022).
  * Schüssler  _et al._ [2020] R. X. Schüssler, H. Bekker, M. Braß, H. Cakir, J. R. C. López-Urrutia, M. Door, P. Filianin, Z. Harman, M. W. Haverkort, W. J. Huang, P. Indelicato, C. H. Keitel, C. M. König, K. Kromer, M. Müller, Y. N. Novikov, A. Rischka, C. Schweiger, S. Sturm, S. Ulmer, S. Eliseev, and K. Blaum, Nature 581, 42 (2020).
  * Kimura  _et al._ [2023a] N. Kimura, Priti, Y. Kono, P. Pipatpakorn, K. Soutome, N. Numadate, S. Kuma, T. Azuma, and N. Nakamura, Communications Physics 6 (2023a).
  * Safronova  _et al._ [2014a] M. S. Safronova, V. A. Dzuba, V. V. Flambaum, U. I. Safronova, S. G. Porsev, and M. G. Kozlov, Phys. Rev. Lett. 113, 030801 (2014a).
  * Safronova  _et al._ [2014b] M. S. Safronova, V. A. Dzuba, V. V. Flambaum, U. I. Safronova, S. G. Porsev, and M. G. Kozlov, Phys. Rev. A 90, 052509 (2014b).
  * Safronova  _et al._ [2014c] M. S. Safronova, V. A. Dzuba, V. V. Flambaum, U. I. Safronova, S. G. Porsev, and M. G. Kozlov, Phys. Rev. A 90, 042513 (2014c).
  * Schippers  _et al._ [2007] S. Schippers, E. W. Schmidt, D. Bernhardt, D. Yu, A. Müller, M. Lestinsky, D. A. Orlov, M. Grieser, R. Repnow, and A. Wolf, Phys. Rev. Lett. 98, 033001 (2007).
  * Träbert  _et al._ [2007] E. Träbert, P. Beiersdorfer, and G. V. Brown, Phys. Rev. Lett. 98, 263001 (2007).
  * Kobayashi  _et al._ [2014] Y. Kobayashi, D. Kato, H. A. Sakaue, I. Murakami, and N. Nakamura, Phys. Rev. A 89, 010501 (2014).
  * Quinet  _et al._ [2010] P. Quinet, V. Vinogradoff, P. Palmeri, and E. Biémont, J. Phys. B: At. Mol. Opt. Phys. 43, 144003 (2010).
  * Johnson  _et al._ [2020] C. A. Johnson, S. D. Loch, and D. A. Ennis, Plasma Phys. Control. Fusion 62, 125017 (2020).
  * Li  _et al._ [2015] W. Li, Z. Shi, Y. Yang, J. Xiao, T. Brage, R. Hutton, and Y. Zou, Phys. Rev. A 91, 062501 (2015).
  * Lu  _et al._ [2019] Q. Lu, J. He, H. Tian, M. Li, Y. Yang, K. Yao, C. Chen, J. Xiao, J. G. Li, B. Tu, and Y. Zou, Phys. Rev. A 99, 042510 (2019).
  * Fu  _et al._ [2010] Y. Fu, K. Yao, B. Wei, D. Lu, R. Hutton, and Y. Zou, J. Instrum. 5 (08), C08011.
  * Nakamura  _et al._ [2008] N. Nakamura, H. Kikuchi, H. A. Sakaue, and T. Watanabe, Rev. Sci. Instrum. 79, 063104 (2008).
  * Sakoda  _et al._ [2011] J. Sakoda, A. Komatsu, H. Kikuchi, and N. Nakamura, Phys. Scr. T144, 014011 (2011).
  * Torretti  _et al._ [2017] F. Torretti, A. Windberger, A. Ryabtsev, S. Dobrodey, H. Bekker, W. Ubachs, R. Hoekstra, E. V. Kahl, J. C. Berengut, J. R. C. López-Urrutia, and O. O. Versolato, Phys. Rev. A 95, 042503 (2017).
  * Kimura  _et al._ [2020] N. Kimura, R. Kodama, Priti, N. Numadate, K. Suzuki, M. Monobe, and N. Nakamura, Phys. Rev. A 102, 032807 (2020).
  * Lu  _et al._ [2021] Q. Lu, C. L. Yan, J. Meng, G. Q. Xu, Y. Yang, C. Y. Chen, J. Xiao, J. G. Li, J. G. Wang, and Y. Zou, Phys. Rev. A 103, 022808 (2021).
  * Yan  _et al._ [2022] C. L. Yan, Q. Lu, Y. M. Xie, B. L. Li, N. Fu, Y. Zou, C. Chen, and J. Xiao, Phys. Rev. A 105, 032820 (2022).
  * Kimura  _et al._ [2023b] N. Kimura, Priti, N. Numadate, R. Kodama, and N. Nakamura, Phys. Rev. A 108, 032818 (2023b).
  * Qiu  _et al._ [2014] M. L. Qiu, R. F. Zhao, X. L. Guo, Z. Z. Zhao, W. X. Li, S. Y. Du, J. Xiao, K. Yao, C. Y. Chen, R. Hutton, and Y. Zou, J Phys B: At Mol Opt Phys 47, 175002 (2014).
  * Windberger  _et al._ [2016] A. Windberger, F. Torretti, A. Borschevsky, A. Ryabtsev, S. Dobrodey, H. Bekker, E. Eliav, U. Kaldor, W. Ubachs, R. Hoekstra, J. R. Crespo López-Urrutia, and O. O. Versolato, Phys. Rev. A 94, 012506 (2016).
  * Xiao  _et al._ [2012] J. Xiao, Z. Fei, Y. Yang, X. Jin, D. Lu, Y. Shen, L. Liljeby, R. Hutton, and Y. Zou, Rev. Sci. Instrum. 83, 013303 (2012).
  * He  _et al._ [2022] Z. He, J. Meng, Y. Li, F. Jia, N. Khan, B. Niu, L. Huang, Z. Hu, J. Li, J. Wang, _et al._ , J. Quant. Spectrosc. Radiat. Transf , 108276 (2022).
  * Kramida  _et al._ [2021] A. Kramida, Yu. Ralchenko, J. Reader, and and NIST ASD Team, NIST Atomic Spectra Database (ver. 5.9), [Online]. Available:[2021, December 21]. National Institute of Standards and Technology, Gaithersburg, MD. (2021).
  * Gu [2008] M. F. Gu, The flexible atomic code, Can. J. Phys. 86, 675 (2008).
  * Gu [2004] M. F. Gu, Phys. Rev. A 70, 062704 (2004).
  * Jönsson  _et al._ [2013] P. Jönsson, G. Gaigalas, J. Bieroń, C. F. Fischer, and I. Grant, Comput. Phys. Commun. 184, 2197 (2013).
  * Froese Fischer  _et al._ [2019] C. Froese Fischer, G. Gaigalas, P. Jönsson, and J. Bieroń, Comput. Phys. Commun. 237, 184 (2019).
  * Li  _et al._ [2022] Y. Li, Y. Wang, J. Fan, R. Si, J. Li, M. Zhang, L. Huang, J. Xiao, Y. Zou, B. Wei, and K. Yao, J. Phys. B: At. Mol. Opt. Phys. 54, 235001 (2022).
  * Wang  _et al._ [2022] Y. Wang, Y. Li, J. Liu, F. Jia, R. Si, M. Zhang, L. Huang, B. Tu, Y. Zou, B. Wei, and K. Yao, J. Quant. Spectrosc. Radiat. Transf 293, 108370 (2022).
  * Kozioł [2020] K. Kozioł, J. Quant. Spectrosc. Radiat. Transf 242, 106772 (2020).
  * Guo  _et al._ [2016] X. L. Guo, R. Si, S. Li, M. Huang, R. Hutton, Y. S. Wang, C. Y. Chen, Y. M. Zou, K. Wang, J. Yan, C. Y. Li, and T. Brage, Phys. Rev. A 93, 012513 (2016).
  * Nandy and Sahoo [2021] D. K. Nandy and B. K. Sahoo, Phys. Rev. A 104, 052812 (2021).
  * Brage and Litzén [2006] T. Brage and U. Litzén, Phys. Scr. 35, 662 (2006).
  * Ivanova [2021] E. Ivanova, At. Data. Nucl. Data Tables 139, 101413 (2021).
  * Safronova  _et al._ [2006] U. Safronova, A. Safronova, S. Hamasha, and P. Beiersdorfer, At. Data. Nucl. Data Tables 92, 47 (2006).


# Topology optimization of 3D flow fields for flow batteries

Tiras Y. Lin, Sarah E. Baker, Eric B. Duoss, & Victor A. Beck  Lawrence
Livermore National Laboratory, Livermore, CA 94550, United States
<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

As power generated from renewables becomes more readily available, the need
for power-efficient energy storage devices, such as redox flow batteries,
becomes critical for successful integration of renewables into the electrical
grid. An important component in a redox flow battery is the planar flow field,
which is usually composed of two-dimensional channels etched into a backing
plate. As reactant-laden electrolyte flows into the flow battery, the channels
in the flow field distribute the fluid throughout the reactive porous
electrode. We utilize topology optimization to design flow fields with full
three-dimensional geometry variation, i.e., 3D flow fields. Specifically, we
focus on vanadium redox flow batteries and use the optimization algorithm to
generate 3D flow fields evolved from standard interdigitated flow fields by
minimizing the electrical and flow pressure power losses. To understand how
these designs improve performance, we analyze the polarization of the reactant
concentration and exchange current within the electrode to highlight how the
designed flow fields mitigate the presence of electrode dead zones. While
interdigitated flow fields can be heuristically engineered to yield high
performance by tuning channel and land dimensions, such a process can be
tedious; this work provides a framework for automating that design process.

  * February 2022

##  1 Introduction

The need for effective, large-scale energy storage continues to grow as the
technology driving renewables continues to mature [1, 2, 3]. Renewable
sources, such as wind and solar, are inherently weather and time dependent.
Unfortunately, this intermittent energy supply is not always aligned with user
needs for electricity, necessitating the use of power plants to meet peak
demand [4, 5]. The mismatch between supply and demand hampers our ability to
completely integrate renewables into the electrical grid. A stable, integrated
electrical grid requires effective, large-scale energy storage – that is,
efficient conversion from electrical energy to potential (e.g. chemical)
energy, and vice versa – to buffer against this mismatch by storing and
dispatching any excess generated electricity as necessary [6, 3].
Additionally, to be industrially relevant, these storage devices must operate
with high power efficiency while remaining economically practical [7]. Redox
flow batteries inherently decouple energy storage from power production,
allowing the former to scale independently by merely increasing the size of
the external storage vessels. As a result, redox flow batteries are an
especially promising energy storage system to meet these needs [8, 9, 10, 11].

A critical problem in designing large-scale electrochemical systems, and
energy storage systems such as redox flow batteries in particular, is the
distribution of reactants to the electrode. The electrode is a porous media
with a large exposed surface area of electrochemically active material [12].
Often, the electrode is made of a stochastic, graphitic or vitreous carbon
material cohered into a bulk material such as a paper or felt [8]. Ideally,
the electrolyte flow to the electrode will minimize excess power losses while
appropriately distributing the reactant to promote reaction and mitigate the
emergence of dead zones (i.e., low concentration regions with insufficient
mass transfer) [13, 14, 15, 16, 17]. This problem is usually addressed by
using sophisticated fluid distribution systems both external and internal to
the electrochemical cell stack. The internal fluid distribution system is
known as a flow field [14, 18], and significant effort has been devoted to
understanding how the flow field architecture impacts device performance [19,
20, 21, 22, 23, 24, 25, 13].

Two of the most prevalently used flow field architectures are serpentine
channel designs, where a single winding flow channel connects the electrolyte
inlet to the outlet [26, 27, 28, 29, 30, 31, 32], and interdigitated channel
designs, composed of a series of interlocking channels spanning the electrode
[33, 34, 35, 36, 37]. A key difference between the two is that interdigitated
flow fields force all of the electrolyte through the electrode while
serpentine designs allow fluid bypass. Thus, for the latter, only a fraction
of the flowing reactant penetrates the electrode [14]. Gerhardt et al. [33]
investigated the interdigitated flow field both computationally and
experimentally. They showed how losses due to fluid pumping and the electrical
overpotential change with the widths of the fluid channels and solid lands.
While the former is generally reduced with larger fluid channel width, the
latter is generally reduced with larger solid land width. Furthermore, for
their specific system, they identified channel and land widths that lead to a
minimum in overall power loss, though these dimensions changed as the
operating conditions were varied. While this, and related, work have
elucidated many of the associated transport mechanisms, there remains much
debate about choosing an optimal flow field geometry [38, 39]. This difficult
question is exacerbated by the vast design space possible for each
architecture, including the recognition that the optimal flow field could
change depending on the chemistry employed, the operating conditions, and the
scale, etc.

Indeed, recent work is now beginning to focus on more widely exploring this
design space by studying flow fields which deviate from standard 2D (i.e., a
geometry that only varies in two dimensions) interdigitated and serpentine
configurations [40, 41, 42, 43]. This includes introducing hierarchy [41],
splitting and reconfiguration [42], and also the inclusion of 3D features such
as ramps and obstructions [43]. A central focus of these efforts is to
increase performance at larger scales where standard flow field designs may be
less effective [38, 41, 43].

Despite their novelty, these architectures remain predominantly two
dimensional, driven in large part by manufacturing constraints. However,
recent advances in additive and advanced manufacturing offer viable techniques
for fabricating complex 3D parts for electrochemical systems including flow
fields [44, 45, 46, 47, 48, 49]. Additively manufactured, porous electrodes
which include 3D features have been demonstrated to lead to improved mass
transfer at higher flow rates [50]. There is also evidence that introducing
three-dimensional features into the flow field itself can provide high
performance at industrially relevant scales, e.g., as demonstrated by the
design of the flow field used in the Toyota Mirai, a fuel cell vehicle [51].

However, even as additive and advanced manufacturing technologies improve, the
community lacks guidance on what 3D features to consider for these 3D flow
fields, i.e., flow fields with significant geometry variation in all three
dimensions. This difficulty can be attributed to the complex nature of the
underlying physical phenomena that couples together simultaneous fluid
transport, mass transfer, and electrostatics. While it would be possible to
iterate on 3D designs using a combination of intuition and simulation, such a
process would be tedious, computationally expensive, and slow. Instead, we
propose to use topology optimization [52] to computationally guide flow field
design. This method has seen significant advancements in recent years, and
researchers have used topology optimization to aid in the design of heat
exchangers [53, 54, 55, 56, 57, 58], structural supports [59, 60], and
microfluidic devices [61, 62]. Additionally, topology optimization has been
applied to the design of two-dimensional flow fields for both flow batteries
[63, 64] and fuel cells [65]. More recently, Beck et al. [66] demonstrated the
application of three-dimensional optimization to electrochemical transport, by
designing an optimized and homogenized variable porosity electrode. Here, we
leverage the ideas of that work, and develop a topology optimization method to
design 3D flow fields.

In this work, we aim to use topology optimization, in tandem with a high-
resolution electrochemical model of a redox flow battery, to design 3D flow
fields. We specifically focus on the well-studied vanadium redox flow battery
[67, 68, 69, 70, 71, 72], the most widely commercialized redox flow battery
[7]. We seek to improve upon the popular interdigitated design, and in
particular, we aim to showcase the pertinent three-dimensional features that
arise. We perform the optimization to reduce the sum of electrical and fluid
flow pumping power losses, and we identify how the qualitative design features
of the optimized design change as the operating conditions of the flow battery
change. We quantify the power efficiencies of the optimized designs, showcase
how different optimal designs appear depending on the initial design provided
to the optimization algorithm, and discuss how the optimized designs more
evenly distribute reactant to the electrode. Since flow fields, and flow
manifolds in general, are necessary for any engineered system requiring fluid
and mass distribution, we also provide suggestions for how our work can be
extended in the future.

##  2 Methods

### Problem overview

Figure 1: (A) Schematic of the negative half-cell of the flow battery setup. A
porous electrode is sandwiched between the membrane and flow field plate.
Shown below are traditional interdigitated flow fields with various land and
channel widths. In this work, we apply topology optimization and
electrochemical modeling to further improve these designs. (B-E) For an
example case, various iterations of the flow field designs are shown during
the optimization process, corresponding to the points indicated on the plot of
power loss as a function of optimization iteration shown in panel (F). To aid
in visualization, in panels (B-E), half of the solid flow field is shown, and
the negative of this, i.e. the fluid channels that form, is shown for the
other half.

We consider the negative half-cell of a vanadium redox flow battery as shown
in figure 1(A). Liquid electrolyte, comprised of a solution of $V^{2+}$ and
$V^{3+}$ at a concentration $C_{\mathrm{inlet}}=1\mathrm{M}$ in 1M sulfuric
acid flows in, which is then guided by the flow field towards the porous
carbon-felt electrode, where the reaction $V^{2+}\rightarrow V^{3+}+e^{-}$
occurs on the surfaces of electrode fibers. The spent fluid then flows from
the outlet. Inspired by the popularity and success of the interdigitated flow
field design (examples of which are shown in the bottom of figure 1(A)), we
utilize a combination of a continuum model of the flow battery half-cell
following porous electrode theory [66, 73] – incorporating fluid mechanics,
mass transfer, and electrostatics – and topology optimization to improve upon
this design. An example of the optimization process used in this work is shown
in figures 1(B-E). Specifically, the flow field designs at various iterations
in the optimization process are shown, corresponding to points indicated on
the plot of power loss plotted against iteration number in figure 1(F). To aid
in visualization, in figures 1(B-E), half of the solid flow field is shown,
and the negative of the solid part, i.e. the fluid channels within the flow
field, are shown for the other half. In this work, we focus on a half-cell
with experimentally relevant dimensions: the square porous carbon electrode
has side length of $L=2$ cm, and thus an area of $A=4$ $\mathrm{cm^{2}}$, and
is $500$ $\mathrm{\mu m}$ thick, and the flow field is 3 mm thick.

We aim to design the flow field to minimize the excess power losses $P$ that
occur for a specified discharge current $I$ and electrolyte flow rate $Q$. We
consider in this work current densities $I/A=\\{50,100,150\\}$
$\mathrm{mA/cm^{2}}$ and flow rates $Q=\\{2.4,12,24,36\\}$ $\mathrm{mL/min}$
[71, 67, 33, 69, 72]. The total power loss is defined as

$P=\int\limits_{\mathrm{membrane}}\eta\frac{I}{A}\,\rmd\mathbf{x}+\frac{1}{\Psi_{\mathrm{pump}}}\int\limits_{\mathrm{inlet}}p\bar{u}_{\mathrm{inlet}}\,\rmd\mathbf{x},$
(1)

where $\eta$ is the half-cell overpotential, $I/A$ is the discharge current
density, $\Psi_{\mathrm{pump}}$ is the pump efficiency, which is assumed to be
$0.5$ in this work, $p$ is the fluid pressure, and
$\bar{u}_{\mathrm{inlet}}=Q/A_{\mathrm{inlet}}$ is the specified averaged
inlet fluid velocity. The first term represents the electrical power loss due
to the half-cell overpotential, while the second term represents the fluid
pumping power loss. As highlighted by Gerhardt et al. [33], a flow field
design with more solid present will typically have a reduced Ohmic electrical
power loss due to the high conductivity of the solid, and a flow field design
with more void will typically have a reduced fluid pumping power loss due to
the reduced hydraulic resistance – this paradox highlights why optimization of
this component is nontrivial. Equivalently, our goal is to maximize the
overall power efficiency of the half-cell, which is defined as

$\Xi=1-\frac{P}{IU_{0}},$ (2)

where $U_{0}$ is the equilibrium half-cell potential. As defined, $\Xi=1$ for
a perfectly efficient cell and can be $<0$ for a very poorly designed cell.

### Continuum simulation

The equations governing the continuum electrochemical model are as follows.
The fluid velocity $\mathbf{u}$ is described by the steady incompressible
Navier-Stokes equation with a Darcy drag,

$\displaystyle\rho\mathbf{u}\cdot\nabla\mathbf{u}$
$\displaystyle=-\mu\alpha^{-1}\mathbf{u}-\nabla p+\mu\nabla^{2}\mathbf{u},$
(3) $\displaystyle\nabla\cdot\mathbf{u}$ $\displaystyle=0,$ (4)

where $\rho$ is the fluid density, $\mu$ is the viscosity, $p$ is the
pressure, and $\alpha$ is the local permeability of the material. As described
below, the topology optimization technique models the presence of solid and
liquid material by smoothly varying the parameters using a position dependent
design variable, $\gamma(\mathbf{x})$. In the solid, $\alpha^{-1}$ takes on a
large value to prevent fluid flow, in the liquid regions, $\alpha^{-1}=0$, and
in the fixed porous region representing the electrode, $\alpha=\alpha_{E}$.
This parameter thus controls the location of solid and liquid within the flow
field and approximates the no-slip condition for the fluid velocity. At the
inlet, we specify the velocity to be
$\mathbf{u}=-(Q/A_{\mathrm{inlet}})\mathbf{n}$, and the fluid exits to zero
pressure; no-slip is assumed at all other boundaries.

The mass transfer of the species concentration $C_{i}$, where
$i\in\\{V^{2+},V^{3+}\\}$, is described by the steady advection-diffusion-
reaction equation, given by,

$\mathbf{u}\cdot\nabla C_{i}=\nabla\cdot(D_{i}\nabla
C_{i})+aj_{n,i}\mathbf{1}_{\mathbf{x}\in\mathrm{electrode}},$ (5)

where $D_{i}$ is the design variable dependent diffusivity, $a$ is the exposed
surface area per volume of the porous electrode, and
$\mathbf{1}_{\mathbf{x}\in\mathrm{electrode}}$ is an indicator function that
returns $1$ in the electrode and $0$ elsewhere. The reaction $j_{n,i}$ within
the homogenized electrode is given by
$j_{n,i}=k_{\mathrm{m}}(C_{i}^{\mathrm{s}}-C_{i})$, where $k_{\mathrm{m}}$ is
a mass transfer coefficient and $C_{i}^{\mathrm{s}}$ is the concentration at
the surface of the fibers of the electrode. At all boundaries, the no-flux
condition is applied for $C_{i}$.

The electrostatics of the solid and liquid are described by the following
Poisson equations,

$\displaystyle\nabla\cdot\left(-\sigma\nabla\Phi_{1}\right)$
$\displaystyle=-ai_{n}\mathbf{1}_{\mathbf{x}\in\mathrm{electrode}},$ (6)
$\displaystyle\nabla\cdot\left(-\kappa\nabla\Phi_{2}\right)$
$\displaystyle=ai_{n}\mathbf{1}_{\mathbf{x}\in\mathrm{electrode}},$ (7)

where $\Phi_{1}$ and $\Phi_{2}$ are electric potentials of the liquid and
solid, respectively, and $\sigma$ and $\kappa$ are the design variable
dependent conductivities of the liquid and solid, respectively. The exchange
current $i_{n}$ is given by the Butler-Volmer equation,

$i_{n}=\frac{i_{0}}{C^{\mathrm{ref}}}\left(C_{V^{2+}}^{s}e^{\beta\Delta\Phi}-C_{V^{3+}}^{s}e^{-\beta\Delta\Phi}\right),$
(8)

where $C^{\mathrm{ref}}=0.001\mathrm{M}$ is an arbitrary reference
concentration to maintain dimensional consistency, $\Delta\Phi$ is the local
overpotential, $\beta=0.5F/RT$, where $F$ is Faraday’s constant, $R$ is the
ideal gas constant, $T$ is the temperature, and a charge transfer coefficient
of $0.5$ has been assumed. The current $i_{n}$ is proportional to the mass
transfer reaction term as $i_{n}=Fj_{n,V^{2+}}=-Fj_{n,V^{3+}}$. A current is
discharged through the membrane via the liquid, so that
$-\kappa\nabla\Phi_{2}\cdot\mathbf{n}=I/A$ at the membrane, and the current
collector is attached to the solid part of the flow field so that $\Phi_{1}=0$
at the current collector; all other boundaries are no-flux boundaries for the
potential.

### Topology optimization

We follow the optimization method described by Beck et al. [66], with some
important differences. In contrast to [66], who performed optimization of a
porous electrode via a homogenized model, we utilize topology optimization to
generate flow field designs with well-defined solid and liquid boundaries
separated by a thin, smoothly varying interface [52]. Our goal is to solve,

$\displaystyle\min_{\gamma\in[0,1]}P$ (9) $\displaystyle\mathrm{s.t.}\,$
$\displaystyle\mathbf{R}_{\\{\mathbf{u},p,\Phi_{1},\Phi_{2},C_{i}\\}}=\mathbf{0},$
(10)

where $\mathbf{R}_{\\{\mathbf{u},p,\Phi_{1},\Phi_{2},C_{i}\\}}=\mathbf{0}$ are
the forward partial differential equations (PDEs) given in equations (3)-(7),
and $\gamma(\mathbf{x})$ is the design variable within the flow field that is
equal to $1$ where the flow field is pure liquid and equal to $0$ where the
flow field is pure solid.

Within the flow field, the design variable $\gamma(\mathbf{x})$ is used to
smoothly interpolate the material properties from that of pure solid to that
of pure liquid. Thus, the permeability, diffusivities, and conductivities in
equations (3)-(7) are all functions of the design variable. To ensure that
optimized designs are “black-white,” i.e. comprised of pure solid and pure
liquid with minimal material in between, the specific functional forms of the
interpolations are chosen carefully (e.g. see [53]) and attain the appropriate
values at the end points, $\gamma=0$ and $\gamma=1$. These functions are
described in the Supplementary Information.

To perform the optimization, the continuous adjoint approach is employed [74,
66]. Briefly, the Lagrangian of the problem is

$L=P+\int\mathbf{S}\cdot\mathbf{R}\,\rmd\mathbf{x}.$ (11)

The adjoint PDEs $\mathbf{S}$ corresponding to each PDE given in equations
(3)-(7) are derived analytically by setting the variation of $L$ with respect
to each dependent variable to 0, i.e. $\delta_{f_{i}}L=0$. From this, the
sensitivity is then determined by calculating $\delta L=\delta_{\gamma}L$.
These equations are solved in OpenFOAM, and the sensitivities are used within
the Method of Moving Asymptotes to update the values of $\gamma$ in each
computational cell [75]. For the simulations in this work, we use hexahedral
meshes with $4.3-12.6$M cells to solve the forward and adjoint PDEs. The
Helmholtz filter with a filter radius $r$ of $2\times$ the mesh size used in
the flow field, $r=167$ $\mu\mathrm{m}$, is used to regularize the solution,
and a Heaviside projection is used to push the solution to a black-white
design [76]. This process is described in more detail in the Supplementary
Information.

### Estimation of electrochemical parameters

Estimated parameters and correlations used in our electrochemical simulation
are shown in table 1. Notably, the Bruggeman correlation [71] is used to
describe properties within the porous electrode. Specifically, in the
electrode,

$\displaystyle D$ $\displaystyle=D_{0}\epsilon^{3/2},$ (12)
$\displaystyle\kappa$ $\displaystyle=\kappa_{0}\epsilon^{3/2},$ (13)
$\displaystyle\sigma$ $\displaystyle=\sigma_{0}(1-\epsilon)^{3/2},$ (14)

and in the flow field, interpolations relating the design variable to solid
and liquid properties are described in the Supplementary Information.

Table 1: Parameters used for electrochemical model

  * Description | Parameter | Value | Units | Ref.
---|---|---|---|---
Solid conductivity | $\kappa_{0}$ | $10^{4}$ | $\mathrm{S/m}$ | Est.
Electrode porosity | $\epsilon$ | $0.68$ | - | [71]
Electrode permeability | $\alpha_{E}$ | $5.53\cdot 10^{-11}$ | $\mathrm{m^{2}}$ | [71]
Electrode area/volume | $a$ | $8\cdot 10^{4}$ | $\mathrm{m^{-1}}$ | Est.
Electrolyte conductivity | $\sigma_{0}$ | $40$ | $\mathrm{S/m}$ | Est.
Electrolyte viscosity | $\mu$ | $8.9\cdot 10^{-4}$ | $\mathrm{Pa\cdot s}$ | Est.
Electrolyte density | $\rho$ | $10^{3}$ | $\mathrm{kg/m^{3}}$ | Est.
Diffusivity | $D_{0}$ | $2.4\cdot 10^{-10}$ | $\mathrm{m^{2}/s}$ | [77]
Mass transfer coeff. | $k_{m}$ | $4.4\cdot 10^{-4}|\mathbf{u}|^{0.4}$ | $\mathrm{m/s}$ | [78]
Exchange current | $i_{0}$ | $0.016$ | $\mathrm{A/m^{2}}$ | [71]
Equilibrium potential | $U_{0}$ | $-0.25$ | $\mathrm{V}$ | [67]

##  3 Results and discussion

### Optimized flow field designs

Figure 2: The fluid fraction of the optimized flow field, colored by fluid
pressure, is shown for various flow rates $Q$ and current densities $I/A$. The
porous electrode is shown as a translucent sheet. Fluid flows in through the
higher pressure channels and up into the porous electrode, reacts, and flows
back out through the lower pressure channels. An initial channel and land
width of $2$ $\mathrm{mm}$ is used. (A-C) $I/A=50$ $\mathrm{mA/cm^{2}}$, (D-F)
$I/A=100$ $\mathrm{mA/cm^{2}}$, (G-I) $I/A=150$ $\mathrm{mA/cm^{2}}$. Figure
3: Visualization of the fluid streamlines, colored by reactant concentration,
through the solid flow field for the example case where $I/A=100$
$\mathrm{mA/cm^{2}}$ and $Q=24$ $\mathrm{mL/min}$. Panel (A) shows a view of
the bottom half of the flow field that interfaces with the current collector.
Ramps appear at the ends of the channels to guide flow into and out of the
electrode. Panel (B) shows a view of the top half of the flow field that
interfaces with the membrane. Small secondary channels appear to evenly spread
reactant across the electrode. In both panels, the solid flow field is
slightly transparent and shaded to aid in visualization.

In figure 2, we show the liquid fraction of optimized flow fields for a range
of operating conditions where the optimization is initialized with channel and
land widths of 2 mm. In figures S1 and S2 of the Supplementary Information,
results are also shown where the optimization is initialized with channel and
land widths of 4 and 1.5 mm, respectively. The actual solid part is thus the
negative of what is shown in the figure. The porous electrode is represented
as a translucent sheet above the visible flow structures. As the operating
current and flow rate change, the dominant physical phenomena change. This, in
turn, changes the ratio of electrical power losses to pressure power losses,
leading to qualitatively different optimal flow field designs. When the flow
rate is low, the optimized fluid channels are thin and the flow field is
comprised mostly of solid material, suggesting that an interdigitated design
with thin channels would perform well. At even lower flow rates, we observe no
clear three-dimensional structures that form in the nearly fully-solid flow
field, and these designs are shown in figure S3 of the Supplementary
Information. In contrast, as the flow rate increases and more reactant is
supplied to the electrode, smaller secondary channels emerge from the primary
fluid channels beneath – resembling a tree-root structure – and guide fluid up
into the electrode. Additionally, the larger primary fluid distribution
channels grow larger, presumably to avoid large fluid resistances and
increased flow power losses.

The observed features of the optimized designs are fully three-dimensional,
providing evidence that there is room for improvement in more popular two-
dimensional designs. Remarkably, the optimization algorithm automatically
identified features which have previously been devised through
experimentation. At flow rates of $Q=$12 mL/min and 24 mL/min, the optimized
designs show a ramped channel design as in the work of Akuzum et al. [43],
helping to guide both unreacted reactant up into the electrode and used
reactant down out of the electrode. The optimized design, however, does
further taper the channel in addition to ramping. Additionally, at high flow
rates and currents, the optimization algorithm not only re-identified the use
of hierarchy, but it further suggests that the density of branching should
increase at lower operating currents and higher flow rates [41]. The
optimization procedure not only suggests interesting design candidates, but
because it uses a model and simulation for the system behavior to determine
the geometries, it also provides strategies to most effectively employ the 3D
features. To provide further intuition on the design features that form and
how fluid flows through the solid part, in figure 3, we show a visualization
of the fluid streamlines through the flow field for the example case where
$I/A=100$ $\mathrm{mA/cm^{2}}$ and $Q=24$ $\mathrm{mL/min}$. These streamlines
are colored by reactant concentration $C_{V^{2+}}$ to highlight the depletion
of concentration as fluid leaves the flow field. In figure 3(A), we observe
the previously described ramps that appear at the ends of the primary fluid
distribution channels, and in figure 3(B), we observe the smaller secondary
channels that form at the flow field-electrode interface that work to evenly
distribute fluid across the electrode.

In our optimization procedure, it is important to note that the length scale
of the smallest feature size can, in principle, be controlled to meet minimum
feature size constraints imposed by a specific additive manufacturing
technique. As described in §2, in these optimization simulations, a Helmholtz
filter [76] is used to regularize the flow field design, which as a result
sets the minimum feature size that can arise here. In this work, we have
chosen a filter radius that results in a minimum feature size of approximately
$0.5$ mm, however, this filter radius could be modified. Current additive
manufacturing techniques have resolutions that range from sub-micron to
millimeter scale [44].

### Power losses and efficiencies

Figure 4: (A) Electric and flow power losses of the unoptimized and optimized
flow fields for various flow rates $Q$, current densities $I/A$, and initial
conditions. Current densities $I/A=\\{50,100,150\\}$ $\mathrm{mA/cm^{2}}$,
flow rates $Q=\\{2.4,12,24,36\\}$ $\mathrm{mL/min}$, and flow fields
initialized with land and channel widths of $\\{1.5,2,4\\}$ $\mathrm{mm}$ are
considered. (B) Corresponding power efficiencies $\Xi=1-P/(IU_{0})$. Asterisks
indicate the flow rate yielding the highest power efficiency. In the
individual bars, the light portion of the bar shows the value before the
optimization procedure and the dark portion the value after, as indicated by
the arrows.

The qualitative design features observed in figure 2 are now quantitatively
supported. As implied by the variety of designs shown in figure 2, the optimal
flow field design is achieved by balancing the electrical and flow power
losses. This balance is similar to what is observed in the optimized
homogenized porous electrodes studied by Beck et al. [66], where lower
operating flow rates lead to less porous optimized structures and vice versa.
The associated electric and flow power losses of the optimized designs are
shown in figure 4(A), along with the power losses obtained with the
traditional, unoptimized interdigitated flow field. The power efficiencies
$\Xi$, calculated from equation (2), corresponding to these power losses are
shown in figure 4(B). In all cases, the optimized designs lead to greater
power efficiency than the interdigitated designs.

Figure 5: Optimized fluid distribution channels. Slices of the liquid fraction
of the flow field are shown at $y=1$ mm (green) and close to the flow field-
electrode interface ($y=2.9$ mm). Current densities $I/A=\\{50,100,150\\}$
$\mathrm{mA/cm^{2}}$ are considered, and flow rates $Q$ are chosen to
correspond to the optimal power efficiency as shown in figure 4(B). The
different panels show different starting interdigitated designs; specifically
the initial land and channels widths are (A,D,G) 1.5 mm (B,E,H) 2 mm and
(C,F,I) 4mm. As the initial condition changes, the optimization procedure
leads to qualitatively different designs: secondary channels that are largely
perpendicular to the primary flow distribution channels appear when the
initial channel and land widths are thick, whereas they are parallel when the
initial channel and land widths are thin.

At a particular current, traditional flow field designs yield power losses
where the electric losses are maximized when the supplied flow rate is
minimized, and the flow losses are maximized when the supplied flow rate is
maximized. Thus, the optimization prioritizes reduction of the appropriate
power loss contribution based on the operating conditions. This is evident in
figure 4(A) from the large reduction in the electric power loss for designs
generated for low flow rates, and, conversely, in the large reduction in flow
power loss for designs generated for high flow rates. For the particular
system studied here, the electrical power losses typically dominate, with the
flow power losses typically contributing $<10\%$, due to the fact that the
dominant source of flow resistance is only associated with the thin porous
electrode. For larger area electrodes, the flow losses would be larger, and
this optimization process would appropriately reprioritize the loss reduction.
Since the specific choice of channel and land widths of traditional
interdigitated flow fields are known to affect performance, we also consider
simulations initialized with channel and land widths of 1.5 and 4 mm, and
observe the same trends for each initial condition.

For a given current, the power efficiency of the optimized design reaches
approximately the same value across all of the flow rates and initial
conditions tested. For example, when $I/A=100$ $\mathrm{mA/cm^{2}}$, the
optimized power efficiency $\Xi$ is always $\approx 0.82$. In contrast, at the
same current density, the standard interdigitated designs have efficiency
$\Xi\approx 0.25$, when $Q=2.4$ $\mathrm{mL/min}$ and the channel and land
widths are 4mm, while $\Xi\approx 0.78$, when $Q=36$ $\mathrm{mL/min}$ and the
channel and land widths are 1.5 mm. The improvement from optimization at this
current density thus ranges from $5-228\%$. This demonstrates that topology
optimization can be an effective tool to identify a high performance geometry
regardless of the desired operating conditions. Alternatively, a single
standard unoptimized design shows strong variations in performance across
operating conditions. While it is certainly possible to use a traditional
interdigitated design to yield high performance, a significant amount of
testing is required to find the optimal fluid channel and solid land widths
for the specific system used; this process needs to be repeated whenever these
conditions change. For example, Gerhardt et al. [33] find a channel and land
width that yields a maximum pumping-corrected voltage efficiency, a metric
similar to our definition of $\Xi$, of $\approx 80\%$ for their ferrocyanide-
ferricyanide flow battery system operating at $250$ $\mathrm{mA/cm^{2}}$, but
also discuss how the optimal dimensions change as the operating conditions
change. A key takeaway from the present work is that this iterative process
can be automated through the use of topology optimization methods such as the
one presented here.

### Analysis of 3D flow field design features

All of the flow field designs presented in figure 2 have been obtained by
conducting simulations with the same initial condition: an interdigitated flow
field with channel and land width of $2$ $\mathrm{mm}$. However, as evidenced
by the power efficiencies presented in figure 4(B), the optimization scheme is
able to yield a design with similar performance, regardless of the initial
channel and land widths chosen, even though the ultimate design depends on the
specific initial conditions. To further understand how the optimization
procedure is selecting design features, in figure 5, we show slices of the
fluid fraction of the optimized flow fields from figures 2, S1, and S2, where
the flow rate $Q$ is selected based on the optimal power efficiency as
observed in figure 4(B). In doing so, we highlight how the optimized three-
dimensional structures change for various currents and starting initial
conditions. In green, the primary underlying fluid distribution channels are
shown by plotting the fluid fraction along a slice within the flow field
($y=1$ mm), and in blue and red, the secondary fluid channels are shown by
plotting the fluid fraction along a slice near the flow field-electrode
interface ($y=2.9$ mm). The fresh reactant-laden fluid thus enters the flow
field through the green channels on the bottom, enters the electrode in the
out-of-page $y$ direction through the red high pressure regions, exits the
electrode through the blue low pressure regions, and out of the device through
the green channels on the top.

The structure of the secondary channels that form provide clues for what
constitutes an optimal design. In figure 5(C), where the current is relatively
low and the initial interdigitated design has relatively thick channel and
land widths, secondary channels appear perpendicular to the primary flow
channels. In contrast, in figure 5(G), where the opposite is true, secondary
channels appear parallel and coincident to the primary flow channels. In other
cases, a mixed regime appears, where both parallel and perpendicular secondary
channels appear. This provides evidence that the initial channel and land
widths in figure 5(C) are too thick for the system, and the optimization
algorithm works to create secondary channels to more evenly distribute
reactant. Similarly, this shows that the initial channel and land widths in
figure 5(G) are close to optimal, and the interdigitated design under these
operating conditions already performs relatively well. In figure 5(A), we
observe a design resembling that of figure 5(G), but with small perpendicular
channels appearing: since the operating current in this case is lower, there
is less need for the optimization to prioritize electrical contact with the
flow field.

To further understand and quantify how evenly the reactant is distributed, we
consider the dimensionless root-mean-squared-deviation of the oxidant
concentration within the electrode,

$\mathrm{RMSD}(C_{V^{3+}})=\sqrt{\left<\left(\frac{C_{V^{3+}}}{C_{\mathrm{inlet}}}-1\right)^{2}\right>},$
(15)

where $\left<\cdot\right>$ denotes an average over the volume of the porous
electrode. This quantity thus gives a measure of the concentration
polarization within the electrode. Similarly, we also consider the
dimensionless root-mean-squared-deviation of the exchange current $i_{n}$
within the electrode,

$\mathrm{RMSD}(i_{n})=\sqrt{\left<\left(\frac{ai_{n}V}{I}-1\right)^{2}\right>},$
(16)

where $a$ is the exposed surface area per volume of the electrode, $V$ is the
electrode volume, and $i_{n}$ is the exchange current given by the Butler-
Volmer relation in equation (8). It is evident that the lower these quantities
are, the more uniform the distribution of reactant and Faradaic reaction. In
the limit of perfect uniformity, both approach a value of zero. In figure 6(A)
and (B), we show both $\mathrm{RMSD}(C_{V^{3+}})$ and $\mathrm{RMSD}(i_{n})$
as a function of flow rate for various currents and initial channel and land
widths. In all cases, we observe an appreciable reduction in these quantities
for the optimized designs, with the most significant improvement observed at
low flow rate.

While both of these quantities give a global measure of how well the reactant
is distributed within the electrode, it is also insightful to examine the
local quantities. In figure 6(C)-(H), we compare the vertical fluid velocity,
local concentration polarization, and local exchange current deviation within
a slice of the electrode between the unoptimized and optimized flow field; we
consider the particular case where $I/A=100$ $\mathrm{mA/cm^{2}}$, $Q=24$
$\mathrm{mL/min}$, and the initial channel and land width is $2$
$\mathrm{mm}$. This highlights how the optimized design is able to more evenly
utilize the entire electrode and avoid electrode dead zones. In figure 6(C),
we observe that there is a large area of the electrode where the fluid is
essentially stagnant. In contrast, the optimized design is able to deliver
fluid to and from the electrode nearly uniformly (figure 6(F)), thus reducing
the concentration and exchange current polarization (figure 6(G,H)).

Figure 6: Dead zones of the electrode. Panels (A) and (B) show the root-mean-
squared deviation of the species concentration and the exchange current,
respectively, within the electrode for both the unoptimized and optimized flow
fields. Panels (C-H) show the local vertical velocity, concentration
polarization and exchange current polarization at the center of the porous
electrode for the specific case of $I/A=100$ $\mathrm{mA/cm^{2}}$, $Q=24$
$\mathrm{mL/min}$, and an initial channel and land width of $2$ $\mathrm{mm}$.
The images to the left show the solid fraction of the corresponding flow
field. The optimized flow field more evenly distributes reactant, leading to
more uniform reaction and a reduction in the extent of dead zones in the
electrode. In the color legend, numbers on the top correspond to panels (C-E)
and numbers on the bottom correspond to panels (F-H).

##  4 Conclusions

In this work, we have shown that topology optimization, coupled with a model
of the fluid mechanics, mass transfer, and electrostatics describing a flow
battery, can be used to improve the performance of a flow field. Specifically,
we use topology optimization to design flow fields with full three-dimensional
variation, i.e. 3D flow fields. We have improved upon the popular two-
dimensional interdigitated flow field design and identified pertinent three-
dimensional features, such as ramps and hierarchy, that have previously been
devised heuristically. Additionally, we have shown that our optimized designs
lead to better overall reactant and reaction distribution within the
electrode. Notably, while using a traditional interdigitated design can
certainly yield high performance, this requires tedious testing to find the
optimal land and channel widths for the particular system; our method has the
potential to automate that process. Understanding and improving fluid and mass
distribution in electrochemical systems such as flow batteries is critical,
and this work gives preliminary evidence that investing effort into developing
advanced manufacturing techniques to fabricate such complex components may be
fruitful. Our ongoing work includes experimentally demonstrating that the flow
fields like the ones described in the present work can provide real-world
benefit.

There are a number of other considerations that may need to be taken into
account in an industrial setting for a specific application. For example, one
may be interested in maximizing the single-pass conversion efficiency of the
reactant, instead of the power loss, when considering a system-level analysis
of a flow battery stack. Next, while we have focused on a particular flow
battery size in this work, it would be of interest to explore how the benefits
of optimization are affected as the device is scaled up. Certainly, the ratio
of electrical losses to flow losses would be expected to change with system
size. Furthermore, it would be of interest to understand how the optimized
flow fields in this work perform within a full cell, where the membrane would
be expected to make a nontrivial contribution to the performance. In a real-
world application, it would also be important to understand whether or not
these optimized designs would perform well in the presence of obstacles such
as bubbles or unwanted debris that can clog flow channels. Finally, while flow
batteries have been subject of the present work, the distribution of reactants
remains to be a key problem in a range of engineered systems, such as
$\mathrm{CO_{2}}$ [79] and water electrolyzers, fuel cells, and bioreactors.
Given a particular system chemistry and desired quantity to be minimized, the
process outlined in this work can be repeated, leading to the potential to
improve mass distribution in such systems.

T.Y.L. thanks Drs. A. Wong, J. Davis, B. S. Jayathilake, and S. Chandrasekaran
for insightful conversations about interdigitated flow fields and additive
manufacturing. This work was performed under the auspices of the U.S.
Department of Energy by the Lawrence Livermore National Laboratory under
contract DE-AC52-07NA27344. This work was also supported by the Lawrence
Livermore National Laboratory LDRD 19-SI-005 and 19-ERD-035. This work was
performed under the auspices of a cooperative research and development
agreement between the Lawrence Livermore National Laboratory, Stanford
University, and TOTAL American Services, Inc. (affiliate of TOTAL SE) under
agreement number TC02307. LLNL release number: LLNL-JRNL-832021-DRAFT.

## References

## References

  * [1] Steven Chu and Arun Majumdar.  Opportunities and challenges for a sustainable energy future.  nature, 488(7411):294–303, 2012.
  * [2] Steven Chu, Yi Cui, and Nian Liu.  The path towards sustainable energy.  Nature materials, 16(1):16–22, 2017.
  * [3] Turgut M Gür.  Review of electrical energy storage technologies, materials and systems: challenges and prospects for large-scale grid storage.  Energy Environ. Sci., 11(10):2696–2767, 2018.
  * [4] Anya Castillo and Dennice F Gayme.  Grid-scale energy storage applications in renewable energy integration: A survey.  Energy Conversion and Management, 87:885–894, 2014.
  * [5] Erik Ela, Michael Milligan, and Brendan Kirby.  Operating reserves and variable generation.  Technical report, National Renewable Energy Lab.(NREL), Golden, CO (United States), 2011.
  * [6] Haisheng Chen, Thang Ngoc Cong, Wei Yang, Chunqing Tan, Yongliang Li, and Yulong Ding.  Progress in electrical energy storage system: A critical review.  Progress in natural science, 19(3):291–312, 2009.
  * [7] Gareth Kear, Akeel A Shah, and Frank C Walsh.  Development of the all-vanadium redox flow battery for energy storage: a review of technological, financial and policy aspects.  International journal of energy research, 36(11):1105–1120, 2012\.
  * [8] Adam Z Weber, Matthew M Mench, Jeremy P Meyers, Philip N Ross, Jeffrey T Gostick, and Qinghua Liu.  Redox flow batteries: a review.  Journal of applied electrochemistry, 41(10):1137–1164, 2011.
  * [9] Grigorii L Soloveichik.  Flow batteries: current status and trends.  Chemical reviews, 115(20):11533–11558, 2015.
  * [10] Wei Wang, Qingtao Luo, Bin Li, Xiaoliang Wei, Liyu Li, and Zhenguo Yang.  Recent progress in redox flow battery research and development.  Advanced Functional Materials, 23(8):970–986, 2013.
  * [11] Mike L Perry and Adam Z Weber.  Advanced redox-flow batteries: a perspective.  Journal of The Electrochemical Society, 163(1):A5064, 2015.
  * [12] Gang Qiu, Abhijit S Joshi, CR Dennison, KW Knehr, EC Kumbur, and Ying Sun.  3-d pore-scale resolved model for coupled species/charge/fluid transport in a vanadium redox flow battery.  Electrochimica Acta, 64:46–64, 2012.
  * [13] E Knudsen, P Albertus, KT Cho, AZ Weber, and A Kojic.  Flow simulation and analysis of high-power flow batteries.  Journal of Power Sources, 299:617–628, 2015.
  * [14] Xinyou Ke, Joseph M Prahl, J Iwan D Alexander, Jesse S Wainright, Thomas A Zawodzinski, and Robert F Savinell.  Rechargeable redox flow batteries: flow fields, stacks and design considerations.  Chemical Society Reviews, 47(23):8721–8743, 2018.
  * [15] V Pavan Nemani and Kyle C Smith.  Uncovering the role of flow rate in redox-active polymer flow batteries: simulation of reaction distributions with simultaneous mixing in tanks.  Electrochimica Acta, 247:475–485, 2017.
  * [16] Purna C Ghimire, Arjun Bhattarai, Rüdiger Schweiss, Günther G Scherer, Nyunt Wai, and Qingyu Yan.  A comprehensive study of electrode compression effects in all vanadium redox flow batteries including locally resolved measurements.  Applied Energy, 230:974–982, 2018.
  * [17] Andrew A Wong, Shmuel M Rubinstein, and Michael J Aziz.  Direct visualization of electrochemical reactions and heterogeneous transport within porous electrodes in operando by fluorescence microscopy.  Cell Reports Physical Science, 2(4):100388, 2021.
  * [18] Luis F Arenas, Carlos Ponce de León, and Frank C Walsh.  Redox flow batteries for energy storage: their promise, achievements and challenges.  Current Opinion in Electrochemistry, 16:117–126, 2019.
  * [19] S Kumar and S Jayanti.  Effect of flow field on the performance of an all-vanadium redox flow battery.  Journal of Power Sources, 307:782–787, 2016.
  * [20] Mirko Messaggi, Patrizio Canzi, R Mereu, A Baricci, F Inzoli, A Casalegno, and M Zago.  Analysis of flow field design on vanadium redox flow battery performance: Development of 3d computational fluid dynamic model and experimental validation.  Applied energy, 228:1057–1070, 2018.
  * [21] Sandip Maurya, Phong Thanh Nguyen, Yu Seung Kim, Qinjun Kang, and Rangachary Mukundan.  Effect of flow field geometry on operating current density, capacity and performance of vanadium redox flow battery.  Journal of Power Sources, 404:20–27, 2018.
  * [22] T Jyothi Latha and S Jayanti.  Hydrodynamic analysis of flow fields for redox flow battery applications.  Journal of Applied Electrochemistry, 44(9):995–1006, 2014.
  * [23] Jacob Houser, Jason Clement, Alan Pezeshki, and Matthew M Mench.  Influence of architecture and material properties on vanadium redox flow battery performance.  Journal of Power Sources, 302:369–377, 2016.
  * [24] Robert M Darling and Mike L Perry.  The influence of electrode and channel configurations on flow battery performance.  Journal of The Electrochemical Society, 161(9):A1381, 2014.
  * [25] Xianguo Li and Imran Sabir.  Review of bipolar plates in pem fuel cells: Flow-field designs.  International journal of hydrogen energy, 30(4):359–371, 2005.
  * [26] Qian Xu, TS Zhao, and Cheng Zhang.  Performance of a vanadium redox flow battery with and without flow fields.  Electrochimica Acta, 142:61–67, 2014.
  * [27] Ravendra Gundlapalli and Sreenivas Jayanti.  Effect of channel dimensions of serpentine flow fields on the performance of a vanadium redox flow battery.  Journal of Energy Storage, 23:148–158, 2019.
  * [28] Qian Xu, TS Zhao, and PK Leung.  Numerical investigations of flow field designs for vanadium redox flow batteries.  Applied energy, 105:47–56, 2013.
  * [29] Xinyou Ke, J Iwan D Alexander, Joseph M Prahl, and Robert F Savinell.  Flow distribution and maximum current density studies in redox flow batteries with a single passage of the serpentine flow channel.  Journal of Power Sources, 270:646–657, 2014.
  * [30] Xinyou Ke, J Iwan D Alexander, Joseph M Prahl, and Robert F Savinell.  A simple analytical model of coupled single flow channel over porous electrode in vanadium redox flow battery with serpentine flow channel.  Journal of Power Sources, 288:308–313, 2015.
  * [31] Xinyou Ke, Joseph M Prahl, J Iwan D Alexander, and Robert F Savinell.  Redox flow batteries with serpentine flow fields: Distributions of electrolyte flow reactant penetration into the porous carbon electrodes and effects on performance.  Journal of Power Sources, 384:295–302, 2018.
  * [32] Xinyou Ke, Joseph M Prahl, J Iwan D Alexander, and Robert F Savinell.  Mathematical modeling of electrolyte flow in a segment of flow channel over porous electrode layered system in vanadium flow battery with flow field design.  Electrochimica Acta, 223:124–134, 2017.
  * [33] Michael R Gerhardt, Andrew A Wong, and Michael J Aziz.  The effect of interdigitated channel and land dimensions on flow cell performance.  Journal of The Electrochemical Society, 165(11):A2625, 2018.
  * [34] Cong Yin, Yan Gao, Shaoyun Guo, and Hao Tang.  A coupled three dimensional model of vanadium redox flow battery for flow field designs.  Energy, 74:886–895, 2014.
  * [35] Shohji Tsushima and Takahiro Suzuki.  Modeling and simulation of vanadium redox flow battery with interdigitated flow field for optimizing electrode architecture.  Journal of The Electrochemical Society, 167(2):020553, 2020.
  * [36] Cong Yin, Yan Gao, Guangyou Xie, Ting Li, and Hao Tang.  Three dimensional multi-physical modeling study of interdigitated flow field in porous electrode for vanadium redox flow battery.  Journal of Power Sources, 438:227023, 2019.
  * [37] Ravendra Gundlapalli and Sreenivas Jayanti.  Performance characteristics of several variants of interdigitated flow fields for flow battery applications.  Journal of Power Sources, 467:228225, 2020.
  * [38] Xinyou Ke, Joseph M Prahl, J. Iwan D Alexander, Jesse S Wainright, Thomas A Zawodzinski, and Robert F Savinell.  Rechargeable redox flow batteries: flow fields, stacks and design considerations.  Chem. Soc. Rev., 47(23):8721–8743, 2018.
  * [39] Oladapo Christopher Esan, Xingyi Shi, Zhefei Pan, Xiaoyu Huo, Liang An, and T S Zhao.  Modeling and Simulation of Flow Batteries.  Adv. Energy Mater., 22:2000758–42, June 2020.
  * [40] Kleber Marques Lisboa, Julian Marschewski, Neil Ebejer, Patrick Ruch, Renato Machado Cotta, Bruno Michel, and Dimos Poulikakos.  Mass transport enhancement in redox flow batteries with corrugated fluidic networks.  Journal of Power Sources, 359:322–331, 2017.
  * [41] Yikai Zeng, Fenghao Li, Fei Lu, Xuelong Zhou, Yanping Yuan, Xiaoling Cao, and Bo Xiang.  A hierarchical interdigitated flow field design for scale-up of high-performance redox flow batteries.  Applied Energy, 238:435–441, 2019.
  * [42] Ravendra Gundlapalli and Sreenivas Jayanti.  Effective splitting of serpentine flow field for applications in large-scale flow batteries.  Journal of Power Sources, 487:229409, 2021.
  * [43] Bilen Akuzum, Yigit Can Alparslan, Nicholas C Robinson, Ertan Agar, and E Caglan Kumbur.  Obstructed flow field designs for improved performance in vanadium redox flow batteries.  Journal of Applied Electrochemistry, 49(6):551–561, 2019.
  * [44] Adriano Ambrosi, Raymond Rong Sheng Shi, and Richard D Webster.  3d-printing for electrolytic processes and electrochemical flow systems.  Journal of Materials Chemistry A, 8(42):21902–21929, 2020.
  * [45] Greig Chisholm, Philip J Kitson, Niall D Kirkaldy, Leanne G Bloor, and Leroy Cronin.  3d printed flow plates for the electrolysis of water: an economic and adaptable approach to device manufacture.  Energy & Environmental Science, 7(9):3026–3032, 2014.
  * [46] Jesse R Hudkins, Danika G Wheeler, Bruno Peña, and Curtis P Berlinguette.  Rapid prototyping of electrolyzer flow field plates.  Energy & Environmental Science, 9(11):3417–3423, 2016.
  * [47] Gaoqiang Yang, Jingke Mo, Zhenye Kang, Yeshi Dohrmann, Frederick A List III, Johney B Green Jr, Sudarsanam S Babu, and Feng-Yuan Zhang.  Fully printed and integrated electrolyzer cells with additive manufacturing for high-efficiency water splitting.  Applied Energy, 215:202–210, 2018.
  * [48] Gaoqiang Yang, Shule Yu, Jingke Mo, Zhenye Kang, Yeshi Dohrmann, Frederick A List III, Johney B Green Jr, Sudarsanam S Babu, and Feng-Yuan Zhang.  Bipolar plate development with additive manufacturing and protective coating for durable and high-efficiency hydrogen production.  Journal of Power Sources, 396:590–598, 2018.
  * [49] Justin C Bui, Jonathan T Davis, and Daniel V Esposito.  3d-printed electrodes for membraneless water electrolysis.  Sustainable Energy & Fuels, 4(1):213–225, 2020.
  * [50] Victor A Beck, Anna N Ivanovskaya, Swetha Chandrasekaran, Jean-Baptiste Forien, Sarah E Baker, Eric B Duoss, and Marcus A Worsley.  Inertially enhanced mass transport using 3d-printed porous flow-through electrodes with periodic lattice structures.  Proceedings of the National Academy of Sciences, 118(32), 2021.
  * [51] Toshihiko Yoshida and Koichi Kojima.  Toyota mirai fuel cell vehicle and progress toward a future hydrogen society.  The Electrochemical Society Interface, 24(2):45, 2015.
  * [52] Ole Sigmund and Kurt Maute.  Topology optimization approaches.  Structural and Multidisciplinary Optimization, 48(6):1031–1055, 2013.
  * [53] Sumer B Dilgen, Cetin B Dilgen, David R Fuhrman, Ole Sigmund, and Boyan S Lazarov.  Density based topology optimization of turbulent flow heat transfer systems.  Structural and Multidisciplinary Optimization, 57(5):1905–1918, 2018.
  * [54] Hiroki Kobayashi, Kentaro Yaji, Shintaro Yamasaki, and Kikuo Fujita.  Freeform winglet design of fin-and-tube heat exchangers guided by topology optimization.  Applied Thermal Engineering, 161:114020, 2019.
  * [55] Lukas Christian Høghøj, Daniel Ruberg Nørhave, Joe Alexandersen, Ole Sigmund, and Casper Schousboe Andreasen.  Topology optimization of two fluid heat exchangers.  International Journal of Heat and Mass Transfer, 163:120543, 2020\.
  * [56] Florian Feppon, Grégoire Allaire, Charles Dapogny, and Pierre Jolivet.  Body-fitted topology optimization of 2d and 3d fluid-to-fluid heat exchangers.  Computer Methods in Applied Mechanics and Engineering, 376:113638, 2021.
  * [57] Miguel A Salazar De Troya, Daniel A Tortorelli, and Victor A Beck.  Two dimensional topology optimization of heat exchangers with the density and level-set methods.  In 14th WCCM-ECCOMAS Congress, volume 1300. Lawrence Livermore National Lab.(LLNL), Livermore, CA (United States), 2021.
  * [58] Miguel A Salazar de Troya, Daniel A Tortorelli, Julian Andrej, and Victor A Beck.  Three-dimensional topology optimization of heat exchangers with the level-set method.  arXiv preprint arXiv:2111.09471, 2021.
  * [59] Martin P Bendsøe.  Optimal shape design as a material distribution problem.  Structural optimization, 1(4):193–202, 1989.
  * [60] Tyler E Bruns and Daniel A Tortorelli.  Topology optimization of non-linear elastic structures and compliant mechanisms.  Computer methods in applied mechanics and engineering, 190(26-27):3443–3459, 2001.
  * [61] Allan Gersborg-Hansen, Ole Sigmund, and Robert B Haber.  Topology optimization of channel flow problems.  Structural and multidisciplinary optimization, 30(3):181–192, 2005\.
  * [62] Thomas Borrvall and Joakim Petersson.  Topology optimization of fluids in stokes flow.  International journal for numerical methods in fluids, 41(1):77–107, 2003.
  * [63] Kentaro Yaji, Shintaro Yamasaki, Shohji Tsushima, Takahiro Suzuki, and Kikuo Fujita.  Topology optimization for the design of flow fields in a redox flow battery.  Structural and multidisciplinary optimization, 57(2):535–546, 2018\.
  * [64] Chih-Hsiang Chen, Kentaro Yaji, Shintaro Yamasaki, Shohji Tsushima, and Kikuo Fujita.  Computational design of flow fields for vanadium redox flow batteries via topology optimization.  Journal of Energy Storage, 26:100990, 2019.
  * [65] Reza Behrou, Alberto Pizzolato, and Antoni Forner-Cuenca.  Topology optimization as a powerful tool to design advanced pemfcs flow fields.  International Journal of Heat and Mass Transfer, 135:72–92, 2019\.
  * [66] Victor A Beck, Jonathan J Wong, Charles F Jekel, Daniel A Tortorelli, Sarah E Baker, Eric B Duoss, and Marcus A Worsley.  Computational design of microarchitected porous electrodes for redox flow batteries.  Journal of Power Sources, 512:230453, 2021.
  * [67] Liyu Li, Soowhan Kim, Wei Wang, M Vijayakumar, Zimin Nie, Baowei Chen, Jianlu Zhang, Guanguang Xia, Jianzhi Hu, Gordon Graff, et al.  A stable vanadium redox-flow battery with high energy density for large-scale energy storage.  Advanced Energy Materials, 1(3):394–400, 2011.
  * [68] Qiong Zheng, Xianfeng Li, Yuanhui Cheng, Guiling Ning, Feng Xing, and Huamin Zhang.  Development and perspective in vanadium flow battery modeling.  Applied energy, 132:254–266, 2014.
  * [69] Trevor J Davies and Joseph J Tummino.  High-performance vanadium redox flow batteries with graphite felt electrodes.  C, 4(1):8, 2018.
  * [70] Yun Wang and Sung Chan Cho.  Analysis and three-dimensional modeling of vanadium flow batteries.  Journal of The Electrochemical Society, 161(9):A1200, 2014.
  * [71] AA Shah, MJ Watt-Smith, and FC Walsh.  A dynamic performance model for redox-flow batteries involving soluble species.  Electrochimica Acta, 53(27):8087–8100, 2008.
  * [72] Xiangkun Ma, Huamin Zhang, and Feng Xing.  A three-dimensional model for negative half cell of the vanadium redox flow battery.  Electrochimica Acta, 58:238–246, 2011.
  * [73] John Newman and William Tiedemann.  Porous-electrode theory with battery applications.  AIChE Journal, 21(1):25–41, 1975.
  * [74] Carsten Othmer.  A continuous adjoint formulation for the computation of topological and surface sensitivities of ducted flows.  International journal for numerical methods in fluids, 58(8):861–877, 2008.
  * [75] K Svandberg.  The method of moving asymptotes - a new method for structural optimization.  International Journal for Numerical Methods in Engineering, 24:359–373, 1987.
  * [76] Boyan Stefanov Lazarov and Ole Sigmund.  Filters in topology optimization based on helmholtz-type differential equations.  International Journal for Numerical Methods in Engineering, 86(6):765–781, 2011.
  * [77] Tomoo Yamamura, Nobutaka Watanabe, Takashi Yano, and Yoshinobu Shiokawa.  Electron-transfer kinetics of np3+/ np4+, npo2+/ npo2 2+, v2+/ v3+, and vo2+/ vo2+ at carbon electrodes.  Journal of The Electrochemical Society, 152(4):A830, 2005.
  * [78] D Schmal, J Van Erkel, and PJ Van Duin.  Mass transfer at carbon fibre electrodes.  Journal of applied electrochemistry, 16(3):422–430, 1986.
  * [79] Danika G Wheeler, Benjamin AW Mowbray, Angelica Reyes, Faezeh Habibzadeh, Jingfu He, and Curtis P Berlinguette.  Quantification of water transport in a co 2 electrolyzer.  Energy & Environmental Science, 13(12):5126–5134, 2020.