Topology optimization of 3D flow fields for flow batteries

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

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

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.

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,

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,

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,

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,

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.

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,

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

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.

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,

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

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].

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.

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.

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,

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,

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)).