Analysis of the Effect of Operating Conditions on the Performance of a Direct Ammonia Fuel Cell Using Multiphysics Modelling

Zhao et al. ²³ assumed, and provided supporting evidence, that the rate-determining step (RDS) of the AOR was the recombination of two adsorbed NH₂ molecules to N₂H₄ on the catalyst surface

$$\begin{eqnarray}&&2{}^{\ast }{\rm{N}}{{\rm{H}}}_{2}\rightleftharpoons {}^{\ast }{\rm{N}}_{2}{{\rm{H}}}_{4}\end{eqnarray} \tag{ 4 }$$

Taking into account the surface coverage of NH₂, this led to a current-overpotential relationship that differs from the more common Butler-Volmer or Tafel equations, the derived expression being

$$\begin{eqnarray}&&i=\frac{A{\left[{\text{NH}}_{3}\right]}^{2}{e}^{\frac{2{nF}}{{RT}}\left({\varphi }_{a}-{\varphi }_{1}^{0}\right)}}{{\left(1+\left[{\text{NH}}_{3}\right]{e}^{\frac{{nF}}{{RT}}\left({\varphi }_{a}-{\varphi }_{1}^{0}\right)}\right)}^{2}}\end{eqnarray} \tag{ 5 }$$

The equilibrium potential of ^*NH₂ is ${\varphi }_{1}^{0},$ the applied potential is ${\varphi }_{a}$ , A is a lumped kinetic parameter that corresponds to the maximum AOR activity, and n (=1) is the number of electrons transferred in the oxidation of NH₃ to NH₂, which is the previous reaction step that is assumed to be in Nernst equilibrium. The ^*NH₂ equilibrium potential is about 0.40–0.45 V vs RHE, decreasing linearly with increasing temperature (values at 60 °C and 95 °C specified in Table S1), corresponding to about 0.4 V positive of the AOR theoretical potential. Zhao et al. ²³ found a temperature coefficient of ca. −0.580 mV K⁻¹ (168 J mol⁻¹ K⁻¹ entropy with n = 3) to best fit their AOR kinetic data.

$$\begin{eqnarray}&&{\left(\frac{\partial \varphi }{\partial T}\right)}_{p}=\frac{-{\rm{\Delta }}S}{{nF}}\end{eqnarray} \tag{ 6 }$$

$$\begin{eqnarray}&&{\varphi }_{1}^{0}\left(T\right)={\varphi }_{1}^{0}\left({T}^{0}\right)-0.000580\frac{V}{K}\left(T-{T}^{0}\right)\end{eqnarray} \tag{ 7 }$$

Here the reference temperature T⁰ is 95 °C and the potential at the reference temperature ${\varphi }_{1}^{0}\left({T}^{0}\right)$ is 0.416 V vs RHE. The potential as a function of temperature is plotted together with the other theoretical potentials in the next section.

The kinetics given by Eq. 5 at 75 °C with our parameters (we later use 75 °C as baseline temperature for our DAFC simulations) is shown below in Fig. 1, where the vertical dashed line indicates ${\varphi }_{1}^{0}.$ At potentials more negative than the ^*NH₂ potential, the jV curve behaves similarly to the Tafel equation with a charge transfer coefficient of 2, i.e., about 30 mV/decade Tafel slope at room temperature, but eventually the current density saturates to A at more positive potentials.

Zoom In Zoom Out Reset image size

While Zhao et al. ²³ provided the necessary parameter values at 60 °C and at 95 °C, the temperature dependency was given for only some parameters. For kinetics (${j}_{0,\text{ORR}}$ and A), the Arrhenius equation was assumed, and a linear temperature dependency was shown for the ^*NH₂ potential ${\varphi }_{1}^{0}.$ For the ORR, we used the given parameters directly. However, for the AOR we noticed that a numerical value given in a table (in supporting information, 35 A g⁻¹ at 60 °C) did not quite match the kinetics data shown in the figures, most likely due to a typographical error, so we extracted the value and activation energy from the figure showing the catalytic activity parameter as a function of temperature (45 A g⁻¹ at 60 °C and 33.89 kJ mol⁻¹).

For the ORR, we use the theoretical voltage of water electrolysis together with the 0.9 V vs RHE potential to extrapolate the exchange current density of the ORR. A linear fit to the NIST-JANAF thermochemical tables ³⁹ gave for the ORR

$$\begin{eqnarray}&&\begin{array}{c}{E}_{\text{ORR}}^{0}\left(T\right)={V}_{{\text{H}}_{2}\text{O}}^{0}\left(T\right)\\ =1.2289\text{V}-0.00083465\frac{\text{V}}{\text{K}}\left(T-298.15\text{K}\right)\end{array}\end{eqnarray} \tag{ 8 }$$

This corresponds also to the theoretical voltage of a hydrogen fuel cell (or a water electrolyser). A similar linear fit yields the theoretical AOR potential as follows

$$\begin{eqnarray}&&{E}_{\text{AOR}}^{0}\left(T\right)=0.05654\text{V}-0.00035374\frac{\text{V}}{\text{K}}\left(T-298.15\text{K}\right)\end{eqnarray} \tag{ 9 }$$

Therefore, the theoretical DAFC voltage as the potential difference of the half reactions is

$$\begin{eqnarray}&&{V}_{\text{DAFC}}\left(T\right)=1.17235\text{V}-0.0004809\frac{\text{V}}{\text{K}}\left(T-298.15\text{K}\right)\end{eqnarray} \tag{ 10 }$$

A linear fit directly to the thermodynamic data of the total reaction yielded a temperature slope whose magnitude was about 2.5 $\mu$V K⁻¹ smaller, so we consider this difference negligible. The voltage and potentials are shown in Fig. 2, including also the ^*NH₂ potential ${\varphi }_{1}^{0}$ (Eq. 7).

Zoom In Zoom Out Reset image size

For the ionic conductivity, based on the measurements of Schwämmlein et al., ⁴⁰ we assumed that the product of conductivity and temperature follows the Arrhenius equation.

$$\begin{eqnarray}&&\sigma \left(T\right)T=\sigma \left({T}^{0}\right){T}^{0}{e}^{-\frac{{E}_{\text{A}}}{R}\left(\frac{1}{T}-\frac{1}{{T}^{0}}\right)}\end{eqnarray} \tag{ 11 }$$

Fitting the values for the AEM and the catalyst layer conductivity yielded activation energies of 16.6 kJ mol⁻¹ and 26.8 kJ mol⁻¹, respectively. Also in this case, the reference temperature T⁰ is 95 °C.

The parameter values given by Zhao et al. ²³ that we used are given in Table S1. For the AOR and ORR standard potentials we used Eqs. 8 and 9, and the ORR and AOR activities are assumed to follow the Arrhenius equation.

Our model will be discussed in more detail in the following sections, and here we give a brief, general description of the cell operation together with a scheme of the cell cross-section, including the most important physical phenomena, couplings, and boundary conditions.

The most important physical phenomena are shown in Fig. 4 together with the associated boundary conditions. The CLs, where the electrochemical reactions occur, involve the largest number of phenomena and their interactions: Electric conductivity is needed to supply electrons to and extract them from the electrode, and hydroxide ions (and water) and gas need to be moved to and from the catalysts in the ionomer and empty space, respectively. The reactions, AOR and ORR, convert electric and ionic current to each other and couple the different forms of transport to each other. The AEM between the electrodes is responsible for connecting the reactions to each other through ionic transport, while blocking electron transport through the cell. In CCMs the gas that is transported into and out of the CLs proceeds through the PTLs, as well as the electron transport in the electrically conductive substrate (electric resistance, ohmic losses). As mentioned earlier, ammonia crossover through the AEM is not included in the model.

Zoom In Zoom Out Reset image size

The contribution of the partial pressures was added to the theoretical AOR and ORR potentials using the Nernst equation

$$\begin{eqnarray}&&{E}_{\text{AOR}}={E}_{\text{AOR}}^{0}\left(T\right)-\frac{{RT}}{6F}\mathrm{ln}\left(\frac{{p}_{{\text{NH}}_{3}}^{2}}{{p}_{{\text{N}}_{2}}{p}^{0}}\cdot {a}_{\text{w},\text{m}}^{-6}\right)\end{eqnarray} \tag{ 12 }$$

$$\begin{eqnarray}&&{E}_{\text{ORR}}={E}_{\text{ORR}}^{0}\left(T\right)+\frac{{RT}}{4F}\mathrm{ln}\left({a}_{\text{w},\text{m}}^{2}\cdot \frac{{p}_{{\text{O}}_{2}}}{{p}^{0}}\right)\end{eqnarray} \tag{ 13 }$$

We assumed that water absorbed in the ionomer participated in the reaction, so the water activity was the ionomer water activity a_w,m (defined later, Eq. 19), which was assumed to be in equilibrium with the gas phase humidity. The reference pressure for O₂, N₂ and NH₃ is ${p}^{0}=1\text{bar}.$ The activity of the hydroxide ions in the ionomer and the AEM was assumed to be constant 1 and therefore it did not affect the potentials. The temperature-dependent reference potentials are given by Eqs. 8 and 9.

The main expression of the AOR kinetics is given by Eq. 5 and the effects of catalyst loading and gas transport amounted to varying the catalytic activity parameter A (or rather m_catalyst A/L like in Eq. S2) and the ammonia concentration, respectively. Catalytic activity was assumed to depend linearly on the catalyst loading (per thickness or volume). The catalytic activity and the ^*NH₂ reference potential ${\varphi }_{1}^{0}$ depended on temperature, as specified earlier in section "Literature model."

For the ORR, we extended the kinetics by assuming that the exchange current density depends on the square root of the O₂ partial pressure, based on Neyerlin et al. ⁴⁵

$$\begin{eqnarray}&&{i}_{0}\left(T,p\right)={i}_{0}^{{\prime} }\left(T,{p}^{* }\right){\left(\frac{{p}_{{\text{O}}_{2}}}{{p}^{* }}\right)}^{0.5}\end{eqnarray} \tag{ 14 }$$

Here the reference pressure $\left({p}^{* }\right)$ corresponds to 3 bar absolute pressure, because Zhao et al. applied 2 bar gauge pressure to their cathode. ²³ The local exchange current density depends on the catalyst loading and CL thickness (L_CL), similarly to AOR and Eq. S2.

We use effective medium approximation to describe the electrodes. For most properties, we used the Bruggeman approximation

$$\begin{eqnarray}&&{X}_{\text{eff}}={X}_{b}{\varepsilon }^{1.5}\end{eqnarray} \tag{ 15 }$$

The bulk property is marked with subscript b and the power 1.5 of the volume fraction $\varepsilon$ includes both the linear cross-section area reduction ($\varepsilon$) and the consequently increased tortuosity that reduces e.g., ionic conductivity (${\varepsilon }^{0.5}$). However, with the ionomer volume, as discussed later in this section, the main parameter varied in simulations was the fraction of the empty volume left after the catalyst and substrate, and not the fraction of the total volume.

The solid volume fraction of both the SSF and the carbon paper substrates is 20%, ${\varepsilon }_{s}=0.2,$ and in the CL of the CCM the carbon support of the catalyst occupies 25% of the volume, ${\varepsilon }_{s}=0.25.$ Carbon paper is based on ca. 2 g cm⁻³ solid density (amorphous carbon 1.8 − 2.1 g cm⁻³⁴⁶ ), and on ca. 200 $\mu$m thickness and 90 g m⁻² areal weight ⁴⁷ yielding 22.5% solid fraction, rounded down. The SSF properties are based on typical steel felt characteristics, e.g. ⁴⁸ For the carbon support, using the anode properties of Zhao et al. ²³ (4 mg_Pt cm⁻² of PtIr, 340 $\mu$m thickness), 2 g cm⁻³ carbon density, and assuming that 20% of the total catalyst weight is metal, i.e. 80% carbon, yields ca. 24% carbon volume fraction. We rounded this up to 25% and assumed the same volume fraction for the cathode because a similar volume fraction could be produced with reasonable numbers. (At the cathode, 3 mg cm⁻² of catalyst and 120 μm thickness with 30% of catalyst weight being metal would give a 29% carbon volume fraction. Without knowing the specific values, using 25% also for the cathode therefore seemed reasonable.)

When simulating the CCS, the catalyst has no additional support. The catalyst volume fraction depends on the layer thickness L, catalyst loading m_catalyst, and the density of the catalyst material ${\rho }_{\text{catalyst}}$

$$\begin{eqnarray}&&{\varepsilon }_{\text{catalyst}}=\frac{{m}_{\text{catalyst}}}{L{\rho }_{\text{catalyst}}}\end{eqnarray} \tag{ 16 }$$

For PtIr we use the average of the densities of Pt (21.4 g cm⁻³) and Ir (22.56 g cm⁻³), 21.98 g cm⁻³, and for the Acta 4020, an FeCo-based catalyst, ⁴⁹ the density of cobalt iron alloy, 8.6 g cm⁻³. ⁵⁰ As mentioned earlier, in the case of CCM, we assume that the catalyst loading per volume does not change, but the layer thickness depends on the loading, meaning that the catalyst volume fraction is constant (∼0.54% at anode and ∼2.9% at cathode, based on the catalyst loadings and layer thicknesses from Zhao et al. ²³ ), the baseline loading of 4 mg cm⁻² of catalyst corresponding to a thickness of 340 $\mu \text{m}$ at anode and to 120 $\mu \text{m}$ at cathode. In the case of the CCS and 0.5 mm thick SSF, the catalyst volume fraction depends on the catalyst loading, the baseline 4 mg cm⁻² corresponding to ∼0.36% at anode and to ∼0.93% at cathode.

The electric conductivities of the substrate, catalyst, and possible carbon support were added together for the total electric conductivity of the layer. For the metals, we used their bulk values and the Bruggeman approximation, and in the case of the stainless steel felt we also considered the effect of the temperature. The conductivity of the carbon support was based on the data for carbon black from Pantea et al. ⁵¹

For the ionomer volume fraction, we vary the fraction of the non-solid volume (${f}_{\varepsilon ,\text{ionomer}}$) in our simulations to ensure that that the sum of the volume fractions of all components is 1. The absolute ionomer volume fraction is

$$\begin{eqnarray}&&{\varepsilon }_{\text{ionomer}}={f}_{\varepsilon ,\text{ionomer}}\left(1-{\varepsilon }_{\text{s}}-{\varepsilon }_{\text{catalyst}}\right)\end{eqnarray} \tag{ 17 }$$

In practice, the absolute volume fraction is in most cases about 75–80% of the value of f_ε,ionomer. The remaining empty volume where gas is transported is

$$\begin{eqnarray}&&\begin{array}{c}{\varepsilon }_{\text{void}}=1-{\varepsilon }_{\text{s}}-{\varepsilon }_{\text{catalyst}}-{\varepsilon }_{\text{ionomer}}\\ =\left(1-{f}_{\varepsilon ,\text{ionomer}}\right)\left(1-{\varepsilon }_{\text{s}}-{\varepsilon }_{\text{catalyst}}\right)\end{array}\end{eqnarray} \tag{ 18 }$$

The main reason for this parameter choice was to simplify varying the parameter range by guaranteeing that the entire range from 0 to 1 would be feasible for f_ε,ionomer in all cases, meaning that the remaining empty volume fraction ε_void would always be positive.

In this section we discuss our model for the AEM and ionomer, including the water uptake and transport, ionic conductivity and the ionomer volume fractions in the electrode layers.

The water uptake and ionic conductivity of the ionomer and the AEM are based on the parametrization given by Gerhardt et al. ⁵² and we only modify the ionic conductivity by reducing it by constant multipliers to match the values to the ones given by Zhao et al. ²³ The water uptake is

$$\begin{eqnarray}&&\begin{array}{c}\lambda =\left(-0.6{a}_{\text{w},\text{ionomer}}^{3}+0.85{a}_{\text{w},\text{ionomer}}^{2}-0.2{a}_{\text{w},\text{ionomer}}+0.153\right)\left(T-313\right)\\ +39{a}_{\text{w},\text{ionomer}}^{3}-47.7{a}_{\text{w},\text{ionomer}}^{2}+23.4{a}_{\text{w},\text{ionomer}}+0.117\end{array}\end{eqnarray} \tag{ 19 }$$

In practice, we calculate the uptake from the water concentration (and ionomer density and equivalent weight, EW) and then solve the water activity using the formula for the cubic equation (see supporting information for details). ⁵³ From the water volume fraction (f, Eq. 38 in Weber and Newman ⁵⁴ ) the following expression for the water uptake can be derived, by recognizing that the water concentration is equal to the water volume fraction divided by the molar volume of liquid water (c_W = f/V_M,W)

$$\begin{eqnarray}&&\lambda =\frac{{c}_{\text{w}}\cdot {EW}}{{\rho }_{\text{ionomer}}\bullet \left(1-{c}_{W}{V}_{M,W}\right)}\end{eqnarray} \tag{ 20 }$$

The second term in the parentheses corresponds to membrane swelling, and the molar volume is equal to the molar mass divided by density, V_M,W(T) = M_H2O/ρ_H2O(T) (M_H2O = 18.015 g/mol). For the density of liquid water, we used the data from the CRC Handbook of Chemistry and Physics. ⁴⁶ For the EW we use 588.24 g mol⁻¹ and the density is 1.1 g cm⁻³, both corresponding to Tokuyama A201 AEM. ^52,55

The basis for the ionic conductivity, that we reduce by constant multipliers in both AEM and ionomer, is

$$\begin{eqnarray}&&\begin{array}{c}{\kappa }_{{\text{OH}}^{-}}\\ =\frac{20}{100}\left(0.1334-0.0003882T+\left(0.01148T-3.9\right){a}_{\text{w},\text{ionomer}}-\left(0.0669T-23\right){a}_{\text{w},\text{ionomer}}^{2}+\right.\\ \left.\left(0.1227T-42.61\right){a}_{\text{w},\text{ionomer}}^{3}-\left(0.06T-21.8\right){a}_{\text{w},\text{ionomer}}^{4}\right)\end{array}\end{eqnarray} \tag{ 21 }$$

The temperature is in Kelvins and the unit of the conductivity is S cm⁻¹. In the catalyst layer, the effective conductivity is calculated using the Bruggeman model

$$\begin{eqnarray}&&{\kappa }_{{\text{OH}}^{-},\text{eff}}={\varepsilon }_{\text{ionomer}}^{1.5}{\kappa }_{{\text{OH}}^{-}}\end{eqnarray} \tag{ 22 }$$

To match the AEM of Zhao et al. ²³ with ${a}_{\text{w},\text{m}}=1,$ we divide the conductivity by 11.73, which is the average of the multipliers that match their measured values at 60 °C and at 95 °C (${\kappa }_{\text{AEM}}={\kappa }_{{\text{OH}}^{-}}/11.73$). The match is of course not perfect because the temperature dependencies of this parametrized model and the reported values reported by Zhao et al. ²³ differ from each other, but this model forms a good basis for simulating the effect of the DAFC operating conditions on the AEM and ionomer. For the ionomer, assuming ${\varepsilon }_{\text{ionomer}}\approx 0.64$ (substrate volume fraction being ${\varepsilon }_{\text{s}}=0.25$ and ${f}_{\varepsilon ,{ionomer}}=0.85,$ i.e. ionomer occupying 85% of the non-solid volume), the best match at 95 °C would be achieved by dividing the conductivity by 5.23 and lower conductivities improve the match at lower temperature. For most of our simulations, we divided the ionomer conductivity by 7, because this value provided a close match near our 75 °C baseline temperature (as discussed later in "Results"). Here we are matching the product of two model parameters to one measured value, meaning in practice that one parameter needs to be fixed (ionomer volume) before the other (conductivity) is varied for the best fit. If the resulting layer conductivity does not change, the individual parameters values are not important for most of the ionomer volume fraction range (discussed later in "Results—Ionomer Volume Fraction and Gas Transport").

We model the water flux in the ionomer as a sum of the diffusional flux and the electro-osmotic drag. However, due to software limitations, we needed to adapt the electro-osmotic drag term to a form of drift in electric field, which we explain in this section. In steady state

$$\begin{eqnarray}&&{\rm{\nabla }}\cdot \vec{{j}_{\text{w}}}={R}_{\text{w}}\end{eqnarray} \tag{ 23 }$$

And the water flux is

$$\begin{eqnarray}&&\vec{{j}_{\text{w}}}=-{D}_{\text{w}}{\rm{\nabla }}{c}_{\text{w}}-z\mu F{c}_{\text{w}}{\rm{\nabla }}{\varphi }_{\text{l}}\end{eqnarray} \tag{ 24 }$$

The bulk diffusion coefficient is the parameterization given by Gerhardt et al. ⁵² without the term for the fraction of the bicarbonate form of the ionomer, and in the ionomer we apply the Bruggeman model to account for the ionomer volume fraction.

$$\begin{eqnarray}&&\begin{array}{c}{D}_{{\text{H}}_{2}\text{O},b}\\ =2.07\cdot {10}^{-7}\frac{{\text{cm}}^{2}}{\text{s}}\exp \left(3.95{a}_{\text{w},\text{ionomer}}\right)\cdot \exp \left(\frac{17.7\frac{\text{kJ}}{\text{mol}}}{R}\left(\frac{1}{333}-\frac{1}{T}\right)\right)\end{array}\end{eqnarray} \tag{ 25 }$$

$$\begin{eqnarray}&&{D}_{{\text{H}}_{2}\text{O}}={\varepsilon }_{\text{ionomer}}^{1.5}{D}_{{\text{H}}_{2}\text{O},b}\end{eqnarray} \tag{ 26 }$$

As mentioned, the water activity is solved as the inverse of Eq. 19.

The electro-osmotic drag of the water molecules in the ionomer and AEM is included in our model. We assume that the hydroxide ions drag water molecules in the same direction with them and that the fluxes are equal to each other. ^52,56

$$\begin{eqnarray}&&{\vec{j}}_{\text{w},\text{drag}}=-\frac{\vec{i}}{F}\end{eqnarray} \tag{ 27 }$$

$$\begin{eqnarray}&&\vec{i}=-{\sigma }_{\text{l}}{\rm{\nabla }}{\varphi }_{\text{l}}\end{eqnarray} \tag{ 28 }$$

The ionomer/AEM conductivity (${\sigma }_{\text{l}}$) is the conductivity determined by the cell temperature and water activity. In practice, we implement this using the "Migration in electric field" option in the Transport of Diluted Species node of Comsol. The electronic/ionic current is driven by the electrolyte potential gradient (${\rm{\nabla }}{\varphi }_{\text{l}}$), so, for the water flux to match the hydroxide, the water mobility (${\mu }_{\text{w}}$) needs to be expressed as a function of the ionic conductivity and water concentration.

$$\begin{eqnarray}&&{\vec{j}}_{\text{w},\text{drag}}=-{z}_{\text{w}}{\mu }_{\text{w}}F{c}_{\text{w}}{\rm{\nabla }}{\varphi }_{\text{l}}\end{eqnarray} \tag{ 29 }$$

From this and Eq. 28, an expression for the mobility follows, setting ${z}_{\text{w}}=-1$ and ${\sigma }_{\text{l}}={\kappa }_{{\text{OH}}^{-}}$ (Eq. 21)

$$\begin{eqnarray}&&{\mu }_{\text{w}}=\frac{{\kappa }_{{\text{OH}}^{-}}}{{F}^{2}{c}_{\text{w}}}\end{eqnarray} \tag{ 30 }$$

This is treated as a variable, meaning that the value is calculated for each point in the cross-section with the corresponding water uptake and ionomer conductivity. In the catalyst layers both the ionic current density (conductivity) and water mobility are corrected with the volume fraction of the ionomer according to the Bruggeman model

$$\begin{eqnarray}&&{\mu }_{\text{w},\text{eff}}={\mu }_{\text{w}}{\varepsilon }_{\text{ionomer}}^{1.5}\end{eqnarray} \tag{ 31 }$$

Both AOR and ORR are assumed to produce or consume water, respectively, in the ionomer, corresponding to a source term

$$\begin{eqnarray}&&R=-\frac{\nu {i}_{\text{V}}}{{nF}}\end{eqnarray} \tag{ 32 }$$

The current density i_V is the reaction rate per volume, in A m⁻³, and the source term thus has the dimension mol m⁻³. For AOR, the number of electrons is $n=6$ and the stoichiometric coefficient $\nu =-6,$ and for the ORR, $n=4$ and $\nu =-2$ (see reactions (2) and (3)). In fuel cell operation the current density of the AOR is positive (oxidation), hence this source term is also positive, corresponding to water generation. For the ORR, the opposite applies.

The water activities in the ionomer and in the gas are assumed to be in equilibrium, which we in practice enforce with a kinetic source term with a high rate constant

$$\begin{eqnarray}&&{R}_{\text{m},\text{gas}}=-K\left({a}_{\text{w},\text{ionomer}}-{a}_{\text{w},\text{gas}}\right)\end{eqnarray} \tag{ 33 }$$

For the rate constant we use 10⁶ mol m⁻³ s⁻¹ as a large value to enforce equilibrium of water activities in the ionomer and in the gas. Negative values correspond to water evaporating from the ionomer to the gas and positive to the ionomer absorbing humidity from the gas. We assume that the gas phase activity is equal to the relative humidity, the ratio of the water partial pressure to the vapor pressure.

$$\begin{eqnarray}&&{a}_{\text{w},\text{gas}}={\rm{RH}}=\frac{{p}_{{\text{H}}_{2}\text{O}}}{{p}_{\text{w}}}\end{eqnarray} \tag{ 34 }$$

We calculate the flow of the gas mixture in the porous layers using Darcy's law, meaning that the gas velocity ($\vec{u}$) depends on the pressure gradient, dynamic viscosity of the gas ($\mu$) and the permeability ($\kappa$) of the porous layer.

$$\begin{eqnarray}&&{\rm{\nabla }}\cdot \left(\rho \vec{u}\right)={Q}_{\text{m}}\end{eqnarray} \tag{ 35 }$$

$$\begin{eqnarray}&&\vec{u}=-\frac{\kappa }{{\mu }_{v}}{\rm{\nabla }}p\end{eqnarray} \tag{ 36 }$$

The density of the gas is $\rho$ and Q_m is the mass source. The density $\rho$ and dynamic viscosity ${\mu }_{v}$ (subscript v to differentiate from ionic mobility) depend on the gas composition as will be defined later in this section.

For the permeability of the carbon paper we used the value $\kappa =1.0\cdot {10}^{-11}{\text{m}}^{2}$ based on Mangal et al. ⁵⁷ We calculated the permeability of the SSF using the Carman-Kozeny equation

$$\begin{eqnarray}&&\kappa =\frac{{\varepsilon }_{\text{void}}^{3}{d}^{2}}{16{k}_{\text{K}}{\left(1-{\varepsilon }_{\text{void}}\right)}^{2}}\end{eqnarray} \tag{ 37 }$$

We used the value for a fibrous medium, ${k}_{\text{K}}=6,$ ⁵⁸ for the geometry constant and 10 $\mu$m as the fibre diameter (d). In all layers, the void volume fraction was calculated with Eq. 18.

At the CL-AEM boundaries we set a no flow boundary condition for the gas mixture and fixed the pressure at the PTL-gas channel boundary to the total absolute pressure. This is the sum of the pressure of dry gas (N₂, O₂, NH₃) and the water pressure that corresponds to the set relative humidity at the cell temperature. The water partial pressure is calculated as the product of the RH and the vapor pressure given by Eq. S8.

$$\begin{eqnarray}&&{p}_{\text{A}}={p}_{\text{dry}}+\text{RH}{p}_{\text{w}}\end{eqnarray} \tag{ 38 }$$

The dry component consists of NH₃ and N₂ at the anode, and of O₂ and N₂ at the cathode. For both the anode and the cathode gas inflow, we define the molar fractions of these molecules as

$$\begin{eqnarray}&&{f}_{\text{k}}=\frac{{p}_{\text{k}}}{{p}_{\text{dry}}}\end{eqnarray} \tag{ 39 }$$

For the dynamic viscosity and binary diffusion coefficient we used Comsol's built-in values and models, selecting the Fuller-Schettler-Giddings model for the gas diffusivity and Wilke's model (without high pressure correction) for the gas viscosity. Wilke's model for the viscosity of a multi component system is ^59,60

$$\begin{eqnarray}&&{\mu }_{v}=\displaystyle \sum _{\text{j}}^{\text{n}}\frac{{x}_{\text{j}}{\mu }_{\text{j},\text{v}}}{\displaystyle \sum _{\text{i}}^{\text{n}}{x}_{\text{i}}{\varphi }_{\text{ji}}}\end{eqnarray} \tag{ 40 }$$

$$\begin{eqnarray}&&{\varphi }_{\text{ji}}=\frac{{\left(1+\sqrt{\frac{{\mu }_{\text{j},\text{v}}}{{\mu }_{\text{i},\text{v}}}}\cdot {\left(\frac{{M}_{\text{i}}}{{M}_{\text{j}}}\right)}^{1/4}\right)}^{2}}{\sqrt{8\left(1+\frac{{M}_{\text{j}}}{{M}_{\text{i}}}\right)}}\end{eqnarray} \tag{ 41 }$$

The mole fractions x_k depend on the diffusion of the species in the mixture, as discussed below. The vapor phase viscosity coefficients of all species (${\mu }_{\text{j},\text{v}}$) were set by Comsol Multiphysics and depend on the temperature. While not the precise expression used by Comsol, in our temperature range (20–95 °C) the viscosity coefficients increasing as T^0.76 is a reasonable approximation of the values used in the simulations (temperatures in Kelvins)

$$\begin{eqnarray}&&{\mu }_{\text{j},\text{v}}\approx {\mu }_{\text{j},\text{v}}^{* }{\left(\frac{T}{{T}^{* }}\right)}^{0.76}\end{eqnarray} \tag{ 42 }$$

The values at the 75 °C baseline temperature are given in Table S3. We modeled the mixture as ideal gas so they were independent of the cell pressure.

The density and viscosity of the gas were calculated from the gas composition and the diffusion of the molecules in the gas mixture was calculated using Maxwell-Stefan diffusion ⁶¹ in the Transport of Concentrated Species node.

$$\begin{eqnarray}&&{\rm{\nabla }}\cdot {\vec{j}}_{i}+\rho \left(\vec{u}\cdot {\rm{\nabla }}\right){\omega }_{i}={R}_{i}\end{eqnarray} \tag{ 43 }$$

The average velocity of the gas mixture is $\vec{u},$ and its density is $\rho .$ The mass fraction, mass source and mass flux of species i relative to the average velocity are ${\omega }_{i},$ R_i and ${\vec{j}}_{i},$ respectively. The mass flux, neglecting thermal diffusion, is ⁶⁰

$$\begin{eqnarray}&&{\vec{j}}_{i}=-\rho {\omega }_{i}\displaystyle \sum _{k=1}^{Q}{\widetilde{D}}_{{ik}}{\vec{d}}_{k}\end{eqnarray} \tag{ 44 }$$

Here ${\widetilde{D}}_{{ik}}$ are the multicomponent Fick's diffusion coefficients. Without external forces acting on the gas mixture and in the presence of total and partial pressure gradients, the diffusional driving force acting on species k, ${\vec{d}}_{\text{k}},$ (unit m⁻¹) is ⁶⁰

$$\begin{eqnarray}&&{\vec{d}}_{k}={\rm{\nabla }}{x}_{k}+\frac{1}{{p}_{A}}\left[\left({x}_{k}-{\omega }_{k}\right){\rm{\nabla }}{p}_{A}\right]\end{eqnarray} \tag{ 45 }$$

Here p_A is the total absolute pressure (see Eq. 38), x_k the mole fraction and ${\omega }_{k}$ the mass fraction of species k.

$$\begin{eqnarray}&&{x}_{k}=\frac{{\omega }_{k}}{{M}_{k}}{M}_{n}\end{eqnarray} \tag{ 46 }$$

Here M_k and M_n are the molar mass of species k and the effective molar mass of the mixture, respectively,

$$\begin{eqnarray}&&{M}_{n}={\left(\displaystyle {\sum }_{k}\displaystyle \frac{{\omega }_{k}}{{{\rm{M}}}_{{\rm{k}}}}\right)}^{-1}\end{eqnarray} \tag{ 47 }$$

We assumed ideal gas, so the density was calculated from the molar masses (M_k) and mass fractions of the species

$$\begin{eqnarray}&&\rho =\frac{{p}_{\text{A}}{M}_{\text{n}}}{{RT}}\end{eqnarray} \tag{ 48 }$$

The binary diffusion coefficients of the molecules were calculated using the Fuller-Schettler-Giddings model. ^60,62

$$\begin{eqnarray}&&{D}_{\text{AB}}=\frac{1.01325\cdot {10}^{-2}{T}^{1.75}\sqrt{\frac{1}{{M}_{\text{A}}}+\frac{1}{{M}_{\text{B}}}}}{p{\left[{\left({\sum }_{{\rm{A}}}{\rm{\nu }}\right)}^{1/3}+{\left({\sum }_{{\rm{B}}}{\rm{\nu }}\right)}^{1/3}\right]}^{2}}\end{eqnarray} \tag{ 49 }$$

Here the temperature is in Kelvins, pressure in Pa, molar masses M_A and M_B in g mol⁻¹ and the resulting diffusion coefficient is in cm² s⁻¹. Similarly to the other transport processes described in the previous sections, the Bruggeman model was applied to the gas diffusion.

$$\begin{eqnarray}&&{D}_{{AB},\text{eff}}={D}_{{AB}}{\varepsilon }_{\text{void}}^{1.5}\end{eqnarray} \tag{ 50 }$$

For the definition of the void volume fraction, see Eq. 18.

The values of the molar volumes $\left({\sum }_{{\rm{j}}}{\rm{\nu }}\right)$ used by Comsol to calculate the total molecular volume are given in Table S3. ^60,63 In general, the total volume of a molecule is based on the contribution of the contained groups summed together, but in our case the diffusion volumes of the molecules have specific values (that are different from the sum of the diffusion volume increments of the composite atoms).

At the PTL-gas channel boundary, we set a flux boundary condition for all species of the gas mixture that corresponded to the difference of the selected gas inflow composition and the true molar fractions.

$$\begin{eqnarray}&&{j}_{\text{k}}=Q\rho \frac{{M}_{\text{k}}}{{M}_{\text{n}}}\left({x}_{\text{k},\text{in}}-{x}_{\text{k}}\right)\end{eqnarray} \tag{ 51 }$$

The volumetric gas inflow rate per electrode area is Q (which has the dimension of velocity, e.g. ml min⁻¹ cm⁻² = cm min⁻¹), M_k is the molar mass of species k, M_n is the effective molar mass defined in Eq. 47 and x_k and x_k,in are the mole fractions of species k at the boundary and in the specified gas mixture flowing into the DAFC, respectively. Positive values correspond to more of species k flowing into the PTL and negative to the species flux being out of the electrode. The dry molar fractions f_k (that we used as a varied parameter in our study) and the total inflow mole fractions x_k,in differ from each other by a factor that depends on the RH and the absolute dry pressure.

$$\begin{eqnarray}&&{x}_{\text{k},\text{in}}=\frac{{f}_{\text{k}}}{1+\frac{\text{RH}{p}_{\text{w}}}{{p}_{\text{dry}}}}\end{eqnarray} \tag{ 52 }$$

We study the effect of the operating conditions and cell parameters on the DAFC performance by varying each parameter individually, while keeping the rest fixed to their baseline values. To estimate their relative importance, we change the parameters in each case by 20% of the baseline value (in case of temperature 20% of the value in degrees Celsius), which was selected so that all 20% variations would be within the physically realistic range. The parameter list as well as the values for the comparison to Zhao et al. ²³ and for the baseline and the 20% range of the sensitivity analysis are shown in Table I. The cell is assumed to be at uniform temperature, so different from the other parameters, the separate contributions of changes to the anode and cathode side were not evaluated. The gauge pressures are dry pressures and the total absolute pressure is the sum of the absolute dry pressure and the water vapor pressure, given by the RH and temperature (Eqs. 34 and 38). We simulated a wider value range for all parameters shown in the table, but this range was chosen to compare the importance of the parameters with each other, and in the most important cases, we discuss also the effects over a wider parameter value range. The catalyst loadings and gas flow rates were set equal at the anode and at the cathode for simplicity. The cathode flow rate 100 ml min⁻¹ cm⁻² was selected as the baseline simply because it was a round number, and the higher anode catalyst loading was selected to avoid reducing the cell performance too much by increasing the AOR kinetic overpotential. The dry fraction of O₂ at the cathode had to be reduced from 99% (not 100% to avoid numerical issues with 0% N₂ fraction) for the 20% increase to be physically feasible, so we chose to set the reduced fraction to be equal to the anode NH₃ dry fraction.

Table I. Operating conditions and cell parameters studied in the sensitivity analysis.

Parameter, (unit), [symbol]	Literature comparison	Baseline [and range] in the sensitivity analysis
Temperature (°C) [T]	60–95	75 [60–90]
Relative humidity [RH]	1	0.75 [0.6–0.9]
Anode AOR activity (A g−1) [A]	45–144	76.2 [61.0–91.5]
Cathode ORR activity at 0.9 V vs RHE (A g−1) [i0]	7.1–10	8.3 [6.6–9.9]
Anode catalyst loading (mg cm−2) [mcatalyst,A]	4	4 [3.2–4.8]
Cathode catalyst loading (mg cm−2) [mcatalyst,C]	3	4 [3.2–4.8]
Anode (dry) gauge pressure (bar) [pg,A]	0	1 [0.8–1.2]
Cathode (dry) gauge pressure (bar) [pg,C]	2	1 [0.8–1.2]
Ionomer volume fraction [${f}_{\varepsilon ,\text{m}}$]	0.85	0.75 [0.6–0.9]
Anode gas flow rate (ml min−1 cm−2)) [QA]	160	100 [80–120]
Cathode gas flow rate (ml min−1 cm−2)) [QC]	100	100 [80–120]
Anode NH3 molar fraction [${f}_{{\text{NH}}_{3}}$]	0.5	0.5 [0.4–0.6]
Cathode O2 molar fraction [${f}_{{\text{O}}_{2}}$]	0.99	0.5 [0.4–0.6]

The voltage loss over the AEM is calculated simply as the electrolyte potential difference between the electrode-AEM interfaces at the anode and the cathode sides. All other overpotentials are averages, based on the total power loss due to the process in question

$$\begin{eqnarray}&&{\eta }_{k}\left({i}_{A}\right)=\frac{\int {i}_{V}\left(x\right){\eta }_{k}\left(x\right){dx}}{\int {i}_{V}\left(x\right){dx}}\end{eqnarray} \tag{ 53 }$$

This is integrated over the thickness of the related electrode and the current source i_V has the dimension of A m⁻³. The geometric current density i_A (A m⁻² or more commonly mA cm⁻²) is simply the volumetric source integrated over the electrode thickness, which is equal to the electric current density out of the anode

$$\begin{eqnarray}&&{i}_{A}=\int {i}_{V}\left(x\right){dx}\end{eqnarray} \tag{ 54 }$$

The catalyst layers need to balance electron and ionic conduction with gas (and water) transport. Increasing the volume fraction of any conducting material (for gas transport, the empty space) increases the losses associated with the other two processes by reducing the associated conductance (following Eq. 15). From preliminary simulations, we already knew that the optimal ionomer volume fraction (${f}_{\varepsilon ,\text{ionomer}}$) in our model would be over 0.9 (corresponding to about 0.7 absolute volume fraction, Eq. 17). For the literature comparison, we chose to use slightly lower volume fraction to represent a very good, but perhaps not fully optimized catalyst layer. In the case of the sensitivity study, the baseline value had to be lower than 0.83 to enable a 20% increase. To avoid constricting the performance by limiting gas transport, and to match the temperature and RH, we chose the slightly lower 0.75 volume fraction (ca. 0.55–0.60 fraction of the total volume, depending on substrate, catalyst loading etc, Eq. 17). When varying the ionomer volume fraction and conductivity together so that the resulting layer resistance is constant, the ionomer volume has only a small effect on the DAFC polarization curve, except at high volumes (ca. 0.9 or more) where anode RH is first increased, increasing conductivity and improving performance (Fig. S5), and further increase would begin to limit reactant transport reducing the power.

The effect of the ionomer volume fraction on the DAFC power density at 0.25 V voltage is shown in Fig. 8. The maximum powers are ca. 30 mW cm⁻² for the CCM and ca. 24 mW cm⁻² for the CCS. The effects of the ionomer volume are very similar for both electrode types, and the power is more sensitive to the anode than to the cathode volume fraction. At very high volume fractions, the power is reduced because constricted gas transport decreases the Nernst voltage of the cell and increases both the kinetics and ion transport losses in the electrode, as illustrated in Figs. 9 and S4. Importantly for the cell and electrode optimization, the volume fraction of the highest power depends on the electrode structure. With the CCM, the optimal volume fractions are ca. 0.95 at both electrodes, whereas in the case of the CCS the peak is at around 0.9 for both electrodes. Because the transport at the level of microstructure is not included in our model, the total mass transport losses are likely underestimated (see e.g. ⁷⁰ ) and therefore the optimal ionomer volume fractions overestimated. Additionally, our relatively low current densities likely also bias the optimum towards high ionomer volume fractions, as higher current densities would require faster gas transport, i.e. lower ionomer volume fraction. Therefore, the optimal absolute numbers are probably not very important, but the main results are the differences in anode vs cathode and CCM vs CCS comparisons. In the case of CCS, the lower optimal volume fractions are likely due to the entire 0.5 mm electrode thickness becoming more constricted, whereas with CCM the affected CL thickness was 0.16 mm (cathode) or 0.34 mm (anode). The generally higher sensitivity to the anode volume fraction is likely due to the higher anode overpotentials and higher sensitivity to the ammonia concentration than to the oxygen concentration (discussed later).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The voltage losses as a function of the cathode and anode ionomer volume fractions at 75 mA cm⁻² current density are shown in Fig. 9. In this figure, and in the other similar figures, the solid lines and round markers correspond to the anode properties being varied and the dashed lines with open squares to the cathode properties being varied. Up to about 0.95 both the kinetic and ionic transport losses of the varied electrode decrease, but increase with further ionomer volume increase, especially at the anode, as the reduced reactant supply and product accumulation begin to slow the kinetics and even move the reaction farther from the AEM. The higher sensitivity of the anode side losses, and the increase at the highest volume fractions explains the pattern seen in Fig. 8. The thermodynamic voltage of the DAFC decreased with increasing volume fractions, except in the case of the cathode volume fraction of the CCM. The biggest effect came from the concentration polarization at the varied electrode, the CCM cathode being less sensitive than the other electrodes (CCM ca. 8.6 mV vs ca. 1 mV, CCS ca. 13 mV both), and most of the difference between the lowest and highest volume fractions came from the difference between the two highest simulated volume fractions (supporting information).

We assumed that each layer is a homogeneous mixture of the component materials and, as mentioned, neglected transport at smaller length scales. Considering the transport at different length scales would affect the transport in the ionomer-catalyst mixture, especially at high ionomer volume fractions, when considering that (more of) the catalysts would become deeper embedded in the ionomer, or the catalyst particles could become poorly connected with each other increasing electric resistance especially in CCM. This could add further restrictions that would begin to limit the DAFC operation at lower ionomer fractions but, at roughly equal ionic conductivities at anode and at cathode, and with the AOR maximum rate fixed to the rate A, the DAFC would likely remain more sensitive to the anode properties than to the cathode properties. While the optimal ionomer loadings could be lower than indicated here, in a CCS they could be somewhat lower than in a CCM due to the affected layer thicknesses. Additionally, the optimal ionomer loadings could in principle depend on the catalyst activities, but there are several qualifiers related to this statement. Higher activities enable higher current densities, which would require higher mass transport rates, so gas transport could then need more open electrode structures. Compared to HFCs, the optimal DAFC electrodes could therefore have slightly higher ionomer loadings and more closed structure, assuming that this does not negatively impact the electrocatalytic performance. This is of course temporary, until the DAFC current densities reach levels comparable to HFCs, neglecting the possible effect of different reactant gases, and it might not make sense to fully optimize a device for significantly lower performance levels than what is targeted.

The effect of catalyst loading on the power density at 0.25 V voltage is shown in Fig. 10. As shown earlier, at the baseline 4 mg cm⁻² loading the CCM achieves a higher power than the CCS, but the CCS is more sensitive to the catalyst loading and achieves significantly higher power densities at high catalyst loadings. The main reason for this is the anode ion transport overpotential. While the kinetic losses were reduced by increased catalyst loading, the increase in ion transport losses reduced the total gains, similarly to the results of Singh et al. about AEM-HFC cathode. ⁶⁴ Although Hu et al. ²⁸ found experimentally the optimal DAFC cathode catalyst loading to be only 2 mg cm⁻², the relatively low effect of the catalyst loading compared to other varied parameters is similar to our results. Similar behavior has been observed also in PEM-HFC simulations. ⁷¹ The effect of the catalyst loadings on the voltage losses at 75 mA cm⁻² current density is shown in Fig. 11. The CCS anode differs from the other three electrodes in that increasing the catalyst loading decreases both the kinetic and the ionic transport losses in the electrode, explaining its higher sensitivity. The behavior of the other three electrodes was more similar to each other, albeit the changes in the overpotentials were different in their magnitude: with increased catalyst loading the kinetic overpotential was reduced and the transport losses barely changed, and the kinetic overpotential of the other electrode (constant 4 mg cm⁻² loading) remained constant. As an exception, in the case of CCM, increasing catalyst loading increased the ionomer losses at catalyst loadings below 3 mg cm⁻² (Fig. 11a).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

As discussed in the model description, with CCM we assumed that increasing the catalyst loading increases the catalyst layer thickness linearly (constant catalyst volume fraction), whereas with the CCS increased loading corresponds to more catalyst in the same thickness. Therefore, especially with the maximum rate limitation of the AOR kinetics (Eq. 5), increasing the catalyst loading per electrode thickness increases the rate near the AEM, bringing the reaction closer to it, thus reducing the ionic transport distance and losses. In the case of the CCS cathode, i.e., Butler-Volmer kinetics, the ionic transport losses in the electrode are not appreciably affected by the catalyst loading. Therefore, we ascribe the sensitivity of the CCS anode to the mathematical form of the AOR kinetics, i.e., the RDS not being electrochemical by nature.

With both electrode types, varying the catalyst activity produces results that are very similar to the effect of the catalyst load on the CCS: At constant 75 mA cm⁻² current density, the kinetic overpotential of the varied catalyst responds expectedly, with almost all other losses remaining constant, as shown in Fig. S6. However, the anode ionomer losses decrease with increasing anode activity, due to the mathematical form of the AOR kinetics and the eventual saturation to A. Compared to varying the catalyst loading and to the baseline scenario, the relative range is wider, but nevertheless the CCM anode shows clear differences in the response of the ionic transport losses: ionic transport losses increase with increased catalyst loading (i.e. CL thickness), but reduce with increased catalyst activity. In contrast, for the CCS, the responses to catalyst loading and activity are similar because the catalyst loading only affects the catalyst volume fraction which is small (and thus indirectly the absolute ionomer volume) but does not affect the layer thickness.

Increasing the anode loading with CCM expectedly decreased the AOR kinetic losses and, interestingly, increasing the cathode loading also reduced the AOR kinetic losses, which was clearest at over 10 mg cm⁻² anode loadings (Fig. S7). Although the effect was smaller, the same happened to the ORR kinetics at high cathode catalyst loadings when increasing the anode loading. In both cases the difference was visible only at over ca. 50 mA cm⁻². With CCS this effect was practically absent (Fig. S8). The practical effect is quite small, especially at lower catalyst loadings, but this demonstrates that the electrode kinetics and transport losses are not fully independent of each other and for the best accuracy this should be considered. In principle, cell optimization therefore should not be done one electrode at a time, independent of the other, but the effect on practical device optimization is probably negligible.

The enhancement of the unchanged electrode was due to water transport. Increasing the CCM catalyst loading increases the electrode thickness, which then increases the RH and ionomer (and AEM) conductivity in the cell, reducing ionic transport losses. This in turn slightly increases the reaction rate throughout the CL at constant overpotential. The same current density can therefore be achieved with lower total voltage losses when the catalyst loading of the other electrode is increased. Increasing the cathode loading increased the anode RH and reduced the cathode RH (supporting information, Fig. S9), and vice versa with increasing anode loading, due to the changed water transport properties. Increasing electrode thickness made water removal through it more difficult so the RH was increased throughout the cell (maximum RH is shown in Fig. S9) and more water was removed through the electrode whose thickness was not changed, as diffusive transport through it became comparatively easier. The maximum RH is increased, and it might be moved slightly towards the cathode (increased cathode loading), or deeper into anode and farther from the AEM (increased anode loading), but the operation of both electrodes is affected by the changes to one. Although, as will be discussed in the next section, it seems likely that the flooding would practically always begin in the anode, in principle and with a very unusual combination of electrode thicknesses (and other properties), flooding could perhaps be forced to begin in the cathode.

With both electrode structures, the RH of the gas feeds had the largest impact on the cell power. However, with the CCM the cathode side humidity affected the power more than the anode side, whereas the opposite was true with the CCS. The effect of the anode side relative humidity on the power was approximately equal in both cases. The voltage losses at 50 mA cm⁻² current density as a function of the RH are shown in Fig. 12. Figure (a) shows the losses for the CCM and (b) for the CCS. The current density of the comparison was reduced from the 75 mA cm⁻² because reducing the anode RH with CCS limited the maximum current density to below this value. In all cases, increasing RH decreased both kinetic overpotentials. In fact, the only loss component that does not clearly decrease with increasing RH is the cathode ionic transport losses in CCS. Due to the low current density this loss is quite low, so the ca. 10% reduction amounts to only few millivolts. The clearest differences between the electrode types are the sensitivity of the anode ionomer and AEM voltage losses, and the AOR kinetics on the anode and cathode RH. With CCM, both the anode ionomer and the AEM losses are more sensitive to the cathode RH than to the anode RH, whereas with the CCS the anode ionomer losses are more sensitive to the anode RH and there is no clear difference in how the AEM losses respond to the RH. Additionally, the AOR kinetics in CCS are more sensitive to the anode RH than to the cathode RH, with practically no difference between the two in CCM. Therefore, the greater sensitivity of the CCM to the cathode RH is a combination of the increased sensitivity of the following losses to the cathode RH (or reduced sensitivity to the anode RH)

• 1.  
  AOR kinetics
• 2.  
  anode ionomer losses
• 3.  
  AEM losses

Zoom In Zoom Out Reset image size

At least a part of the increased sensitivity to the cathode side RH is that the cathode CL of the CCM is thinner than the other simulated ionomer-containing layers, making it the least resistive to water transport through it. The water flux at both cathodes is in from the O₂ flow, but with CCS the magnitude is higher, so from this it could be expected that CCS would be more strongly coupled to cathode RH through the cathode-side outflux boundary condition (Eq. 51). As the opposite is true for the sensitivity of the power density to the cathode RH, and because the effects on the anode side losses were an important difference, the CL layer thickness and water transport resistance appears a more likely explanation. The reduced AOR losses are probably due to the increased ionomer conductivity increasing the driving kinetic overpotential farther from the AEM. This increases the AOR rate there, meaning that, for the same total current density, the AOR rate near the AEM could be lower, reducing the kinetic overpotential. Due to the functional form of the AOR kinetics, at the ca. 400 mV overpotentials where the current density is near the saturation to A (see Fig. 1, a relatively small reduction in the local reaction rate could reduce the overpotential more than the corresponding increase in the AOR rate elsewhere requires the overpotential to increase.

The spatial distribution of water activity (RH) for the baseline scenario at 93 mA cm⁻² current density is illustrated in Fig. 13. Qualitatively, the CCM profile is in agreement with the neutron radiography results of Omasta et al. ⁷² and the simulated water content profile of Gerhardt et al. ⁵² (in their SI) for AEM HFC. Despite the different FC types (HFC vs DAFC), a direct comparison is valid because the most important physics and boundary conditions of water transport are the same (humidified gaseous feeds, AEM, the same numbers of electrons per water molecule in both reactions). The electro-osmotic drag, generation of water by the AOR, and consumption by the ORR bias the peak RH towards the anode, and back diffusion to the cathode is responsible for the RH remaining slightly higher than the inflow RH near the AEM. While in our case the peak is in the CL side of the anode CL-AEM interface (Omasta et al. had it in the AEM), this is likely influenced by our higher catalyst loadings, thicker CLs, and thinner AEM. Because the water outflow rate depends directly on the difference to the inflow RH, in this case 0.75 (see Eq. 51), the greatest water outflux is at the CCS anode and at the cathodes the flux is in from the gas channel rather than out of the electrodes. The maximum RH in the CCS is higher than in CCM, which could lead to earlier flooding in the cell, and especially CCS cathode might suffer from drying at low cathode inflow RHs. The cathode water influx is larger with CCS so the magnitude of the flux likely does not explain the higher sensitivity of CCM to the cathode RH, which might be due to the CCS electrode being more restrictive to gas transport than the substrate without the ionomer in the case of CCM.

Zoom In Zoom Out Reset image size

As a difference to experimental results that indicate CCS to have lower mass transport resistance than CCM, e.g., ^41,73 in our case the CCS has no clear separate CL, and catalyst and ionomer fill the entire PTL substrate thickness, and this difference would affect mass transport through the electrode. Our catalyst loading is also significantly higher than in HFC literature, about 10-fold difference (or even more), which could affect the comparison. It is not clear how a reduced CL thickness and partially filled porous substrate would affect water transport, but we consider this analysis to be outside the scope of the study. The simulation model could also overestimate the diffusion coefficients inside the porous layers, which could affect the overall mass transport in the CCM and CCS simulations differently due to different layer thicknesses and porosities.

The gas flow rate has an important role in humidity management by controlling the removal rate of the produced water. Compared to PEM HFCs with dry cathode stream, the humidified cathode stream and electro-osmotic drag moving water to the anode reduces the capacity of the cathode stream to accomplish this, similarly to AEM HFCs. In fact, although the simulated AEM was quite thin, in principle helping back diffusion towards the cathode, varying the cathode gas flow rate had practically no effect on the DAFC output power or the RHs at the electrode - gas channel interfaces. The power density as a function of gas flow at 0.25 V and at 50 mA cm⁻² current density is shown in Fig. 14. Except for the sharp decrease below ca. 0.5 ml min⁻¹ cm⁻², the cathode side flow has only a very small effect on the DAFC output power and both the anode and cathode RHs (Fig. S9). With the anode flow, there is a maximum at about 3 ml min⁻¹ cm⁻² due to reduced water removal, i.e. increased RH and ionomer conductivity (Fig. S9), but not yet too diminished ammonia supply. In practice, the larger effect of the changed downstream composition (Fig. 7) would mean that reducing the gas flow rate almost certainly would reduce the DAFC power, instead of increasing it.

Zoom In Zoom Out Reset image size

The peak in power at about 3 ml min⁻¹ cm⁻² is because our model is a point model, when considering the geometric area of the electrodes and thus does not account for any nonuniformities, due to e.g. gas composition change along the flow channel, over the electrode area. Additionally, the downstream effects could make the DAFC more sensitive to the gas composition or flow rate than our simulations suggest. That said, optimization for better humidity retention in the cell could improve the DAFC performance, provided that the transport of the other species would not be severely impeded.

The gas flow rate can be interpreted as either the gas volume flow per electrode area or the gas velocity (as mentioned earlier, ml min⁻¹ cm⁻² is in fact equal to cm min⁻¹). Because we do not consider nonuniformities over electrode area or channel length, the simulated effects probably match small electrode areas the best. With increasing electrode area, the gas velocity typically increases faster than the volume flow per electrode area: The average velocity is the function of the flow field channel layout and the total gas flow rate. Because the gas channel cross section typically does not increase as fast as the electrode area, increasing the electrode area while maintaining a constant gas flow per electrode area would increase the gas velocity in the channels. Therefore, the effects primarily related to the reactant supply rate (flow rate per area) and gas velocity may impact a DAFC differently with different electrode areas.

Figure 15 shows the effect of the fraction of ammonia and oxygen of the dry anode and cathode inflows, respectively (f_k in Eq. 52). In the baseline scenario the fraction of both is 0.5, or 50%. Similarly to most scenarios with the other parameters, in general the CCM produces few mW cm⁻² higher power at all values. The pattern of the Nernst voltage (Fig. S12) at the same cell voltage shows the same pattern with the highest power matching the highest thermodynamic voltage. The reaction kinetics are also affected, with increasing ammonia fraction decreasing the AOR and oxygen fraction ORR losses (Fig. S13), except that AOR losses increase slightly when ammonia fraction increases to over 0.9, and the further increased power is probably due to increased thermodynamic voltage. In both cases the ammonia concentration is more important than the oxygen concentration, higher power densities can be achieved or maintained with low oxygen concentration than with low ammonia concentration. The reaction stoichiometry certainly affects this: three electrons per NH₃ molecule and four per O₂, meaning that at any current density NH₃ is consumed faster than O₂, and the AOR Nernst potential is more sensitive to the NH₃ partial pressure than the ORR potential to the O₂ partial pressure (∼20 mV/decade vs ∼15 mV/decade). The difference may be smaller than in HFCs (2 electrons per H₂, thus ∼30 mV/decade), but the NH₃ supply rate likely still needs to be higher than the O₂ supply rate to avoid gas composition problems along the flow channels.

Zoom In Zoom Out Reset image size

The AOR and ORR concentration overpotentials (Eqs. 12 and 13) for the CCM are shown in Fig. 16. In the case of the AOR, negative concentration polarization corresponds to increased DAFC voltage, whereas in the case of the ORR more positive potentials increase the DAFC voltage. The differences compared to the CCS were negligible. While the contours are neither perfectly horizontal nor vertical, they are very close to it, demonstrating that other factors did not meaningfully contribute to the equilibrium potentials. Comparing the AOR and ORR potentials to each other, the difference between the maximum and minimum for the ORR was about 30 mV compared to about 65 mV of the AOR, quantifying the higher sensitivity to the ammonia concentration. This result was of course affected by our choice of using N₂ as the non-reactant component of the dry gas mixture. Alternatively, N₂ being a very common inert gas means that its contribution to the AOR Nernst potential magnifies the sensitivity of DAFCs to the anode gas composition, when compared to other FC types. Nevertheless, changes to the ammonia pressure corresponded to about 45 mV and N₂ pressure to about 20 mV of the total 65 mV change (Fig. S14), meaning that NH₃ alone affected the thermodynamic voltage about 50% more than O₂. In practice, it may therefore be necessary to operate a DAFC with higher anode side flow rate or/and pressure than the cathode side to minimize voltage losses especially due to composition gradients along the flow channel.

Zoom In Zoom Out Reset image size

While the performance in Fig. 15 correlates positively with the thermodynamic voltage of the cell, in other cases this connection does not necessarily exist. For example, in the case of the catalyst loading, the highest powers corresponded to the lowest thermodynamic voltages (Figs. 10 and S15), and in the case of CCM to the highest peak RH inside the cell (Fig. S10), the match between the peak RH and Nernst voltage being especially striking. Therefore, the best performance corresponded to the best capability to extract power from the available reactants, not to the best gas transport. In the case of the catalyst loading and CCM specifically, the general trend of the thermodynamic voltage appears to be dictated by the anode and the peak RH, as increased cathode RH would increase the voltage. Increased anode RH would, of course, reduce ion transport losses and maybe also reduce the average kinetic overpotential. Similarly to Fig. 16, the anode RH is more important to the Nernst voltage than the cathode RH due to reaction stoichiometries, AOR (reaction (2)) producing a water molecule for each electrode, whereas ORR (reaction (3)) requires 2 electrons for each consumed water molecule, further amplifying the importance of the anode gas composition on the DAFC performance compared to the cathode.

Compared to the temperature and the RH, the effect of increased pressure was small. However, when considering the ultimate, achievable performance of a DAFC, the pressure can be multiplied, whereas the maximum RH and temperature increases are limited closer to the 20% of our sensitivity analysis. (The temperature limit is due to material constraints, and solid oxide fuel cells could be operated at much higher temperatures. ¹¹ ) As shown in Table II, the anode gauge pressure had a larger impact on the power density than the cathode pressure and increasing it to 5 bar (with constant mass flow) increased the power to 29.1 mW cm⁻² with CCM and to 23.7 mW cm⁻² with CCS. Pressurizing both electrodes to 5 bar increased the power densities to 32.8 mW cm⁻² and 26.1 mW cm⁻² with CCM and CCS, respectively. Compared to the 20% changes to the other parameters in Fig. 7, increasing the anode pressure has roughly the same impact as increasing the RH of one electrode, and pressurizing both electrodes to 5 bar has roughly the same impact as increasing the cell temperature. The effect of the gauge pressures is plotted in Fig. 17. Therefore, although minor changes to pressure had only a small impact on the cell performance, if the cell can be properly pressurized, the performance gains are comparable to the effects of increased temperature and RH. The main cause of the enhanced performance were the reduced kinetic overpotentials, especially of AOR with increased anode pressure, as shown in Fig. S16. The other voltage losses were almost not affected at all.

Zoom In Zoom Out Reset image size