Simulating Discharge Curves of an All-Aqueous TRAB to Identify Pathways for Improving System Performance

We modified a previous steady-state model to simulate time-dependent discharge curves of the Cu_aq-TRAB battery using room temperature parameters from the previous model (Table S1) and compiled the equations with COMSOL Multiphysics. ¹⁷ A mesh independent study was conducted and found that there was a 0.5% change in simulated energy output as the number of mesh elements was increased from 13,000 to 16,000, thus the mesh with 16,000 elements was used for this study (Fig. S1). The model used laminar, incompressible Brinkman and continuity equations for solving the flow field,

$$\begin{eqnarray}&&\begin{array}{c}\displaystyle \frac{\rho }{{\epsilon }}\left(\displaystyle \frac{\partial \vec{u}}{\partial t}+\left(\vec{u}\,\cdot {\rm{\nabla }}\right)\mathop{\displaystyle \frac{u}{{\epsilon }}}\limits^{\longrightarrow}\right)\\ =\,-\vec{{\rm{\nabla }}}p+{\rm{\nabla }}\cdot \left[\displaystyle \frac{1}{{\epsilon }}\left\{\mu \left({\rm{\nabla }}\vec{u}+\left({\rm{\nabla }}{\vec{u}}^{T}\right)\right)-\displaystyle \frac{2}{3}\mu \left({\rm{\nabla }}\cdot \vec{u}\right){\rm{I}}\right\}\right]\,\end{array}\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}&&\rho {\rm{\nabla }}\cdot \vec{u}={Q}_{m}\end{eqnarray} \tag{ 2 }$$

where ρ is the fluid density, ε is the electrode porosity, $\vec{u}$ is the velocity vector, t is time, p is the fluid pressure, μ is the fluid viscosity, and Q_m is a source or sink term. ²³ The flow path of the electrolyte was simulated independently of the electrochemical and diffusive transport equations because it was assumed that the change in concentration would not change the flow path. Both electrodes and electrolytes were governed by electroneutrality conditions, conservation of mass, and the convection–diffusion equations,

$$\begin{eqnarray}&&\displaystyle \sum {z}_{i}{c}_{i}=0,\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial {c}_{i}}{\partial t}+\vec{{\rm{\nabla }}}\cdot {\vec{N}}_{i}=0,\end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray}&&{\vec{N}}_{i}=-{D}_{i}\vec{{\rm{\nabla }}}{c}_{i}-\displaystyle \frac{{z}_{i}F}{RT}{D}_{i}{c}_{i}\vec{{\rm{\nabla }}}{\phi }_{e}+\vec{u}{c}_{i},\end{eqnarray} \tag{ 5 }$$

where z_i is charge of the ith ion, c_i is the bulk molar concentration of the ith species, $\vec{N}$ _i is the flux and D_i is the diffusion coefficient of an individual species, F is Faraday's constant, R is the molar gas constant, T was the temperature of the fluid, and ${\phi }_{e}$ was the potential of the electrolyte. Electrode reaction kinetics were modelled using the linearized Butler-Volmer equation,

$$\begin{eqnarray}&&{i}_{loc}={i}_{0}\left(\displaystyle \frac{\left({\alpha }_{a}+{\alpha }_{c}\right)F}{RT}\right)\eta \end{eqnarray} \tag{ 6 }$$

where i_loc is the local current density, i₀ is the exchange current density of the reaction, α_a and α_c are the anodic and cathodic transfer coefficients, and η is the overpotential.

Unique to this new work, the model was constructed in 2D with the simulated flow field being a projection of the three-dimensional flow field (Fig. S2). The membrane was modelled as a fixed ohmic drop and acted as a source/sink for the background electrolyte ions (ammonium and bromide) to conserve electroneutrality. To make the model time-dependent and simulate a recirculating electrolyte, we added an ordinary differential equation to calculate the single pass conversion in the flow cell to change the bulk reservoir concentration and update the inlet concentration of all species in both electrolytes with time. The equation was defined using a surface integral as

$$\begin{eqnarray}&&\displaystyle \frac{V}{L}\displaystyle \frac{d{c}_{i,in}}{dt}=\displaystyle \int \displaystyle {\int }_{{S}_{outlet}}\left({N}_{i}\,\cdot n\right)dS-\displaystyle \int \displaystyle {\int }_{{S}_{inlet}}\left({N}_{i}\,\cdot n\right)dS\end{eqnarray} \tag{ 7 }$$

where V is the volume of the tank (50 mL), L is the length of the electrode (5 cm), c_i is the concentration of the ith species, N_i is the molar flux of the ith species, n is the normal vector from the surface, and S is the surface of interest (either the inlet or outlet flow field boundary). The initial copper concentration was 0.5 M to match the experimental setup used for model validation. ¹¹ The molar flux at each inlet and outlet was extrapolated to the full 3D distance using the out-of-plane thickness of the experimental flow field (5 cm).

To account for potential losses from ammonia in the negolyte crossing over to the posolyte, we modified the Nernst equation of the positive electrode reaction by adding an empirical decay rate, E⁰ _pos,decay [V s⁻¹] to the positive electrode standard potential:

$$\begin{eqnarray}&&{E}_{eq}=\left({E}_{pos}^{0}-{E}_{pos,decay}^{0}t\right)-\displaystyle \frac{RT}{nF}\,\mathrm{ln}\left(Q\right)\end{eqnarray} \tag{ 8 }$$

where E_eq is the equilibrium potential, E⁰ _pos is the positive electrode standard potential, t is time, R is the ideal gas constant, n is the number of electrons transferred, F is Faraday's constant, and Q is the reaction quotient. This decay rate was used as a fitting parameter to replicate experimental discharge curves and account for the self-discharge that occurs when free ammonia in the posolyte complexes with copper species and reduces the standard potential of the Cu(I/II) reaction. ^9,24 By utilizing this decay rate instead of constantly simulating the reaction and diffusion of all the species within the electrolyte, the model becomes significantly simpler while still effectively representing the influence of ammonia crossover, a crucial factor in reproducing Cu_aq-TRAB discharge curves. The use of this single, linear decay term can be justified by previous experiments that observed linear changes in ammonia concentration during membrane crossover experiments (see Fig. S3 of Ref. 11).

The energy efficiency of converting waste heat into electrical energy, η, was calculated using a simplified method relating the electrical energy output, $\hat{U}$ to the thermal energy input, E_heat (107 Wh l^–1), ^5,17

$$\begin{eqnarray}&&\eta =\displaystyle \frac{\hat{U}}{{E}_{heat}}\end{eqnarray} \tag{ 9 }$$

We also calculate a modified energy efficiency, η_* , defined as ^3,4

$$\begin{eqnarray}&&{\eta }_{* }=\eta \displaystyle \frac{{\rm{\Delta }}{T}_{HE}+{\rm{\Delta }}{T}_{* }}{{T}_{H}-{T}_{L}}\end{eqnarray} \tag{ 10 }$$

where T_H is the hot side temperature (100 °C), ^16,17,25,26 T_L is the cold side temperature (25 °C), ΔT_HE is the heat exchanger temperature difference (5 °C), ^3,4 and ΔT_* is the temperature difference across a single process of a multi-effect and heat recovery scheme (MEHRS) (75 °C). ^3,4

Polarization curves were simulated to validate that the model cell resistance was the same as that of the experimental cell resistance and therefore matched previous experimental data. ¹¹ The conductivity values from the validated model polarization curve were then used in the time dependent discharge curve simulations. All simulations validated against experiments used the same parameters as the previously published model on this work and used the same electrolyte concentration as previous experimental studies: 0.5 M CuBr_1/2, 5 M NH₄Br, and 4 M NH₃ (negolyte only). ¹⁷

For simulations of alternative ligands and copper solubility, only the parameters directly relevant to the modification (solubility and electrode standard potential) were changed. Discharge curves were simulated with parameters and decay rate (E⁰ _pos,decay ) matched to experimental data for Selemion CMVN with an applied current density of 50 mA cm^–2. To estimate the performance difference between bromide and chloride, only the standard potential of the positive electrode was changed to E⁰ _Cl = 0.543 V from E⁰ _Br = 0.727 V, as this was the previously measured standard potential for chloride at high concentration. ⁹ We assumed that all other thermodynamic values would be the same, although it's worth noting that variables like solubility and membrane resistance, which are not considered, might exhibit differences. For results demonstrating a larger potential window, it is assumed that thermal energy input for regeneration (E_heat ) is the same as ammonia (107 Wh L^–1). If an alternative to ammonia were used to achieve the larger potential difference, the thermal energy requirement would change, but the purpose of the demonstration is to show the improvement of electrical energy output of the battery as a result of the electrolyte modification.

Validating the time-dependent model for the Cu_aq-TRAB discharge process

A comparison of simulated discharge curves with and without a crossover term confirmed the importance of ammonia crossover in replicating experimental data. Discharge curves that did not account for ammonia crossover overestimated energy and average power output by ∼20% (Fig. 1). This difference showed the limit to which decreasing ammonia transport can improve the energy and power output of the Cu_aq-TRAB. If Selemion CMVN or an alternative membrane could be modified to have significantly slower ammonia flux with the same ohmic resistance, the maximum performance would increase by a maximum of only 20%. Unavoidable polarization losses constitute the rest of the battery performance decrease (approximately 25% compared to the theoretical energy storage density calculated previously ¹⁰ ), and therefore either further chemical (i.e., electrolyte) or material (e.g., electrode, membrane) modifications would be needed to increase the energy and power output of the battery. As the decay rate increased, the discharge time and therefore energy output of the battery both decreased, as expected (Fig. 1).

Zoom In Zoom Out Reset image size

The model was further validated by fitting simulated discharge curves to previous experimental results for different membranes and different applied current densities. The fitted decay rates for each membrane were strongly correlated to previously reported ammonia flux values specific to each membrane (Fig. 2a). The strong correlation was despite large differences in resistance between the different membranes tested, indicating that the ohmic resistance and ammonia diffusion are largely independent of one another. For each membrane modeled, both the energy output and average power density were less than 5% different than the experimental discharge curves (Table S2). In these experiments and simulations, a low current density (10 mA cm^–2) was applied to the cell, thus supporting the claim that ammonia crossover is the primary source of energy loss in the system.

Zoom In Zoom Out Reset image size

Simulations at higher applied current density (50 and 100 mA cm^–2) with Selemion CMVN showed that larger decay rates (values of E⁰ _pos,decay ) were needed as the discharge current density increased (Fig. 2b, Table S3). Larger mass transfer and ohmic losses that occur at higher applied current density caused lower initial cell voltages, but they do not explain the order of magnitude increase in the fitted positive electrode standard potential decay rate. The increase in decay rate with applied current density implies that parasitic crossover increased with current density, but this increase is unlikely to be linked to an increase in ammonia (NH₃) crossover because it is uncharged and not as susceptible to electric field changes as the ionic species. Given that Selemion CMVN is a cation exchange membrane, changes in the flux of positively charged species (NH₄ ⁺ and Cu(NH₃)_x ^1/2+) are most likely to increase with higher applied current. While ammonium (NH₄ ⁺) had a much higher flux than the copper-ammonia complex, ammonium is an unlikely source of energy loss since it would primarily result in only small changes in pH. However, increase in the flux of the copper-ammonia complex would have a similar parasitic effect of ammonia (NH₃) crossover by creating a mixed potential at the positive electrode while simultaneously decreasing the total copper in the negolyte. Additionally, since the copper-ammonia complex is positively charged, the electric migration to the cathode would increase at higher applied current densities. Therefore, the source of faster decay at high polarization is likely increased Cu(NH₃)_x ^1/2+ flux.

The method of using a single decay term to modify the Nernst equation has been shown to be able to accurately simulate discharge curves for the Cu_aq-TRAB and can be extended to other flow battery chemistries as well. The methodology should be applicable to any system where the major contributor to performance decay over time is changes to the electrode thermodynamics. If there are multiple sources of parasitic loss occurring at similar rates, such as changes in surface chemistry or material degradation, then this method is unlikely to model battery discharge as accurately as is shown in this work. However, the method could be modified to non-thermodynamic performance decay in a similar manner such as decreasing exchange current density or membrane conductivity over time.

Examination of methods to improve performance

Given that the model was able to accurately simulate a wide range of materials and operating current densities, we extended the model to quantify the performance impact of different possible electrolyte modifications. First, in a previous technoeconomic analysis, it was demonstrated that the chloride ligand is much cheaper than bromide due to higher abundance and lower molecular weight, ²⁷ so we simulated polarization and discharge curves to estimate the difference in electrochemical performance between these two ligands. Next, since energy efficiency is a key performance metric for TRBs and battery theoretical energy density is linearly related to active species concentration and cell voltage, ^10,28,29 we simulated how improvement of these two variables would impact system performance. Last, we discuss potential methods and considerations for increasing copper solubility and cell voltage to help provide guidance for potential future research on Cu_aq-TRBs.

Quantifying pathways for improved Cu_aq-TRAB performance

Using chloride instead of bromide resulted in a 32% decrease in performance for both polarization and discharge curves (Fig. 3). In the polarization curve, the open circuit potential for chloride was 190 mV lower (1.04 vs. 0.85 V) than bromide due to differences in thermodynamic potential, which resulted in a peak power density of 57 mW cm⁻² instead of 88 mW cm⁻² (Fig. 3a). For the discharge curve, the average power density when using chloride was 18 compared to 27 mW cm^–2 with bromide, and the standard energy efficiency was 3.9% compared to 5.8% (Fig. 3b). The previous economic analysis assumed that the performance losses would be 15%, which our analysis indicates is an underestimation. However, it was reported that the levelized cost of storage of chloride was 80% lower than bromide; ²⁷ therefore, even with a performance decrease of 30% that is predicted here, the chloride ligand would still likelybe economically more favorable.

Zoom In Zoom Out Reset image size

Larger copper concentrations gave diminishing returns on energy output because as the discharge time increased, so did the time for decay of the positive electrode potential (Fig. 4). To achieve twice the energy output compared to the 0.5 M concentration, it was necessary to raise the copper concentration to 1.25 M. However, this change led to an 18% reduction in the average power density, reducing it from 27 to 22 mW cm^–2. The use of the 0.5 M copper concentration would reduce crossover issues and would result in a 44% decrease in energy output compared to the 1.25 M solution, but a 10% increase in average power density. Previously it has been shown that larger copper concentrations did not impact the power output of the cell, thus performance changes should be limited to energy output. ¹⁰ Lower copper concentration could be used in tandem with lower ammonia concentration, which would result in a slower decay rate due to the ammonia concentration gradient being smaller. Lower ammonia concentration has been previously shown to decrease the thermal energy requirement for regeneration (a 20% decrease when moving from 4 M to 1 M NH₃ ¹⁷ ), which would help increase the energy efficiency at low copper concentrations.

Zoom In Zoom Out Reset image size

Larger standard potential differences resulted in increases across all performance metrics with an increase of 0.2% per mV increase. At +150 mV, the peak power density increased to 120 mW cm^–2, the average power to 35 mW cm^–2, and the energy efficiency to 7.6%/5.8% (η/η_*) (Fig. 5). There was little difference in the discharge time despite the larger standard potential difference. The current density at which peak power density occurred increased to 220 mA cm^–2 from 190 mA cm^–2, therefore increasing the possible range of applied current density for discharging the battery. Additionally, the "optimal" current density for discharging would also increase by a similar percentage from 50 to 58 mA cm^–2.

Zoom In Zoom Out Reset image size

Methods and considerations for increasing copper solubility and cell voltage

Copper solubility can be increased by using a higher ligand concentration, higher temperature, or electrolyte additives such as ionic liquids and deep eutectic solvents. Multiple previous studies have shown that increasing ligand concentration enables higher copper solubility because of the higher degree of complexation of the copper-ligand complex. ^10,30–32 Provided that the dissolution reaction of many copper-ligand complexes are endothermic, ³³ increasing temperature would also increase the solubility limit. This is in line with past studies on TRABs that have shown that higher temperatures are beneficial for system performance. ^17,34,35 There has been extensive research on the use of ionic liquid and deep eutectic solvents for increasing the solubility of metal ions in conventional redox flow batteries as a means to increase their energy density. ³⁶ While many studies have focused on vanadium and zinc chemistries, the limited studies that have investigated copper metal have shown that solubilities up to 3 M are possible, ^37–40 but research on solubility of Cu(I) complexes is scarce and has shown that solubility limits tend to be lower. ^10,31,41

The standard potential difference between the electrolytes can be increased by further increasing the ligand concentration or by finding alternative ligands to bromide and ammonia. Using higher ligand concentrations has been previously shown to push the reaction potential further away from the standard electrochemical potential, but there are diminishing returns, as the degree of complexation of the copper ions is limited. ⁹ When considering alternative ligands, it is critical that they stabilize both aqueous oxidation states of copper to match the high coulombic efficiency of the Cu_aq-TRAB. Copper (II) is generally stable with most ligands, but copper (I) is more sensitive, and some ligands cause it to disproportionate, which would not be compatible with the Cu_aq-TRB. For an alternative to bromide, the only weaker field ligand in the spectrochemical series is iodide, but the copper iodide and ammonium iodide salts have low solubility in water. There are many ligands outside the spectrochemical series that could increase the standard potential higher than that of bromide but are relatively uninvestigated in copper-based electrolytes. ⁴²

For alternative negative ligands, there are multiple stronger field ligands in the spectrochemical series, but not all of them are thermally separable below the boiling point of water. It is imperative that the separation of the negative ligand happens before water vaporization because if the water vaporizes first, then the copper and ammonium salts will solidify within the thermal regeneration unit. Other considerations are less imperative for operation, but vital for maintaining and improving the energy efficiency of the Cu_aq-TRB system. Thermal separation between ligand and water should occur in one step. Multi-step separation has previously been shown to have larger thermal energy requirements than single step, ⁴³ resulting in lower system efficiency, higher system footprint, and higher capital costs. While it is possible that these drawbacks will not be true for all ligands, it is highly likely and therefore strongly recommended to avoid azeotropic systems. The ligand and metal-ligand complexes should have high solubility in water. The copper solubility limit of the Cu_aq-TRAB is ∼0.6 M, ¹⁰ which is well below traditional redox flow battery technologies and significantly limits the theoretical energy density and therefore energy efficiency of the system. Metal-ligand solubility limits should be similar or greater than 0.6 M to not drastically decrease system performance. The ligand should have low production cost. Due to its ubiquity in agricultural use as a fertilizer, ammonia is produced in extremely large quantities and is therefore a cheap reagent. While alternatives do not need to be (and are unlikely to be) as cheap as ammonia, an expensive ligand would have to have far superior performance to keep system capital costs down. The molecular size of the alternative ligand should be larger than ammonia. Due to the small molecular size and neutral charge of ammonia, control of parasitic crossover from negative to positive electrolyte is difficult, resulting in energy density losses. ¹¹ A larger ligand would diffuse slower through the membrane due to steric hindrance and therefore parasitic losses would decrease. The new electrolytes can be modified to have similar pHs. Currently, the background electrolyte salt of the Cu_aq-TRAB is ammonium bromide, a weak acid, and the negative ligand is ammonia, a strong base because it comes in the form of ammonium hydroxide. This creates a pH differential between the positive and negative electrolytes of roughly 5–6 pH units. As the battery is discharged and the electrolytes mix, the pH of each electrolyte trends to a more neutral pH, causing issues of thermodynamic stability of the copper species. Changing the ligand such that the positive and negative electrolytes can have similar pHs would avoid this phenomena and lead to less issues with long-duration battery operation.

A time-dependent model of the Cu_aq-TRAB was created to simulate discharge curves without direct simulating of species diffusion through the membrane. This model used an empirically-fitted term to represent the decay of the standard potential of the positive electrode and account for the impact of parasitic species crossover. This model was shown to accurately simulate discharge curves for a range of membranes, and the decay rate was strongly correlated to the flux of ammonia, which was previously shown to be a primary source of energy losses in TRABs. At high applied current densities, larger decay rates were required to match discharge curves, and this was attributed to increased crossover of the Cu(NH₃)^1/2+ species. The model was then used to predict the impact on performance different potential future pathways would have on Cu_aq-TRB development. Attractive strategies to further increase the performance of the Cu_aq-TRB were identified and discussed. The model results showed that changing from bromide to chloride would likely result in approximately a 30% decrease in performance. However, that decrease could be counteracted by increasing copper solubility or finding new methods to increase the standard potential difference, such as the use of an alternative negative ligand to ammonia.

The writing and editorial work of the first author were performed in part under the auspices of the U.S. Department of Energy by Lawrence Livermore National Laboratory under Contract DE-AC52–07NA27344. LLNL release number LLNL-JRNL-859138-DRAFT. This material is based upon work supported by the Department of Energy under Award Number DE-FE0032030.