Zinc Electrode Shape-Change in Secondary Air Batteries: A 2D Modeling Approach

TL;DR

This paper introduces a 2D multi-physics model of zinc-air batteries that simulates zinc shape-change during cycling, validated with tomography, and used to optimize battery performance and lifespan.

Contribution

It presents the first multi-dimensional simulation of zinc-air batteries, including zinc shape-change and validation with experimental tomography data.

Findings

Zinc shape-change impacts cycle-life significantly.

Charging current is the main factor limiting cycle-life.

Optimized initial anode structure improves battery performance.

Abstract

Zinc-air batteries offer large specific energy densities, while relying on abundant and non-toxic materials. In this paper, we present the first multi-dimensional simulations of zinc-air batteries. We refine our existing theory-based model of secondary zinc-air systems. The model comprises thermodynamically consistent multi-species transport in alkaline electrolytes, formation and dissolution of metallic zinc and passivating zinc oxide, as well as multi-phase coexistence in gas diffusion electrodes. For the first time, we simulate zinc shape-change during battery cycling by modeling convection of zinc solids. We validate our model with in-situ tomography of commercial button cells. Two-dimensional volume-averaged simulations of cell voltage and zinc electrode morphology during discharge agree with these measurements. Thus, we can study how electrolyte carbonation limits shelf-life and how zinc shape-change limits cycle-life. The charging current is found to be the major contributor to cycle-life limitations. Finally, we optimize initial anode structure and charge-discharge protocols for improved performance and cycle-ability.

Tables

Table 1. Table 1: Threshold values for reconstructed data (8 Bit) segmentation. To match the correct total mass of zinc and zinc oxide, we adapt the threshold values in each measurement.

Phase	$100\%$ SOC	$47\%$ SOC	$0\%$ SOC
gas	65,5	65,5	65,5
$\mathrm{ZnO}\left({\epsilon }_{\mathrm{ZnO}}^0=0\right)$	121,5	102,5	104,5
$\mathrm{ZnO}\left({\epsilon }_{\mathrm{ZnO}}^0\ne 0\right)$	110,5	99,5	101,5
$\mathrm{Zn}$	121,5	133,5	141,5
passive material	200,5	200,5	200,5

Table 2. Table 2: General setup of the charge-discharge protocols with a galvanostatic operation mode. For the protocol B every 5-th cycle is a recovery cycle with 𝒋 charge recovery subscript superscript 𝒋 recovery charge \boldsymbol{j}^{\mathrm{recovery}}_{\mathrm{charge}} and 𝒋 discharge recovery subscript superscript 𝒋 recovery discharge {\boldsymbol{j}^{\mathrm{recovery}}_{\mathrm{discharge}}} . The number of cycles in protocol E is only limited by a hard limit on the computation time, indicated by a “ > > ”.

protocol number	A	B	C	D	E
${𝒋}_{\mathrm{charge}}$ / $A\cdot m^{-2}$	100	100	100	100	25
${𝒋}_{\mathrm{discharge}}$ / $A\cdot m^{-2}$	100	100	100	25	100
min SOC / %	21	21	40	21	21
max SOC / %	90	90	60	90	90
recovery	-	5	-	-	-
${𝒋}_{\mathrm{charge}}^{\mathrm{recovery}}$ / $A\cdot m^{-2}$	-	10	-	-	-
${𝒋}_{\mathrm{discharge}}^{\mathrm{recovery}}$ / $A\cdot m^{-2}$	-	10	-	-	-
achieved cycles	19	23	76	32	$>$233
utilizable capacity / $A\cdot h$	3.9	3.9	4.5	6.5	$>$48

Equations

$${\partial }_t\Psi =-\nabla \left(\mathbf{N}+\Psi \dot{\mathbf{x}}\right),$$

$$⟨{\Psi }_{\alpha }⟩$$

$${\partial }_t⟨{\Psi }_{\alpha }⟩=-\nabla \cdot ⟨{\mathbf{N}}_{{\Psi }_{\alpha }}⟩-\nabla \cdot ⟨{\Psi }_{\alpha }{\dot{\mathbf{x}}}_{\alpha }⟩+{\mathring{\mathcal{P}}}_{{\Psi }_{\alpha }}^{\Theta },$$

$${\overline{\Psi }}_{\alpha }$$

$${\mathring{\Psi }}_{\alpha }$$

$$⟨{\Psi }_{\alpha }⟩=\frac{\left|V_{\alpha }\right|}{\left|V\right|}{\overline{\Psi }}_{\alpha }={\epsilon }_{\alpha }{\overline{\Psi }}_{\alpha }.$$

$$⟨1⟩=\sum _{\alpha }⟨1_{\alpha }⟩=\sum _{\alpha }{\epsilon }_{\alpha }={\epsilon }_{\mathrm{l}}+{\epsilon }_{\mathrm{Zn}}+{\epsilon }_{\mathrm{ZnO}}+{\epsilon }_{\mathrm{P}}+{\epsilon }_{\mathrm{g}}.$$

$${\partial }_t⟨1_{\mathrm{Zn}}⟩={\partial }_t{\epsilon }_{\mathrm{Zn}}$$

$${\partial }_t⟨1_{\mathrm{ZnO}}⟩={\partial }_t{\epsilon }_{\mathrm{ZnO}}$$

$${\partial }_t⟨1_{\mathrm{P}}⟩={\partial }_t{\epsilon }_{\mathrm{P}}$$

$${\overline{\dot{\mathbf{x}}}}_{\mathrm{s}}=\omega \nabla {\overline{p}}_{\mathrm{s}},$$

$$J_{\mathrm{s}}(s_{\mathrm{s}})={\overline{p}}_{\mathrm{s}}-{\overline{p}}_{\mathrm{g}}.$$

$${\partial }_t⟨{\mathcal{N}}_{\mathrm{pp}}⟩=-\nabla \cdot \left(⟨{\mathcal{N}}_{\mathrm{pp}}⟩{\overline{\dot{\mathbf{x}}}}_{\mathrm{s}}\right).$$

$${\epsilon }_{\mathrm{Zn}}$$

$${\epsilon }_{\mathrm{ZnO}}$$

$$0=\nabla \cdot ⟨{\mathbf{j}}_{\mathrm{s}}⟩+F{\mathring{\mathcal{P}}}_{\mathrm{q}}^{\Theta ,\mathrm{s}},$$

$$⟨{\mathbf{j}}_{\mathrm{s}}⟩\approx -{\overline{\kappa }}_{\mathrm{s}}{\epsilon }_{\mathrm{s}}{\mathbf{T}}_{\mathrm{s}}^{\mathrm{eff}}\cdot \nabla \cdot {\overline{\varphi }}_{\mathrm{s}}.$$

$${\mathbf{T}}_{\alpha }^{\mathrm{eff}}$$

$${\partial }_t⟨c_{{\gamma }_{\mathrm{l}}^{-}}⟩={\partial }_t\left({\epsilon }_{\mathrm{l}}{\overline{c}}_{{\gamma }_{\mathrm{l}}^{-}}\right)$$

$${\partial }_t⟨c_{{\gamma }_{\mathrm{l}}^0}⟩={\partial }_t\left({\epsilon }_{\mathrm{l}}{\overline{c}}_{{\gamma }_{\mathrm{l}}^0}\right)$$

$$0={\partial }_t⟨{\rho }_{\mathrm{l}}^{\mathrm{q}}⟩$$

$$⟨{\mathbf{N}}_{{\gamma }_{\mathrm{l}}^{-}}⟩$$

$$⟨{\mathbf{N}}_{{\gamma }_{\mathrm{l}}^0}⟩$$

$$⟨{\mathbf{j}}_{\mathrm{l}}⟩$$

$$\sum _{{\beta }_{\mathrm{l}}}{M}_{{\beta }_{\mathrm{l}}}^{\text{mol}}⟨{\mathbf{N}}_{{\beta }_{\mathrm{l}}}⟩\equiv 0,\sum _{{\beta }_{\mathrm{l}}}{z}_{{\beta }_{\mathrm{l}}}F⟨{\mathbf{N}}_{{\beta }_{\mathrm{l}}}⟩=⟨{\mathbf{j}}_{\mathrm{l}}⟩$$

$$⟨1_{\mathrm{l}}⟩=\sum _{{\beta }_{\mathrm{l}}}⟨{\nu }_{{\beta }_{\mathrm{l}}}{c}_{{\beta }_{\mathrm{l}}}⟩\approx \sum _{{\beta }_{\mathrm{l}}}{\overline{\nu }}_{{\beta }_{\mathrm{l}}}⟨c_{{\beta }_{\mathrm{l}}}⟩$$

$${\partial }_t⟨1_{\mathrm{l}}⟩={\partial }_t{\epsilon }_{\mathrm{l}}=-\nabla \cdot ⟨{\dot{\mathbf{x}}}_{\mathrm{l}}⟩-\sum _{{\beta }_{\mathrm{l}}}{\overline{\nu }}_{{\beta }_{\mathrm{l}}}\nabla \cdot ⟨{\mathbf{N}}_{{\beta }_{\mathrm{l}}}⟩+\sum _{{\beta }_{\mathrm{l}}}{\overline{\nu }}_{{\beta }_{\mathrm{l}}}{\mathring{\mathcal{P}}}_{{\beta }_{\mathrm{l}}}^{\Theta }.$$

$$⟨c_{{\gamma }_{\mathrm{l}}}{\dot{\mathbf{x}}}_{\mathrm{l}}⟩\approx {\overline{c}}_{{\gamma }_{\mathrm{l}}}{\epsilon }_{\mathrm{l}}{\overline{\dot{\mathbf{x}}}}_{\mathrm{l}}=-{\epsilon }_{\mathrm{l}}{\overline{c}}_{{\gamma }_{\mathrm{l}}}\frac{s_{\mathrm{l}}{B}_{\mathrm{l}}^{\mathrm{in}}}{{\overline{\eta }}_{\mathrm{l}}}\nabla {\overline{p}}_{\mathrm{l}},$$

$$J_{\mathrm{l}}(s_{\mathrm{l}})={\overline{p}}_{\mathrm{c}}={\overline{p}}_{\mathrm{g}}-{\overline{p}}_{\mathrm{l}},$$

$${\mathring{\mathcal{P}}}_{{\beta }_{\alpha }}^{\Theta }=\sum _{{\alpha }^{\prime }}\frac{\left|A_{\alpha ,{\alpha }^{\prime }}\right|}{|V|}\sum _{i\in \mathcal{I}}{\chi }_i^{{\beta }_{\alpha }}{j}_i^{\alpha ,{\alpha }^{\prime }},$$

$${\mathring{\mathcal{P}}}_{\mathrm{q}}^{\Theta ,\mathrm{l}}=-{\mathring{\mathcal{P}}}_{\mathrm{q}}^{\Theta ,\mathrm{s}}=\sum _{{\alpha }^{\prime }}\frac{\left|A_{\alpha ,{\alpha }^{\prime }}\right|}{|V|}\sum _{i\in \mathcal{I}}{\chi }_i^{e^{-}}{j}_i,$$

$$\left(4-2{\epsilon }_{\mathrm{f}}\frac{V_{\mathrm{Zn}}^{\text{mol}}}{V_{\mathrm{ZnO}}^{\text{mol}}}\right)j_{\mathrm{I}}={\left(1-{\epsilon }_{\mathrm{f}}\right)}^{3.5}{\overline{D}}_{{\mathrm{OH}}^{-}}\frac{{\overline{c}}_{{\mathrm{OH}}^{-}}-c_{\mathrm{s},{\mathrm{OH}}^{-}}}{{\delta }_{\mathrm{ZnO}}}\frac{r_{\mathrm{ZnO}}}{r_{\mathrm{Zn}}},$$

Full text

1 Introduction

Electrochemical energy storage plays a key role for the transition to a carbon neutral society [1]. Energy-efficient lithium-ion batteries are optimally suited for portable electronics and local electric mobility. Scalable electrolysers and fuel cells address the needs of long-term grid-scale energy storage. In between, metal-air batteries promise high specific energy densities at low costs for household-level energy storage or range extenders [2, 3, 4]. Possible material choices are, with increasing energy density, iron, zinc, magnesium, aluminum, and lithium. A couple of years ago, research concentrated on lithium-air batteries because of its high theoretical energy densities [5, 6, 7]. Unfortunately, fundamental challenges still prevent their commercialization. The reasons include low electrolyte stability under oxygen evolution conditions, slow growth and dissolution of lithium peroxide, low oxygen solubility in the electrolytes, and shape-change of lithium metal.

Zinc-air batteries are a promising alternative. Zinc is cheap, ecologically friendly, nonflammable, globally abundant, and non-toxic [2, 3, 8]. Most importantly, zinc metal is stable in aqueous electrolytes aside from controllable hydrogen evolution [9]. For aqueous electrolytes, gas diffusion cathodes offer fast oxygen transport [7, 10]. Therefore, primary zinc-air cells are commercially available for decades. However, their cycle-life and shelf-life is limited. Firstly, the standard alkaline electrolytes absorb carbon dioxide from air, which reduces their pH value and limits shelf-life to a few months [11]. Possible mitigation methods include carbon dioxide absorbers [12] and novel electrolytes [13], e.g., near-neutral ammonium chloride solutions [14, 15, 16, 17] or ionic liquids [18, 19, 20]. Secondly, porous zinc anodes alter their shape during battery cycling, which reduces porosity and limits cycle-life [8]. Novel zinc structures are proposed to reduce shape-change [21], e.g., hyper-dendritic zinc foam [22], 3D zinc sponge [23], or mechanical coatings [24, 25].

In this article, we present the first multi-dimensional simulations of zinc-air batteries. Our 2D simulations of button cells are enabled by local volume-averaging theory [26, 27]. We optimize cell design and complement experimental efforts. Our electrochemical model elucidates the intimate correlation between electrode morphology and electrolyte dynamics during battery cycling. We parametrize and validate our model by x-ray tomography data [28, 29]. Recently, we developed a theory-based model of aqueous alkaline zinc-air batteries that consistently couples the dynamics of its solid, liquid, and gaseous phases [11]. 1D simulations predicted the effects of inhomogeneous zinc oxide nucleation and growth as well as carbonate formation in the electrolyte. We validated the model predictions of discharge voltage and shelf-life. Here, we extend our model by convection of solid zinc metal and simulate zinc shape-change in secondary zinc-air batteries.

Preceding zinc-air battery models [29, 30, 31, 32, 33, 34] did not achieve such a level of consistency and predictive powers as discussed in Refs. [4, 11]. These 0D or 1D simulations neglect electrochemical interactions and local material inhomogeneities. However, 3D simulations advance in the field of lithium-ion batteries [35, 36, 37, 38, 39, 40, 41]. A few multi-dimensional models are applied to conversion type batteries, mainly analyzing single electrodes. For lithium-air batteries, our group performs 3D Lattice-Boltzmann simulations of phase coexistence in gas diffusion electrodes [10]. Dendrite formation in lithium metal anodes is measured and simulated in Ref. [42]. We are aware of a single 3D simulation of a lithium-air battery cell, which has a limited resolution and shows only slight 3D effects [43].

We proceed with descriptions of our x-ray experiments, our theory-based model, and our computational implementation in Sec. 2, Sec. 3, and Sec. 4, respectively. We simulate discharge at various currents for comparison with tomographies [28] and 1D simulations [11] (see Sec. 5). Finally, we discuss strategies to reduce zinc shape-change in Sec. 6.

2 Experiment

We simulate the Renata ZA13 button cell [44] (see Fig. 1a). Air flows through the lids into the gas diffusion electrode (GDE), where it dissolves in the electrolyte and the oxygen reaction is catalyzed. Most of the cell volume is occupied by the negative electrode. In the fully charged state, it contains a porous zinc powder. The void space in the top region below the negative pole compensates the volume expansion from zinc to zinc oxide during discharge.

For parametrization and validation, we discharge the Renata battery [44] using a Gamry 3000 potentiostat at an ambient temperature of $T=20$*∘*C. The galvanostatic discharge current per cross section area is 54.2\text{\,}\mathrm{A}\text{{\cdot}}{\mathrm{}}^{-2}\text{{\cdot}}\mathrm{m}. The cross section of the cell on the cathode side is $\sim$4.2\text{$\cdot$}{10}^{-5}\text{,}{\mathrm{m}}^{2}$$. A conventional X-ray tube is used for tomographic measurements in combination with a Hamamatsu X-ray flat panel detector with $2316\times 2316$ pixel, resulting in a pixel size of $9\text{\,}\mathrm{\SIUnitSymbolMicro m}$. Details of the setup are given in LABEL:fig:suppl:misc:setup and Refs. [28, 29, 45]. An accelerating voltage of $120\text{\,}\mathrm{kV}$ is applied while the current at the tungsten anode is adjusted to $84\text{\,}\mathrm{\SIUnitSymbolMicro A}$. A $0.5\text{\,}\mathrm{mm}$ Cu filter is used in order to avoid beam hardening. We take 1200 projections for a full three-dimensional tomography. Each angle step is exposed for three times for $1\text{\,}\mathrm{s}$ in order to increase the signal-to-noise-ratio.

We perform tomographic imaging at various state-of-charges during discharge until the cut-off potential $0.6\text{\,}\mathrm{V}$ is reached (see Fig. 1b). Before each tomography, the applied current is switched off and the voltage is allowed to relax. Each tomography under steady-state condition takes about $1.5\text{\,}\mathrm{h}$. We describe data processing in Section 4.3 below.

3 Theory and Model

In this section, we present a thermodynamically consistent and volume-averaged multi-dimensional model for rechargeable zinc-air cells. It is based on our 1D model developed in Refs. [11, 7]. We model the phase coexistence of solid ($\alpha=\mathrm{s}$), liquid ($\alpha=\mathrm{l}$), and gaseous phases ($\alpha=\mathrm{g}$). For each state of matter $\alpha$, we consider distinct species $\beta_{\alpha}$.

As solids, we take into account zinc ($\beta_{\mathrm{s}}=\mathrm{Zn}$), zinc oxide ($\beta_{\mathrm{s}}=\mathrm{ZnO}$), and passive material ($\beta_{\mathrm{s}}=\mathrm{P}$), e.g., housing, separator. Our recent 1D model discussed nucleation, growth, and dissolution of zinc and zinc oxide [11]. We simulated how zinc oxide passivates zinc particles. In this article, we present the first model of zinc convection and simulate shape-change.

The electrolyte is a concentrated aqueous potassium hydroxide solution. Its constituents are water ($\beta_{\mathrm{l}}=\mathrm{H}_{2}\mathrm{O}$), potassium ($\beta_{\mathrm{l}}=\mathrm{K}^{+}$), hydroxide ($\beta_{\mathrm{l}}=\mathrm{OH}^{-}$), zincate ($\beta_{\mathrm{l}}=\mathrm{Zn(OH)}^{2-}_{4}$), carbonate ($\beta_{\mathrm{l}}=\mathrm{CO}^{2-}_{3}$), and dissolved oxygen ($\beta_{\mathrm{l}}=\mathrm{O}^{\mathrm{l}}_{2}$). We denote the neutral, positive, and negative species with the indices $\beta^{0}_{\mathrm{l}}$, $\beta^{+}_{\mathrm{l}}$, $\beta^{-}_{\mathrm{l}}$, respectively. Our theory-based modeling reveals that local charge and volume constraints reduce the number of independent species [11]. We choose water and potassium as dependent species and reserve the index $\gamma_{\alpha}$ for the independent species. A multi-component volume-constraint induces electrolyte convection.

From the gaseous phase, we model the inflow of oxygen ($\beta_{\mathrm{g}}=\mathrm{O}^{\mathrm{g}}_{2}$) and carbon dioxide ($\beta_{\mathrm{g}}={\mathrm{CO}_{2}}$). The gaseous phase makes possible fast oxygen transport through the GDE [7]. The electrolyte degrades because of carbon dioxide absorption. A gas reservoir on top of the zinc anode buffers the volume expansion from zinc to zinc oxide.

In the following, we present our electrochemical transport and reaction model based on Ref. [11]. In Section 3.1, we discuss the volume-averaging process that makes possible our multi-dimensional simulations. We derive equations for the transport of the solid phases ($\mathrm{Zn}$, $\mathrm{ZnO}$, $\mathrm{P}$) in Section 3.2, of the liquid phase in Section 3.3, and of the gas phase in Section 3.4. Finally, we model the relevant reactions in Section 3.5.

3.1 Volume-Averaged Model

Because on the micro-scale of a zinc-air cell multiple distinct phases are entangled, simulating their dynamics is a challenging task. We develop efficient multi-dimensional simulations based on the approach of volume-averaging as discussed by Whitaker and Chen [26, 27, 46]. Thus, we treat the composite systems as a bulk of mixed materials. In this article, we present a simplified derivation of local volume-averaging and refer to an upcoming publication for further details.

Our transport theory relies on microscopic single-phase dynamics. In the Euler frame, i.e., lab frame, the continuity equation holds

[TABLE]

where $\Psi(t,x)$ represents the conserved variables, e.g., concentration. $\dot{\boldsymbol{x}}$ is the center-of-mass convection velocity and $\boldsymbol{N}$ is the non-convective flux density. Here, we omit source terms $\mathcal{P}^{\Omega}=0$, as we get by without homogeneous bulk reactions.

We choose an averaging domain $V$ which exceeds the typical size of the micro-structure. Firstly, we define the superficial volume-average of a field variable $\Psi_{\alpha}$ in phase $\alpha$ as average over the composite domain

[TABLE]

where the integral is restricted to single phases $V_{\alpha}$ within the domain $V$. The single-phase continuity equation Eq. 1, valid in phase $\alpha$, transforms into the volume-averaged form

[TABLE]

with the surface reaction rates $\mathring{\mathcal{P}}^{\Theta}_{\alpha}$ specified in Section 3.5. In the next three sections, Section 3.2, Section 3.3, and Section 3.4, we express the superficial volume-averages of the transport terms with averages over single phases.

Secondly, we define the intrinsic volume-average $\overline{\Psi}_{\alpha}$ and the local surface-average $\mathring{\Psi}_{\alpha}$ as average over a single phase

[TABLE]

with the surface $A_{\alpha}$ of phase $\alpha$ within the domain $V$. The two volume-averages are connected by the volume fractions $\varepsilon_{\alpha}$

[TABLE]

The volume fractions sum to unity

[TABLE]

In the following, we discuss the transport equations for solid, liquid, and gaseous phases.

3.2 Solid Phase Transport

In metal-air batteries, solid phases grow and dissolve during cycling. On the mesoscopic particle-scale, phase-field-like models can capture the dynamics of phase boundaries [6, 48, 49, 50, 51]. On the macroscopic cell-scale, more computationally efficient approaches would be necessary. We note that modern cell models follow the dynamics of volume fractions locally and neglect the non-local movement of moving phases [52]. In this case, the breathing of solid volume is just compensated by a local porosity change. However, operando tomographies like Fig. 1 show an additional convection in the $\mathrm{Zn}$ anode because of different molar volumes $V^{\text{mol}}_{\mathrm{ZnO}}=1.58~{}V^{\text{mol}}_{\mathrm{Zn}}$ [28]. We are the first to capture this phenomenon and model convection of solids in a cell model.

We distinguish three phases that make up the total solid volume, i.e., $\varepsilon_{\mathrm{s}}=\varepsilon_{\mathrm{Zn}}+\varepsilon_{\mathrm{ZnO}}+\varepsilon_{\mathrm{P}}$. The volume-averaged continuity equations for the solid phases are (see Eq. 3)

[TABLE]

where the solid convection velocity $\overline{\dot{\boldsymbol{x}}}_{\mathrm{s}}$ applies to both phases consistent with our primitive particle model (see Fig. 2). In analogy to laminar Darcy flow of liquids through porous structures [7], it is calculated from a pseudo-pressure $\overline{p}_{\mathrm{s}}$

[TABLE]

with the pseudo-permeability $\omega$.

Dependent on the solid saturation $s_{\mathrm{s}}={\left(\varepsilon_{\mathrm{Zn}}+\varepsilon_{\mathrm{ZnO}}\right)}/{(1-\varepsilon_{\mathrm{P}})}$, we determine the solid pseudo-pressure $\overline{p}_{\mathrm{s}}$ with a phenomenological Leverett function $J_{\mathrm{s}}(s_{\mathrm{s}})$,

[TABLE]

We choose a Leverett function $J_{\mathrm{s}}(s_{\mathrm{s}})$ with a step at the critical saturation $s^{*}_{\mathrm{s}}$. If the saturation is small, $s_{\mathrm{s}}<s^{*}_{\mathrm{s}}$, the Leverett function vanishes, $J_{\mathrm{s}}(s_{\mathrm{s}})=0$. In this first case, the solids grow locally and the solid convection vanishes, i.e., $\overline{\dot{\boldsymbol{x}}}_{\mathrm{s}}=0$. However, if the saturation exceeds the critical saturation, $s_{\mathrm{s}}>s^{*}_{\mathrm{s}}$, the Leverett function increases with the saturation $s_{\mathrm{s}}$. In this second case, the pseudo-pressure pushes the solids away and the solid convection is non-zero, i.e., $\overline{\dot{\boldsymbol{x}}}_{\mathrm{s}}\neq 0$. A mathematical description of $J_{\mathrm{s}}(s_{\mathrm{s}})$ is given in the Supplementary Materials in LABEL:sec:suppl:params:levs. Note that alternative approaches to model $\overline{p}_{\mathrm{s}}$ are available [49, 50].

Now, we give details on our microscopic zinc particle model. X-ray/SEM measurements show that a zinc oxide shell forms around the zinc particles [47]. We assume spherical $\mathrm{Zn}$ particles surrounded by a porous $\mathrm{ZnO}$ shell as shown in Fig. 2. Our model solves for the density of primitive particles $\left\langle\mathcal{N}_{\mathrm{pp}}\right\rangle(x,t)$, too. As primitive particles cannot emerge or vanish, they move with the solid convection velocity only,

[TABLE]

According to Fig. 2, the $\mathrm{Zn}$ radius $r_{\mathrm{Zn}}$ and the $\mathrm{ZnO}$ radius $r_{\mathrm{ZnO}}$ are then calculated from

[TABLE]

Here, $\varepsilon_{\mathrm{f}}=0.95$ is the solid volume fraction of the $\mathrm{ZnO}$ shell. Further details are given in the Supplementary Materials (see LABEL:sec:suppl:params:sa).

Electrons are conducted by the solids in anode and cathode. The spatially-resolved electrode potential is determined by

[TABLE]

with the volume-averaged flux term

[TABLE]

The effective transport tensor $\boldsymbol{T}^{\mathrm{eff}}_{\alpha}$ has to be modeled (see for example Refs. [27, 53, 54]). We model a macroscopic and isotropic medium which reduces $\boldsymbol{T}^{\mathrm{eff}}_{\alpha}$ to a scalar tortuosity factor

[TABLE]

with the Bruggemann coefficient $\mathrm{b}$ and the identity matrix $\boldsymbol{I}$.

The parameters for solid phase transport are stated in the Supplementary Materials in LABEL:sec:suppl:params:solid [11, 47, 55, 56].

3.3 Liquid Phase Transport

In this subsection, we summarize the consistent multi-species electrolyte model developed by Stamm et al. [11]. Our transport theory ensures a strictly positive entropy production and is thermodynamic consistent [35, 36]. Since we describe transport in the center-of-mass frame and consider the electro-neutrality of the electrolyte on macroscopic length scales, it suffices to calculate the motion of anions ($\mathrm{OH}^{-}$, $\mathrm{Zn(OH)}^{2-}_{4}$, $\mathrm{CO}^{2-}_{3}$) and dissolved oxygen $\mathrm{O}^{\mathrm{l}}_{2}$. We state Eq. 3 for the independent concentrations $c_{\gamma_{\mathrm{l}}}$ and the vanishing charge density $\rho^{\mathrm{q}}_{\mathrm{l}}=0$ in the porous electrode and separator,

[TABLE]

The volume-averaged center-of-mass flux densities for ions $\boldsymbol{N}_{\gamma^{-}_{\mathrm{l}}}$, neutral solutes $\boldsymbol{N}_{\gamma^{0}_{\mathrm{l}}}$, and charges $\boldsymbol{j}_{\mathrm{l}}$ are

[TABLE]

We use the electric potential $\phi_{\mathrm{l}}$, diffusion coefficients $D_{\gamma_{\mathrm{l}}}$, electrolyte conductivity $\kappa_{\mathrm{l}}$, the transference numbers $t_{\gamma^{-}_{\mathrm{l}}}$, chemical potentials $\mu_{\gamma_{\mathrm{l}}}$, and the charge numbers $z_{\gamma_{\mathrm{l}}}$. We approximate the averaged coefficients as $\overline{D(c)}\approx D(\overline{c})$. The effective transport tensor for an isotropic medium is given in Eq. 18. The flux densities of $\mathrm{H}_{2}\mathrm{O}$ and $\mathrm{K}^{+}$ in the center-of-mass frame follow from the frame condition and the electrical current

[TABLE]

with the partial molar volumes $M^{\text{mol}}_{\beta_{\mathrm{l}}}$.

We take into account the volumetric electrolyte equation of state [7]

[TABLE]

with the partial molar volumes $\overline{\nu_{\beta_{\mathrm{l}}}(c)}\approx\overline{\nu}_{\beta_{\mathrm{l}}}(\overline{c})$. For pressure-independent partial molar volumes, Eq. 26 realizes a multi-component incompressibility constraint. As shown in Ref. [11], Eq. 26 yields the temporal evolution of porosity $\varepsilon_{\mathrm{l}}$ in electrodes and separator,

[TABLE]

We note that Eq. 27 generalizes the standard incompressibility condition ${\boldsymbol{\nabla}\cdot{\dot{\boldsymbol{x}}}_{\mathrm{l}}=0}$ known for single-component and single-phase systems. [36].

For realistic discharge rates of zinc-air batteries, the flow of the alkaline electrolyte is laminar. Hence, the center-of-mass velocity ${\overline{\dot{\boldsymbol{x}}}_{\mathrm{l}}}$ is the gradient of a scalar potential. We apply Darcy’s law and determine the convective flux density $\left\langle c_{\gamma_{\mathrm{l}}}{\dot{\boldsymbol{x}}}_{\mathrm{l}}\right\rangle$ with the liquid pressure $\overline{p}_{\mathrm{l}}$ [7],

[TABLE]

where $B_{\mathrm{l}}^{\mathrm{in}}$ is the intrinsic permeability and $\overline{\eta}_{\mathrm{l}}$ is the intrinsic viscosity of the liquid phase. The liquid saturation $s_{\mathrm{l}}={\varepsilon_{\mathrm{l}}}/{(1-\varepsilon_{\mathrm{s}})}$ serves as relative permeability and takes into account the multi-phase character of the system [57]. The liquid pressure $\overline{p}_{\mathrm{l}}$ is determined by the phenomenological Leverett function $J_{\mathrm{l}}(s_{\mathrm{l}})$ [10, 58]. It holds

[TABLE]

with the capillary pressure $\overline{p}_{\mathrm{c}}$.

The corresponding transport parameters are given in the Supplementary Materials in LABEL:sec:suppl:params:electrolyte [59, 11, 60, 61, 62, 63, 64, 65, 66, 67, 68, 30, 69, 60, 70, 59, 58, 10].

3.4 Gas Phase Transport

In a zinc-air button cell, the gas phase is present in anode, separator, and cathode. The gaseous volume fraction $\varepsilon_{\mathrm{g}}$ follows from the constraint on phase volumes in Eq. 7 as $\varepsilon_{\mathrm{l}}$, $\varepsilon_{\mathrm{Zn}}$, $\varepsilon_{\mathrm{ZnO}}$, and $\varepsilon_{\mathrm{P}}$ are known from Eq. 27, Eq. 8, Eq. 9, and Eq. 10, respectively.

As shown in Ref. [7], the partial pressures in the gas phase remain almost constant during discharge of metal-air batteries with aqueous alkaline electrolytes. Therefore, we keep the partial gas pressures $\overline{p}_{\mathrm{O}^{\mathrm{g}}_{2}}$ and $\overline{p}_{\mathrm{CO}_{2}}$ constant. The total gas pressure is the atmospheric pressure $\overline{p}_{\mathrm{g}}=p_{\mathrm{atm}}^{\ast}$.

3.5 Reaction Rates

We depict the reactions in the zinc-air cell in Fig. 2 and denote them with the indices $\mathcal{I}=\left\{\mathrm{I},\mathrm{II},\mathrm{III},\mathrm{IV},\mathrm{V}\right\}$. The total surface source terms $\mathring{\mathcal{P}}_{\beta_{\alpha}}^{\Theta}$ are calculated as sums of products of the total reaction rates per surface area $j_{i}$ with their specific surface areas ${\left|A_{\alpha,\alpha^{\prime}}\right|}/{\left|V\right|}$,

[TABLE]

with the stoichiometric coefficients $\chi^{\beta_{\alpha}}_{i}$ of species $\beta_{\alpha}$ in reaction $i$. The sources for ionic charge $\mathring{\mathcal{P}}_{\mathrm{q}}^{\Theta,\mathrm{l}}$ and $\mathring{\mathcal{P}}_{\mathrm{q}}^{\Theta,\mathrm{s}}$ are

[TABLE]

with the stoichiometric coefficients $\chi^{e^{-}}_{i}$ of electrons $e^{-}$ in reaction $i$.

The individual reaction rates $j_{i}$ are stated in the Supplementary Materials (see LABEL:sec:suppl:params:reac). Note that we summarize carbon dioxide absorption and carbonate formation in a single surface reaction as shown by Stamm et al. [11]. Thus, all volume source terms vanish in our model.

We introduced our microscopic model of spherical particles in Section 3.2. The $\mathrm{ZnO}$ shell around the $\mathrm{Zn}$ core hinders the ion transport from the bulk electrolyte to the $\mathrm{Zn}$ surface. Hydroxide is the main reactant in zinc dissolution at the $\mathrm{Zn}$ surface. Thus, we model the diffusion of hydroxide through the $\mathrm{ZnO}$ shell assuming spherically symmetric $\mathrm{Zn}$ particles and $\mathrm{ZnO}$ shells (see Fig. 2). We calculate the surface concentration $c_{\mathrm{s},\mathrm{OH}^{-}}$ by a radial diffusion from the bulk to the $\mathrm{Zn}$ surface [11]

[TABLE]

with the reaction rate $j_{\mathrm{I}}$, its thickness $\delta_{\mathrm{ZnO}}=r_{\mathrm{ZnO}}-r_{\mathrm{Zn}}$, and the inner and outer radius of the passivation layer $r_{\mathrm{Zn}}$ and $r_{\mathrm{ZnO}}$. The elevated Bruggemann coefficient is motivated by low porosities of the $\mathrm{ZnO}$ shell $1-\varepsilon_{\mathrm{f}}=0.05$. The factor $\left(4-2\varepsilon_{\mathrm{f}}V^{\text{mol}}_{\mathrm{Zn}}/V^{\text{mol}}_{\mathrm{ZnO}}\right)$ is the number of hydroxide ions passing the $\mathrm{ZnO}$ shell per dissolved zincate ion. Here, we assume that the dissolved $\mathrm{Zn}$ phase is completely replaced by $\mathrm{ZnO}$. Thus, hydroxide ions are created below the passivation shell, too.

Details on reactions and transport parameters are stated in the Supplementary Materials in LABEL:sec:suppl:params:reac [71, 72, 11, 30, 68, 33, 73, 7, 74, 67, 75, 76, 77, 78, 79, 80, 81].

4 Computational Details

4.1 Parameters

Parameters represent the Renata ZA13 button cell [44]. Solid phase, liquid phase, and reaction parameters for the aqueous alkaline electrolyte (32 weight percent KOH) are discussed by Stamm et al. [11]. They are stated in the Supplementary Materials in LABEL:sec:suppl:params:solid, in LABEL:sec:suppl:params:electrolyte, and in LABEL:sec:suppl:params:reac, respectively.

Regarding the structure of the zinc electrode during cycling, the initial radius of the primitive zinc particles was measured with x-ray tomography [45], the formation of zinc oxide shells was observed with scanning electron microscopy [47], and the solid phase movement during discharge was imaged with x-ray tomography [28, 29]. We model the dissolution and growth of primitive $\mathrm{Zn}$/$\mathrm{ZnO}$ particles accordingly (see LABEL:sec:suppl:params:sa). Our phenomenological Leverett function $J_{\mathrm{s}}(s_{\mathrm{s}})$ is adjusted to yield the observed solid convection during discharge (see LABEL:sec:suppl:params:levs).

4.2 Numerics

Our in-house tool, Battery and Electrochemistry Simulation Tool (BEST), provides a framework for multi-dimensional finite-volume discretized electrochemical simulations [35, 36]. BEST was originally developed at Fraunhofer ITWM for simulation of micro-structures on the meso-scale. We employ a backward differentiation formula for the time derivative and a Crank-Nicolson step for the time integration [82]. Convection is implemented with an upwind scheme [82]. The resulting nonlinear systems of equations are solved with an iterative Newton method.

We perform 2D simulations on a radially symmetric cylindrical grid which resembles the button cell geometry. We use a resolution of $(70,3,9)\times 20$ ((anode, separator, cathode) x-dir $\times$ radial discretization units) in Section 5.1 and $(50,1,3)\times 12$ in all other sections. Our resolutions are sufficiently accurate as demonstrated in the Supplementary Materials in LABEL:sec:suppl:ad. In LABEL:sec:suppl:val1D, we compare our new 2D simulations with the existing 1D code used by Stamm et al. [11] and demonstrate accurate implementation.

The multi-component incompressibility constraint in Eq. 27 stiffens the system of equations. We improve the numerical time integration to yield stable simulations. A detailed description and error analysis of this algorithm will be presented in an upcoming publication.

4.3 Construction of Initial Data

After measurements, tomographic images are reconstructed using the Octopus software [83]. Firstly, we subtract the passive cell housing. Secondly, we apply a median filter to the reconstructed data (8 Bit gray values) reducing measurement noise. Thirdly, we perform segmentation of different phases inside the 3D volume by gray value thresholding. Finally, we add the passive cell housing again. The reconstructed data has a resolution of $491\times 816\times 816$ pixels and is segmented, i.e., each voxel is associated to one phase.

Gray value thresholding relies on different beam attenuation of different phases. Dark pixels have less beam attenuation, which correspond to gas volumes. Bright values mean high absorption, which represents $\mathrm{Zn}$ here. Electrolyte and $\mathrm{ZnO}$ have similar medium attenuation coefficients, that is why they are hardly separable. A better segmentation could be achieved by stochastic methods using Fiji tools and WEKA segmentation [84, 85]. The gray value thresholds are determined at all state-of-charges from the expected $\mathrm{Zn}$ and $\mathrm{ZnO}$ amounts. By assuming that zinc is completely transformed to zinc oxide at the end-of-discharge, we estimate the initial amount of $\mathrm{Zn}$ as 5.59\text{\cdot}{10}^{-3}\text{\,}\mathrm{mol} and of $\mathrm{ZnO}$ as 6.14\text{\cdot}{10}^{-4}\text{\,}\mathrm{mol}. From the discharge conditions in Ref. [28], we then calculate the amounts of $\mathrm{Zn}$ and $\mathrm{ZnO}$ at each state-of-charge. We perform simulations with and without $\mathrm{ZnO}$ in order to determine if $\mathrm{ZnO}$ is present in a fresh cell (see Section 5.1). The resulting threshold values are presented in Table 1.

We construct the 2D initial data for the simulations by volume-averaging the 3D reconstructed data. The volume-average is performed in cylindrical coordinates, including a radial average. We average all phases within the constructed discretization units (see Figs. 1 and 4).

5 Validation and Results

In this section, we present 2D simulation results for the discharge and charge of zinc-air button cells. Firstly, we validate our 2D simulations by comparing with measurements during cell discharge (see Section 5.1). Secondly, we evaluate the effect of electrolyte carbonation on cycle-life (see Section 5.2). Finally, we simulate cycling at various current densities (see Section 5.3).

5.1 Validation During Cell Discharge

For model validation during galvanostatic discharge, we compare our simulation results with experimentally measured cell voltages and tomographic images (see Fig. 3a and Fig. 4). The voltage curve contains information about $\mathrm{ZnO}$ nucleation and growth. X-ray tomography reveals the dynamics of $\mathrm{Zn}$ electrode morphology.

It is debated whether $\mathrm{ZnO}$ is admixed to the anode initially. In Ref. [28], simulations were initialized with a small amount of solid $\mathrm{ZnO}$; in Ref. [11], we concluded that the analyzed primary zinc-air button cell contains a pure $\mathrm{Zn}$ anode initially. Therefore, we carry out two simulations. Based on tomographic images of the fresh cell, we initialize one with $\varepsilon^{0}_{\mathrm{ZnO}}=0$ and the other one with $\varepsilon^{0}_{\mathrm{ZnO}}\neq 0$ (see Table 1).

Fig. 3 compares measured and simulated discharge voltages. After $\text{SOC}=90\%$, both simulations and the experiment agree nicely, apart from a minor deviation at the end-of-discharge. The voltage curve is relatively flat as expected for a conversion-type battery chemistry [11]. At the beginning of discharge, however, the two simulations deviate from each other. The simulation with admixed $\mathrm{ZnO}$ predicts a dip of $\sim$50\text{,}\mathrm{mV}$$ which is absent in the simulations without admixed $\mathrm{ZnO}$. The inset of Fig. 3a highlights that the experiment shows a dip in cell voltage. State-of-charge and voltage amplitude of this dip very well agree to the predictions of the simulation with $\varepsilon^{0}_{\mathrm{ZnO}}=0$.

As discussed by Stamm et al. [11], the voltage dip indicates $\mathrm{ZnO}$ nucleation. If $\mathrm{ZnO}$ is not admixed to the $\mathrm{Zn}$ paste, i.e., $\varepsilon^{0}_{\mathrm{Zn}}=0$, a critical $\mathrm{Zn(OH)}^{2-}_{4}$ supersaturation has to be reached before $\mathrm{ZnO}$ can continuously precipitate. The resulting peak in the $\mathrm{Zn(OH)}^{2-}_{4}$ concentration is reflected in the dip in cell voltage. Thus, we choose a pure $\mathrm{Zn}$ paste without admixed $\mathrm{ZnO}$ as initial condition for the simulations in the remainder of this section.

For further validation, we compare the simulated $\mathrm{Zn}$ and $\mathrm{ZnO}$ distribution within the anode with tomography data during galvanostatic discharge. In Fig. 4, we show the phase distribution for the initial, the half-discharged, and the fully-discharged state. In Fig. 4a-c, we observe in simulations and experiments that $\mathrm{Zn}$ preferentially dissolves close to the separator. Therefore, at the end-of-discharge the remaining $\mathrm{Zn}$ is located near the lid at the top of the button cell. $\mathrm{ZnO}$ forms inside the $\mathrm{Zn}$ structure as shown in Fig. 4d-f. The tomography records an inhomogeneous $\mathrm{ZnO}$ distribution at intermediate state-of-charges and an accumulation near the center of the cell. Our simulation results in a slightly flatter distribution of $\mathrm{ZnO}$ with a peak close to the separator. In the simulations, $\mathrm{ZnO}$ does not precipitate in a tiny region directly near the separator (see Fig. 4g-i). The experiment, instead, reveals a larger region, with a nonzero, but significantly reduced $\mathrm{ZnO}$ volume fraction. The overall $\mathrm{ZnO}$ distributions of experiment and simulation have a qualitatively similar shape in the fully-discharged state. On closer inspection, experiments show a slight movement of $\mathrm{Zn}$ and $\mathrm{ZnO}$ to the lid, which is absent in simulations.

Let us begin with discussing the zinc distribution. We observe that our simulation is in very good agreement with the measurement. In both, simulation and experiment, we observe a slightly faster dissolution near the separator during discharge and some $\mathrm{Zn}$ remaining near the lid at the end-of-discharge. This is caused by a hydroxide concentration gradient from the cathode to the anode. The slight movement to the top of the button cell in the experiment results from the larger molar volume of $\mathrm{ZnO}$ than of $\mathrm{Zn}$. Since the anode is a $\mathrm{Zn}$ paste without any supporting structure, this volume expansion induces a movement. The separator on the bottom acts as a boundary, while there is enough space near the lid. Thus, the overall movement of the solid phase is faced towards the lid, where the gas is compressed or extruded. Our novel model for solid convection in Eq. 11 does not reproduce this slight movement during a single discharge. It induces a solid phase movement only if the solid volume saturation rises above a certain threshold. During cycling, however, our model captures the continuous movement of solid phase as shown below in Section 5.3. We highlight that this is the first cell model capturing this effect.

With respect to the zinc oxide phase, simulations also predict the overall shape of the $\mathrm{ZnO}$ block correctly. Due to its negative charge, zincate drifts towards the lid. Therefore, $\mathrm{ZnO}$ shows a clear peak structure in experiment and simulation as highlighted in Fig. 4h-i. Nevertheless, there are deviations to the tomographies in a detailed view. As our previous 1D simulations [11], our 2D simulations predict a small region close to the separator, where no $\mathrm{ZnO}$ precipitates. The gap in this work is smaller than in our previous work because of a different initial phase distribution and a different parametrization of the nucleation model. We note that the size of this region depends on the details of the nucleation model. Our tomographic images do not show a region without any $\mathrm{ZnO}$ precipitation, but a large region with low $\mathrm{ZnO}$ volume fraction close to the separator. We note that the reconstructed data of $\mathrm{ZnO}$ is insufficient for a concluding comparison because $\mathrm{ZnO}$ and electrolyte phases cannot be clearly distinguished. Higher resolution tomography with electron microscopy would allow to improve the nucleation model and to unambiguously determine the phase distribution (compare Ref. [47]).

To sum up the comparisons of the cell voltage profile and the tomography data, we find that or model is in a good agreement with the measurements, but that modeling of $\mathrm{ZnO}$ nucleation is challenging. The fresh primary cell contains a pure $\mathrm{Zn}$ anode without admixed $\mathrm{ZnO}$.

5.2 Carbonation of Electrolyte During Cycling

We simulate the effect of $\mathrm{CO}_{2}$ in the feed gas during battery cycling. In this case, $\mathrm{CO}_{2}$ is absorbed by the electrolyte and reacts to carbonate $\mathrm{CO}^{2-}_{3}$ [11]. In Fig. 3, we compare simulations for ambient air with $\mathrm{CO}_{2}$ and synthetic air without $\mathrm{CO}_{2}$. After a first discharge to $21\text{\,}\mathrm{\char 37\relax}$ state-of-charge, we galvanostatically cycle the cell between $90\text{\,}\mathrm{\char 37\relax}$ and $21\text{\,}\mathrm{\char 37\relax}$ state-of-charge. In the case of carbonation, we observe a clear drop of performance from cycle to cycle. The charging voltage increases slightly, while the discharge voltage decreases significantly. After 10 cycles ($\sim 38$ days), the cell voltage diverges at the end of discharge and the simulation terminates. In the case of pure oxygen supply, the voltage levels remain almost constant. After 24 cycles ($\sim 92$ days) at the end of the charging, the cell voltage diverges.

In Ref. [11], we discussed in detail how carbonation is limiting the shelf-life of alkaline zinc-air batteries during storage. Here, we summarize this effect for battery cycling. Electrolyte carbonation due to reaction (V) decreases the overall $\mathrm{OH}^{-}$ concentration and reduces electrolyte conductivity. Because $\mathrm{OH}^{-}$ is required for $\mathrm{Zn}$ dissolution in reaction (I), the anode overpotential increases. Finally at the end of the 10th cycle, the supply of $\mathrm{OH}^{-}$ to the $\mathrm{Zn}$ metal surface is insufficient to sustain the cell current and the cell voltage diverges. Without carbon dioxide in the feed gas, electrolyte stability is improved. The slow continuous decrease in efficiency due to a reshaping of the anode is discussed in the subsequent Section 5.3.

In order to evaluate long-term battery operation, we simulate cells with synthetic air without $\mathrm{CO}_{2}$ in the gas feed in the remainder of this article.

5.3 Simulation of Cycle-Life at Various Current Densities

Subsequent to model validation, we simulate repeated discharge and charge cycles of zinc-air cells for various current densities and analyze the zinc shape-change effect. Because ambient air with $\mathrm{CO}_{2}$ limits battery life-time [11], we assume synthetic air without $\mathrm{CO}_{2}$ as feed gas. The initial $\mathrm{Zn}$ and $\mathrm{ZnO}$ distribution in the cell is visualized in Figs. 4a,d,g. It constitutes our reference for a quantification of the shape-change effect. The applied galvanostatic current densities range from 25\text{\,}\mathrm{A}\text{{\cdot}}{\mathrm{m}}^{-2} to 175\text{\,}\mathrm{A}\text{{\cdot}}{\mathrm{m}}^{-2}. We cycle the cells between $21\text{\,}\mathrm{\char 37\relax}$ and $90\text{\,}\mathrm{\char 37\relax}$ state-of-charge.

We show in Fig. 5a-c that cycle-life is shorter if greater current densities are applied. Note that available computation time limits the simulation at the lowest current density 25\text{\,}\mathrm{A}\text{{\cdot}}{\mathrm{m}}^{-2} and we do not observe the full cycle-life. At the end-of-life, cell voltage diverges. During a single discharge, cell voltage generally decreases apart from the nucleation dip; during a single charge, cell voltage is U-shaped. The overpotential for charge and discharge is increasing from cycle to cycle.

The cell voltage profile during a single discharging is explained above in Section 5.1. During charging, two concurrent effects result in the U-shaped voltage profile. On the one hand, deposition of $\mathrm{Zn}$ increases $\mathrm{Zn}$ volume fraction and specific surface area. This decreases the cell potential. On the other hand, dissolution of $\mathrm{ZnO}$ decreases $\mathrm{ZnO}$ volume fraction and specific surface area. Thus, $\mathrm{Zn(OH)}^{2-}_{4}$ concentration in the electrolyte decreases. This increases the cell potential. During repeated discharge and charge cycles, solid phases move towards the lid and agglomerate as discussed below. Therefore, the increase in cell voltage towards the end-of-charge becomes more prominent from cycle to cycle. Cycle-life is ultimately limited by agglomeration of solid phases and depletion of $\mathrm{Zn(OH)}^{2-}_{4}$ in the electrolyte.

The phase composition in the fully discharged state just before the end-of-life is shown in Fig. 5d-i. We observe that $\mathrm{Zn}$ and $\mathrm{ZnO}$ are shifted away from the cathode. For low current densities, the solids accumulate at the lid. For the greatest examined current density 175\text{\,}\mathrm{A}\text{{\cdot}}{\mathrm{m}}^{-2}, they do not reach the lid.

The loss of efficiency and limit of cycle-life is a consequence of the movement and reorganization of the solid phases. Firstly, the farther the solids move away from the cathode, the longer the transport paths in the electrolyte become. This increases the transport losses. The second, even more important, reason is the shrinkage of the total active surface area. Agglomeration of $\mathrm{ZnO}$ reduces the total surface area, which implies higher overpotentials. Eventually, the solids hit the housing. Finally, the total $\mathrm{ZnO}$ surface becomes too small to provide enough zincate ions during charging, as required by reaction I. Therefore, the zincate concentration depletes, leading to a diverging cell potential. For low enough current density, the zincate concentration does not deplete. The greater the current density, the earlier $\mathrm{Zn(OH)}^{2-}_{4}$ depletion is reached (see Fig. 5). For the large current density 175\text{\,}\mathrm{A}\text{{\cdot}}{\mathrm{m}}^{-2}, the shrinking $\mathrm{ZnO}$ surface area causes a zincate depletion already after three cycles, even though the solid phases do not touch the lid.

The observed zinc shape-change results from the relation of $\mathrm{Zn(OH)}^{2-}_{4}$ concentration and solubility limit $c_{\mathrm{sat}}$ as shown in Fig. 6a. The $\mathrm{Zn(OH)}^{2-}_{4}$ saturation concentration is inhomogeneous in through direction. Its profile is steeper for larger current densities. The sign of the gradient is different during charging and discharging. The saturation concentration is determined by the the potassium concentration in the electrolyte (see LABEL:eq:suppl:params:reac:II:c_sat). During discharging, the solubility is lower at the lid than at the cathode and $\mathrm{Zn(OH)}^{2-}_{4}$ accumulates at the lid due to the electric field. This leads to faster precipitation of $\mathrm{ZnO}$ near the lid (see reaction II). During charging, the saturation concentration is larger at the lid than at the cathode. A large solubility results in faster dissolution, providing more $\mathrm{Zn(OH)}^{2-}_{4}$ near the lid. Thus, the $\mathrm{Zn}$ deposition rate at the lid increases (see reaction I). The combination of both mechanisms tends to shift the solid phases away from the cathode, resulting in an overall $\mathrm{Zn}$ and $\mathrm{ZnO}$ movement during cycling.

LABEL:fig:suppl:misc:movement in the Supplementary Materials depicts the reaction rate of reaction (I). Positive values correspond to dissolution and negative values correspond to formation of $\mathrm{Zn}$. We observe that $\mathrm{Zn}$ dissolves and forms predominantly near the cathode. Comparing dissolution and deposition, however, $\mathrm{Zn}$ formation occurs slightly closer to the lid than $\mathrm{Zn}$ dissolution. Therefore, charging and discharging do not exactly balance each other and a phase-shift to the lid remains.

In order to examine the movement characteristics, we plot the center-of-volume for the solid phase in Fig. 6b as defined in LABEL:eq:position. It shows the central position of the active solid phases over time for different galvanostatic current densities. All simulations result in a zigzag course of the central position with a tendency towards the lid. Independent from the current, the center-of-volume remains above $\sim$1.3\text{,}\mathrm{mm}$$. Simulations with high current densities show a termination already before the $\mathrm{ZnO}$ reaches the cell housing. Higher current densities cause a faster shift of the solid phases in the same number of cycles.

An important observation is the zigzag course of the center-of-solid-volume. This implies that during charging and discharging the shift has partially opposite directions. However, it does not cancel each other, but the discharge outweighs with a movement towards the lid. We also find that the backwards shift during charging is more dominant for lower currents. Using this characteristic helps to reduce the shape-change with an adjusted charge-discharge protocol. A detailed analysis on this is given in Section 6.2.

The analysis of the central position of the active solid phases in Fig. 6b agrees with the results in Fig. 5. A higher current causes faster shape-change and reduced cycle-life. At some point, the movement is limited because the block of solids hits the lid. As shown in Fig. 5, the solids do not reach the cell housing for high current densities. In this case, agglomeration occurs in the center of the cell as shown in the Supplementary Materials in LABEL:fig:suppl:currentVarConcSat.

We summarize that $\mathrm{Zn}$ electrode shape-change limits cycle-life of zinc-air batteries and that lower cell currents significantly reduce shape-change.

6 Discussion of Cycle-Life

We discuss strategies to reduce the zinc shape-change effect and improve cycle-life. In Section 6.1, we evaluate various $\mathrm{Zn}$ anode morphologies including 3D current collectors and inhomogeneous $\mathrm{Zn}$ distributions. The cycling protocol has a very strong influence on cycle-life as discussed in Section 6.2. Finally, in Section 6.3, we analyze how better reaction and transport kinetics could improve cycle-life.

6.1 Zinc Electrode Design Variation

Here, we discuss how variations of the initial distribution of the solids in the anode affect cycle-life. We use the Renata ZA13 cell [44] as reference. Its initial phase distribution is shown in Figs. 4a,d,g. We simulate 23 cycles for a galvanostatic cycling current density of 100\text{\,}\mathrm{A}\text{{\cdot}}{\mathrm{m}}^{-2} (see Fig. 5b). We add $2\text{ volume percent}$ $\mathrm{ZnO}$ to the $\mathrm{Zn}$ paste for avoiding inhomogeneous $\mathrm{ZnO}$ nucleation as demonstrated in Section 5.1.

Firstly, we assume a homogeneous initial $\mathrm{Zn}$ distribution and achieve 17 cycles (see Figs. 7a-c). The solids cannot move freely without the gas reservoir below the lid which reduces cycle-life. In a second step, we add a supporting structure or 3D current collector in order to suppress solid convection (see Figs. 7d-i). To model this feature, we turn off solid convection (see LABEL:sec:suppl:params:levs). For simplicity, we assume that this structure consists of inert $\mathrm{Zn}$ material with a volume fraction of $1.92\text{\,}\mathrm{\char 37\relax}$. Our ideal supporting structure guarantees a minimum surface area for $\mathrm{Zn}$ deposition and $\mathrm{Zn}$ precipitation throughout the cell and an electron pathway in the solids through the electrode. Thus, we promote the more homogeneous formation of solid phases as shown in Figs. 7c,f,i. Nevertheless, cycle-life is even reduced slightly in this case. We note that our simple model for solid convection is agnostic of the particle nature of the $\mathrm{Zn}$ electrode and tends to overestimate cycle-life. Future battery models could make use of hydrodynamic models for granular flow and include effects such as shear bands and anisotropic pressure distributions [86]. In Figs. 7g-i, we use a $\mathrm{Zn}$ distribution increasing linearly towards the cathode. The adjusted initial profile counteracts zinc shape-change and improves cycle-life up to 19 cycles.

6.2 Charge-Discharge Protocol Optimization

Smart charge-discharge protocols can improve cycle-life. Here, we discuss the effect of cycling range, applied charge/discharge currents, and recovery cycles. Table 2 summarizes the applied protocols. The initial $\mathrm{Zn}$ distribution for all simulations is depicted in Fig. 7h. As baseline reference, we use protocol A, whose simulated cycling behavior is shown in Figs. 7g-i. We show simulated cell voltages during cycling in Fig. 8a and summarize the cycle-life in Table 2. The protocols B-E show a positive effect on the maximum number of cycles. However, not all extend the total utilizable capacity.

Introducing recovery cycles B, in which the cell is charged and discharged very slowly, provides no advantage. This protocol increases the number of maximum cycles exactly by the number of recovery cycles yielding the same utilizable capacity. Since recovery cycles take extra time and waste energy they are even disadvantageous.

Cycling a smaller capacity range increases the number of cycles by a factor of four (see protocol C). Due to the smaller range, however, each discharge provides less capacity. In total, the utilizable capacity increases by $15\text{\,}\mathrm{\char 37\relax}$. Nevertheless, this protocol has a major drawback. It utilizes only $20\text{\,}\mathrm{\char 37\relax}$ of the total capacity and reduces the total energy density by a factor of four.

Reducing the discharging current to 25\text{\,}\mathrm{A}\text{{\cdot}}{\mathrm{m}}^{-2} in D increases the utilizable capacity and the number of cycles by $66\text{\,}\mathrm{\char 37\relax}$. However, it is not always possible to reduce the discharging current, as it depends on the power specifications of the application, in which the battery is used.

The most impressive protocol is E. By reducing the current density for charging, we increase the number of cycles and the total capacity by more than twenty-fold. This shows that cycling zinc-air batteries is limited by the charging current and that low charging currents are beneficial for cycle-life.

Finally, we discuss the movement of the center-of-solid-volume for above protocols. Fig. 8b shows the influence of the protocol on the movement of the solid phases. The most distinct shape-change effect is observed for the protocols A and C, whereas protocol E shows the smallest one. Very small capacity ranges are unfavorable as the inert electrolyte can prevent the back-shift of solids during a short charging period (see protocol C). Recovery cycles in protocol B indeed reduce the movement, but have no further special benefit. The best results are achieved for low charging currents as in protocol E due to the solid movement characteristics during charging and discharging. A large shift occurs during discharge, but the back-shift during charge compensates this shift for low charging currents. This demonstrates that the dissolution of $\mathrm{ZnO}$ is the critical process creating inhomogeneous $\mathrm{Zn}$ deposition and irreversible agglomeration.

6.3 Reaction and Transport Kinetics Optimization

Adjusting the intrinsic reaction and transport kinetics offers further options for improving cycling performance. Experimentally, this can for example be achieved by electrolyte additives, which promote $\mathrm{ZnO}$ dissolution, or additives to the $\mathrm{Zn}$ paste, which promote homogeneous $\mathrm{Zn}$ deposition. In the Supplementary Materials in LABEL:sec:suppl:ic, we discuss this with two simulation approaches. Designing thicker anodes leads to a larger total $\mathrm{ZnO}$ surface area, improves $\mathrm{ZnO}$ dissolution, and avoids early $\mathrm{Zn(OH)}^{2-}_{4}$ depletion. By doubling the anode thickness, we achieve $1.4$ times more cycles. Increasing the temperature to {T=343.15\text{,}\mathrm{K}} enhances the transport processes in the electrolyte. This makes it possible to achieve cycle numbers similar to protocols D and E without decreasing the charging current (see Section 6.2).

7 Conclusion and Outlook

In this work, we present the first multi-dimensional model for conversion type batteries and apply it to zinc-air button cells. Local volume-averaging is the theoretical concept that enables us to model dissolution, growth, and movement of the porous zinc structure. Thus, we develop the most in-depth simulation study of inhomogeneous metal electrode structures so far.

Our model is successfully validated with x-ray tomography of commercial zinc-air cells. We demonstrate how electrolyte carbonation limits shelf-life and how zinc oxide nucleation limits rechargeability. The strong dependence of cycle-life on cell current is a major model prediction. In order to optimize cycling performance, we study various anode structures and examine the impact of charge-discharge protocols. Our simulations show that smaller cycling ranges and recovery cycles yield minor improvements. The best results are achieved for low current densities, especially low charging currents. Thus, dissolution of zinc oxide seems to be the critical process for cycle-life. As a consequence, accelerated intrinsic reaction and transport kinetics would enhance cycle-life as demonstrated by thicker anodes and higher operating temperature.

One-dimensional simulations, as by Stamm et al. [11], remain relevant because the radial distribution of zinc and zinc oxide in coin cells has little influence on our simulation results. Nevertheless, it is well established in the field of lithium-ion batteries that multi-dimensional simulations are indispensable for studying degradation in inhomogeneous electrodes, for example, lithium plating [39]. Thus, our novel simulation methodology paves the way for future insights into degradation of conversion batteries, for example, due to hydrogen evolution.

Conflicts of Interest

There are no conflicts to declare.

Acknowledgments

The authors thank Simon Hein and Max Schammer for fruitful discussions. This work was supported by the German Ministry of Education and Research (BMBF) (project LUZI, BMBF: 03SF0499E). Further support was provided by the bwHPC initiative and the bwHPCC5 project through associated compute services of the JUSTUS HPC facility at the University of Ulm. This work contributes to the research performed at CELEST (Center for Electrochemical Energy Storage Ulm-Karlsruhe).