Three-Dimensional Modeling of Anion Exchange Membrane Electrolysis: A Two-Phase Flow Approach

Abstract

Anion exchange membrane water electrolyzers (AEMWEs) are attracting growing interest as a green hydrogen production technology. Unlike proton exchange membrane (PEM) systems, AEMWEs operate in an alkaline environment, allowing one to use less expensive, non-noble materials as catalysts for the reactions and non-fluorinated anion exchange polymer membranes. However, the performance and stability of AEMWEs strongly depend on the alkaline electrolyte concentration. In this work, a three-dimensional multi-physics model considering two-phase flow effects is applied to understand the impact of KOH electrolyte concentration and its flow rate on AEMWE performance, as well as on the current and gas volume fraction distributions. The numerical results were compared to experimental data published in the literature. For current densities above 1 A/cm², a strongly non-uniform H₂ and O₂ gas volume distribution could be evidenced by the 3D simulations. Increasing the KOH electrolyte flow rate from 10 to 100 mL/min noticeably improves cell performance for current densities above 1 A/cm². These results show the importance of accounting for the three-dimensional geometry of an AEMWE and two-phase flow effects to accurately describe its operation and performance.

1. Introduction

As a high-potential energy carrier, hydrogen can be used for both storing and regulating electricity generated from renewable sources. Hydrogen has also attracted attention since it can replace fossil energies for various industrial applications that significantly contribute to greenhouse gas (GHG) emissions, such as the steel industry and fertilizer production. Green hydrogen produced through water electrolysis from renewable electricity is a promising solution to decarbonizing these sectors. Among the low-temperature water electrolysis technologies, membrane electrolyzers have received substantial attention, allowing for high purity hydrogen production at high current densities, compact cell design, and dynamic operation at high differential pressure [1,2].

There are two types of membrane electrolyzers: proton exchange membrane water electrolyzers (PEMWEs) and anion exchange membrane water electrolyzers (AEMWEs). PEM technology has benefited from intensive research for several years. Its operation in an acidic environment allows high current densities to be reached (>2 A/cm² for applied voltages below 2 V) [3]. In addition, PEM technology is well adapted to load variations, which is particularly suitable for renewable energy storage. PEMs are predominantly made from a polyfluoroalkyl polymer (PFAS), which exhibits good mechanical and thermal stability due to its very stable carbon-fluorine bonds. However, PFAS is problematic as it is highly toxic and polluting with low degradation and high soil/water contamination [4]. Thus, growing interest has been focused on AEMWE technology operating in an alkaline environment, which allows for the use of non-precious and/or less expensive metals as catalysts, diffusion layers, and bipolar plates. Among these, nickel (Ni) and its iron (Fe)-, cobalt (Co)-, or aluminum (Al)-based alloys are predominantly used. However, these catalysts do not achieve the desired performance [5,6]. Therefore, it is necessary to use thicker catalytic layers and a supply with an alkaline electrolyte, usually KOH or NaOH, to achieve similar current densities [7]. Furthermore, unlike PEM systems, there is a wide variety of non-fluorinated anion exchange polymer membranes (AEMs) with lower environmental impact than the ones used in acidic technologies [8]. Despite recent progress, there is a development lag of AEMWE systems compared to others [1]. It is necessary to notably increase their lifespan and improve their ability to operate at high currents. A significant increase in electrode overpotentials is observed in high-current domains, lowering energy efficiency and potentially compromising the durability of electrolyzers. It is indeed under these high-current conditions that we observe significant gas evolution of dihydrogen and dioxygen, which can lead to a limitation of electrochemical reactions by mass transport [2]. It is therefore essential to study two-phase flows in order to better understand the bubble formation process and the transport of H₂ and O₂ bubbles and their impact on reaction mechanisms. For example, it has been shown that bubbles dispersed in the electrolyte of AEMWEs hinder the proper migration of hydroxide ions (OH⁻), thereby reducing the effective conductivity of the electrolyte [9]. In addition to preventing ions from reaching the reaction surface, bubbles attached to it lead to a non-uniform distribution of current density, affecting the charge transfer process in these regions [10,11]. Effective bubble removal is thus linked to the wettability and hydrophilicity properties of the materials used. Electrochemical kinetics at interfaces are thus directly influenced by flows in the layer, which can cause heterogeneities and strong local gradients if mismanagement of these flows is observed. The formation, release, and transport of bubbles are complex phenomena. They depend not only on gas interactions with the liquid phase but also with the solid matrix (support/catalyst/ionomer), which can locally modify the microstructure and porosity of the layer, highlighting the importance of the nature of the materials used and their properties. To date, there is a lack of understanding of these phenomena [12], especially in AEM electrolyzers [13]. The strong attractiveness of this technology (which has been growing rapidly in recent years) necessitates improving our understanding of its operation. Moreover, it remains challenging to experimentally observe these phenomena [14,15,16]. Characterizing two-phase flows is therefore fundamental for studying reaction processes; the strong couplings between charge transfer kinetics at interfaces and mass transport make the study complex [17].

To date, a lot of physics-based models have been developed for alkaline water electrolyzers [18] to analyze their performance, but only a few are dedicated to AEMWEs. An et al. [19] and Lui et al. [20] developed a through-plane 1D model that couples mass transport, charge transport, and electrochemical reactions in an AEMWE. Very recently, Lawand et al. [21] developed a 3D model of a membrane electrode assembly (MEA) composed of Ni Raney at the cathode and NiFe₂O₄ at the anode. The influence of several parameters, such as the electrolyte pH, the membrane, and electrode thicknesses, was investigated.

In this work, a modeling approach of a full AEMWE cell is proposed to evaluate the impact of gas production at both electrodes. Our 3D model is based on the previous 1D model from Liu’s work [20] and coupled multiphysics transports with electrochemical reactions. It provides a more accurate description of the actual geometry of the electrolyzer, the electrodes, flow channels, and membrane, and an improved description of the mass transport processes (flow pattern and concentration gradient). Nonuniformities in the current density, O₂ and H₂ production, and concentration can also be considered in 3D simulations. A two-phase flow approach allows one to investigate the impact of the transport of gas in the channel and the porous media, namely, the catalyst layer where hydrogen and oxygen gas evolution take place. Finally, the influence of KOH concentration and the electrolyte flow rate is investigated.

2. Model Description

2.1. Mathematical Modeling

A 3D, continuum, two-phase, isothermal model was developed to study the impact of KOH concentration and flow rate on AEM electrolysis performance. In the following sections, we introduce all of the physics implemented in the model. The AEM electrolyzer model is developed by considering the following assumptions:

• The fluid mixture is assumed to have Newtonian properties;
• A homogeneous flow is considered in the channels, which means that the gas and liquid velocities are equal;
• No-slip condition is assumed at the free-porous interface on the free-flow side, which means that the tangential component of the mixture velocity is set to zero;
• The oxygen (O₂) and hydrogen (H₂) produced are considered to be present as pure gas bubbles dispersed within the continuous liquid electrolyte;
• The gas and KOH aqueous electrolyte flows are assumed to be laminar due to the low flow velocity and small pressure gradient;
• Hydrogen and oxygen crossover through the membrane are neglected;
• The porous media are isotropic, and the transport phenomena occurring within them are characterized using effective properties;
• The electrochemically active surface area is assumed to be in contact with the ionomer phase and therefore to be constant whatever the KOH concentration.
• The ionic conductivity of the ionomer is assumed to depend on the KOH concentration.
• An isothermal model is assumed.

2.2. Multiphase Flow

The multiphase flow is studied for both liquid and gas in the channels and porous media. The mass transport of the phase k is described for each phase k ∈ l, g.

$$\nabla ·\left({\mathsf{\rho }}_{\mathrm{k}}{\mathrm{u}}_{\mathrm{k}}\right)={\mathrm{R}}_{\mathrm{k}}$$

(1)

where ${\mathsf{\rho }}_{\mathrm{k}}$ and ${\mathrm{u}}_{\mathrm{k}}$ are the density and velocity of phase k, respectively. Here, ${\mathrm{R}}_{\mathrm{k}}$ represents the mass source term of phase k, which is zero in all of the layers except in the catalyst layers (

Table 1. Source terms for mass conservation (Equation (1)).

Component	Gas	Liquid
ACL 1 (positive electrode)	$+\frac{{\mathrm{i}}_{\mathrm{v}}}{4\mathrm{F}}{\mathrm{M}}_{{\mathrm{O}}_2}$	$+{\mathrm{R}}_{\mathrm{ml}} {\mathrm{M}}_{{\mathrm{H}}_2\mathrm{O}}$
CCL 2 (negative electrode)	$-\frac{{\mathrm{i}}_{\mathrm{v}}}{2\mathrm{F}}{\mathrm{M}}_{{\mathrm{H}}_2}$	$+{\mathrm{R}}_{\mathrm{ml}} {\mathrm{M}}_{{\mathrm{H}}_2\mathrm{O}}$

¹ ACL: anode catalyst layer, ² CCL: cathode catalyst layer.

).

The liquid electrolyte is assumed to drive the gas in the channels. Consequently, the gas and liquid velocities are considered to be equal by neglecting any drift velocity and considering that both phases share the same pressure. A mixture model is then implemented for coupling the mass phase transport (Equation (1)) and the total mass (Equation (2)) and momentum (Equation (3)) transport of the mixture as a single-phase continuum flow with macroscopic properties.

$$\nabla ·\left({\mathsf{\rho }}_{\mathrm{m}}{\mathrm{u}}_{\mathrm{m}}\right)=0$$

(2)

$${\mathsf{\rho }}_{\mathrm{m}}\left({\mathrm{u}}_{\mathrm{m}}·\nabla \right){\mathrm{u}}_{\mathrm{m}}=\nabla ·\left[-\mathrm{pI}+{\mathsf{\mu }}_{\mathrm{m}}\left(\nabla {\mathrm{u}}_{\mathrm{m}}+{\left(\nabla {\mathrm{u}}_{\mathrm{m}}\right)}^{\mathrm{T}}\right)-\frac{2}{3}{\mathsf{\mu }}_{\mathrm{m}}\left(\nabla ·{\mathrm{u}}_{\mathrm{m}}\right)\mathrm{I}\right]$$

(3)

where p and I are the total pressure and identity matrix, respectively. Here, ${\mathsf{\rho }}_{\mathrm{m}}$ and ${\mathsf{\mu }}_{\mathrm{m}}$ represent the mixture density and dynamic viscosity determined by volume-averaging (Equations (4) and (5)).

$${\mathsf{\rho }}_{\mathrm{m}}={\displaystyle \sum }_{\mathrm{k}}{\mathrm{s}}_{\mathrm{k}}{\mathsf{\rho }}_{\mathrm{k}}$$

(4)

$${\mathsf{\mu }}_{\mathrm{m}}={\displaystyle \sum }_{\mathrm{k}}{\mathrm{s}}_{\mathrm{k}}{\mathsf{\mu }}_{\mathrm{k}}$$

(5)

where ${\mathrm{s}}_{\mathrm{k}}$ is the volume fraction of phase k.

The flow velocity ${\mathrm{u}}_{\mathrm{k}}$ of each phase k is performed by the extended Darcy’s law in the porous media.

$${\mathrm{u}}_{\mathrm{k}}=-\frac{{\mathrm{k}}_{\mathrm{r},\mathrm{k}}}{{\mathsf{\mu }}_{\mathrm{k}}}{\mathrm{k}}_{\mathrm{sat}} \nabla {\mathrm{p}}_{\mathrm{k}}$$

(6)

where ${\mathrm{k}}_{\mathrm{r},\mathrm{k}}$ and ${\mathrm{k}}_{\mathrm{sat}}$ are the relative and absolute (or saturated) permeabilities. The relative permeability is given by the third power of the volume fraction of each phase.

$${\mathrm{k}}_{\mathrm{r},\mathrm{k}}={\mathrm{s}}_{\mathrm{k}}^3$$

(7)

The absolute permeability is based on the Kozeny–Carman equation for porous media.

$${\mathrm{k}}_{\mathrm{sat}}=\frac{1}{16} {\mathrm{d}}_{\mathrm{CK}}^2\frac{{\mathsf{\epsilon }}^3}{{\mathrm{k}}_{\mathrm{CK}}{\left(1-\mathsf{\epsilon }\right)}^2}$$

(8)

where ${\mathrm{d}}_{\mathrm{CK}}$, ${\mathrm{k}}_{\mathrm{CK}}$, and $\mathsf{\epsilon }$ represent the characteristic diameter, an empirical geometric constant, and the porosity, respectively (

Table 2. Components, phases present, and material properties taken from [20].

			AEM	ACL	CCL	APTL 2	CPTL 3	CHANNEL
Phases present 1	-	-	I	i, g, l, e	i, g, l, e	g, l, e	g, l, e	g, l
Porosity	$\mathsf{\epsilon }$	-	0	0.4	0.4	0.4	0.8	1
Ionomer volume fraction	${\mathsf{\epsilon }}_{\mathrm{i}}$	-	1	0.15	0.15	0	0	0
Electronically conductive volume fraction	${\mathsf{\epsilon }}_{\mathrm{e}}$	-	0	0.45	0.45	0.6	0.2	0
Characteristic diameter	${\mathrm{d}}_{\mathrm{CK}}$	μm	0	0.25	0.25	7.6	7.6	0
Empirical constant	${\mathrm{k}}_{\mathrm{CK}}$	-	0	9.375	9.375	4.06	4.06	0

¹ i: ionic, e: electronic, g: gas, l: liquid. ² APTL: anode porous transport layer, ³ CPTL: cathode porous transport layer.

).

The pressure difference between liquid and gas in the porous media was determined by computing the capillary pressure, ${\mathrm{p}}_{\mathrm{c}}$, according to the Leverett J-function, $\mathrm{J}\left({\mathrm{s}}_{\mathrm{l}}\right)$.

$${\mathrm{p}}_{\mathrm{c}}=\mathsf{\sigma }\cos \left(\mathsf{\theta }\right) {\left(\frac{\mathsf{\epsilon }}{{\mathrm{k}}_{\mathrm{sat}}}\right)}^{\mathrm{½}}\mathrm{J}\left({\mathrm{s}}_{\mathrm{l}}\right)$$

(9)

$$\mathrm{J}\left({\mathrm{s}}_{\mathrm{l}}\right)=\{\begin{array}{cc}1.417 \left(1-{\mathrm{s}}_{\mathrm{l}}\right)-2.120 {\left(1-{\mathrm{s}}_{\mathrm{l}}\right)}^2+1.263 {\left(1-{\mathrm{s}}_{\mathrm{l}}\right)}^3\hfill & \mathrm{for} \mathsf{\theta }<90°\hfill \\ 1.417 {\mathrm{s}}_{\mathrm{l}}-2.120 {{\mathrm{s}}_{\mathrm{l}}}^2+1.263 {{\mathrm{s}}_{\mathrm{l}}}^3\hfill & \mathrm{for} \mathsf{\theta }>90°\hfill \end{array}$$

(10)

where $\mathsf{\sigma }$ and $\mathsf{\theta }$ are the surface tension and contact angle, respectively (

Table 3. Material and electrochemical properties used in the model taken from [20,22].

Symbol	Parameter	Values	Unit
$\mathsf{\sigma }$	surface tension	66.24	mN/m
$\mathsf{\theta }$	contact angle	50	°
${\mathsf{\sigma }}_{\mathrm{e}}^{\mathrm{bulk}}$	bulk electronic conductivity [20]	120	S/cm
$\mathsf{\xi }$	electroosmotic coefficient [20]	1	-
$\mathsf{\alpha }$	charge transfer coefficient for reduction	OER: 0.5 HER: 0.5	-
$\mathsf{\beta }$	charge transfer coefficient for oxidation	OER: 0.5 HER: 0.5	-
$\mathrm{IEC}$	ion exchange capacity	2.353	mmol/g
${\overline{\mathrm{V}}}_{\mathrm{l}}$	molar volume of liquid water	18.307	mL/mol
${\mathsf{\mu }}_{\mathrm{l}}$	viscosity of liquid water [20]	$0.00175\exp \left(-\frac{16,000\frac{\mathrm{J}}{\mathrm{mol}}}{\mathrm{R}}\left(\frac{1}{{\mathrm{T}}_{\mathrm{t}}}-\frac{1}{\mathrm{T}}\right)\right)$	Pa·s
${\mathsf{\mu }}_{\mathrm{g}}$	viscosity of gas	${\mathrm{O}}_2:$ $2.026\times {10}^{-5}$ ${\mathrm{H}}_2:$ $8.789\times {10}^{-6}$	Pa·s
${\mathsf{\rho }}_{\mathrm{l}}$	density of liquid water	983.2	kg/m3
${\mathsf{\rho }}_{\mathrm{g}}$	density of gas	$\frac{{\mathrm{p}}_{\mathrm{g}}}{\mathrm{RT}}{\mathrm{M}}_{\mathrm{i}}$ i: ${\mathrm{O}}_2$ (anode), ${\mathrm{H}}_2$ (cathode)	kg/m3
${\mathrm{c}}_{{\mathrm{OH}}^{-},\mathrm{ref}}$	reference hydroxide concentration	0.1	M
${\mathrm{z}}^{-}$	charge of hydroxide ions	−1	-
${\mathrm{p}}_{\mathrm{ref}}$	reference pressure	1	atm
${\mathrm{p}}_{\mathrm{t}}$	triple point pressure [22]	611.2	Pa
T	operating temperature	60	°C
${\mathrm{T}}_{\mathrm{t}}$	triple point temperature [22]	273.16	K
${\mathrm{H}}_{\mathrm{l}}$	enthalpy of liquid water [22]	0	J/g
${{\mathrm{C}}_{\mathrm{p}}}_{\mathrm{l}}$	heat capacity of liquid water [22]	4.217	J/(g·K)
${\mathrm{M}}_{{\mathrm{H}}_2\mathrm{O}}$	molar mass of water	18	g/mol
${\mathrm{M}}_{{\mathrm{O}}_2}$	molar mass of dioxygen	32	g/mol
${\mathrm{M}}_{{\mathrm{H}}_2}$	molar mass of dihydrogen	2	g/mol
${\mathrm{k}}_{\mathrm{ml}}$	constant for absorption and desorption of liquid water [22]	1	g·mol/(cm3·s·J)

).

2.3. Charge Transport

Electron transport in the electronically conducting phase as the ion transport in the ionomer phase (membrane and catalyst layer) was studied using Ohm’s law. The resulting conservation equations for the charge transport are presented in Equations (11)–(13).

$$\begin{gathered} \nabla ·{\mathrm{i}}_{\mathrm{i}}=\nabla ·\left(-{\mathsf{\sigma }}_{\mathrm{i}}^{\mathrm{eff}}\nabla {\mathsf{\Phi }}_{\mathrm{i}}-\frac{{\mathsf{\sigma }}_{\mathrm{i}}^{\mathrm{eff}}\mathsf{\xi }}{{\mathrm{z}}^{-}\mathrm{F}}\nabla {\mathsf{\mu }}_{{\mathrm{H}}_2\mathrm{O},\mathrm{i}}\right)={\mathrm{i}}_{\mathrm{v}} \mathrm{CLs} \\ \nabla ·{\mathrm{i}}_{\mathrm{e}}=\nabla ·\left(-{\mathsf{\sigma }}_{\mathrm{e}}^{\mathrm{eff}}\nabla {\mathsf{\Phi }}_{\mathrm{e}}\right)=-{\mathrm{i}}_{\mathrm{v}} \mathrm{CLs} \end{gathered}$$

(11)

$$\nabla ·{\mathrm{i}}_{\mathrm{i}}=\nabla ·\left(-{\mathsf{\sigma }}_{\mathrm{i}}\nabla {\mathsf{\Phi }}_{\mathrm{i}}-\frac{{\mathsf{\sigma }}_{\mathrm{i}}\mathsf{\xi }}{{\mathrm{z}}^{-}\mathrm{F}}\nabla {\mathsf{\mu }}_{{\mathrm{H}}_2\mathrm{O},\mathrm{i}}\right)=0 \mathrm{membrane}$$

(12)

$$\nabla ·{\mathrm{i}}_{\mathrm{e}}=\nabla ·\left(-{\mathsf{\sigma }}_{\mathrm{e}}^{\mathrm{eff}}\nabla {\mathsf{\Phi }}_{\mathrm{e}}\right)=0 \mathrm{PTLs}$$

(13)

where $\mathrm{i}$ is the current density. The subscripts i and e denote ionic and electronic, respectively. $\mathsf{\Phi }$ and $\mathsf{\sigma }$ represent the potential and conductivity, respectively. Here, $\mathsf{\xi }$ and ${\mathrm{z}}^{-}$ are the electroosmotic coefficient and charge of hydroxide ions (Table 3) and ${\mathrm{i}}_{\mathrm{v}}$ is the volumetric reaction current. The kinetics will be discussed in more detail in the following section. The second term of the conservation equation of the ion transport corresponds to the contribution of the chemical potential of water absorbed in the ionomer phase, ${\mathsf{\mu }}_{{\mathrm{H}}_2\mathrm{O},\mathrm{i}}$ [20]. Finally, the conductivities in the porous layers were corrected using the Bruggeman correlation to consider the tortuosity.

$${\mathsf{\sigma }}_{\mathrm{i}}^{\mathrm{eff}}={\mathsf{\sigma }}_{\mathrm{i}}^{\mathrm{bulk}} {\mathsf{\epsilon }}_{\mathrm{i}}^{1.5}$$

(14)

$${\mathsf{\sigma }}_{\mathrm{e}}^{\mathrm{eff}}={\mathsf{\sigma }}_{\mathrm{e}}^{\mathrm{bulk}}{ \mathsf{\epsilon }}_{\mathrm{e}}^{1.5}$$

(15)

Here, the superscripts eff and bulk denote the effective and bulk conductivity, respectively. ${\mathsf{\epsilon }}_{\mathrm{i}}$ and ${\mathsf{\epsilon }}_{\mathrm{e}}$ are the volume fraction of the ionomer phase and electronically conductive phase (Table 2).

2.4. Electrochemical Reactions

Water is reduced during the hydrogen evolution reaction (HER) at the cathode while te oxygen evolution reaction (OER) occurs at the anode. The OER and HER kinetics were implemented into the model according to the Butler–Volmer equation.

$${\mathrm{i}}_{\mathrm{v},\mathrm{a}}={\mathrm{s}}_{\mathrm{l}} {\mathrm{a}}_{\mathrm{v}} {\mathrm{i}}_0^{\mathrm{OER}}\left[{\left(\frac{{\mathrm{c}}_{\mathrm{t}}}{{\mathrm{c}}_{{\mathrm{OH}}^{-},\mathrm{ref}}}\right)}^1\exp \left(\frac{\mathsf{\beta }\mathrm{F}}{\mathrm{RT}}\mathsf{\eta }\right)-{\left(\frac{{\mathrm{p}}_{{\mathrm{O}}_2}}{{\mathrm{p}}_{\mathrm{ref}}}\right)}^1\exp \left(-\frac{\mathsf{\alpha }\mathrm{F}}{\mathrm{RT}}\mathsf{\eta }\right)\right] \mathrm{at\; the\; anode}$$

(16)

$${\mathrm{i}}_{\mathrm{v},\mathrm{c}}={\mathrm{s}}_{\mathrm{l}} {\mathrm{a}}_{\mathrm{v}} {\mathrm{i}}_0^{\mathrm{HER}}\left[{\left(\frac{{\mathrm{c}}_{\mathrm{t}}}{{\mathrm{c}}_{{\mathrm{OH}}^{-},\mathrm{ref}}}\right)}^1{\left(\frac{{\mathrm{p}}_{{\mathrm{H}}_2}}{{\mathrm{p}}_{\mathrm{ref}}}\right)}^{0.5}\exp \left(\frac{\mathsf{\beta }\mathrm{F}}{\mathrm{RT}}\mathsf{\eta }\right)-\exp \left(-\frac{\mathsf{\alpha }\mathrm{F}}{\mathrm{RT}}\mathsf{\eta }\right)\right] \mathrm{at\; the\; cathode}$$

(17)

where the subscripts a and c denote anodic and cathodic, respectively. $\mathsf{\alpha }$ and $\mathsf{\beta }$ represent the charge transfer coefficient for oxidation and reduction, both assumed to be equal to 0.5. $\mathrm{F}$ is Faraday’s constant, R—universal gas constant, and T—temperature. Note that the volumetric current is positive at the anode and negative at the cathode. ${\mathrm{a}}_{\mathrm{v}}$ refers to the active catalyst surface area per unit volume of the catalyst layer and ${\mathrm{i}}_0$, the exchange current densities. Both the anodic and cathodic exchange current densities are the ones measured for commercial IrO₂ (Premion^®) and PtRu/C (HiSPEC^® 12100) catalysts [20] and depend on the hydroxide concentration, ${\mathrm{c}}_{{\mathrm{OH}}^{-}}$ (

Table 4. Parameter values as a function of the hydroxide concentrations taken from [20].

Symbol	Parameter	Values					Unit
${\mathrm{c}}_{{\mathrm{OH}}^{-}}$	Hydroxide concentration	0.01	0.02	0.1	0.5	1	mol/L
${\mathsf{\sigma }}_{\mathrm{i}}$	Ionic ionomer conductivity	23.8	32.2	40.8	51.1	57.4	mS/cm
${\mathrm{i}}_0^{\mathrm{OER}}$	Exchange current density for OER	0.00242	0.00354	0.0129	0.0172	0.0256	mA/cm2
${\mathrm{i}}_0^{\mathrm{HER}}$	Exchange current density for HER	0.3439	0.4523	0.6067	0.7894	0.8885	mA/cm2
${\mathrm{a}}_{\mathrm{v}}$	Electrocatalyst specific surface area 1	1.3 × 107	1.3 × 107	1.3 × 107	1.3 × 107	1.3 × 107	m2/m3

¹ Fitted for ${\mathrm{c}}_{{\mathrm{OH}}^{-}}=1 \mathrm{mol}/\mathrm{L}$ for validation by comparing the numerical polarization curves with experimental data obtained in ref. [20] at high current.

). ${\mathrm{c}}_{\mathrm{t}}$ is the concentration of the fixed positive charge in the ionomer, which determines the hydroxide concentration to maintain charge electroneutrality.

$${\mathrm{c}}_{\mathrm{t}}=\frac{{\mathsf{\rho }}_{\mathrm{i}}}{\mathrm{EW}}{\mathsf{\epsilon }}_{\mathrm{i}}\hspace{1em}\mathrm{with}\hspace{1em}\mathrm{EW}=1/\mathrm{IEC}$$

(18)

$${\mathrm{z}}^{-}{\mathrm{c}}_{{\mathrm{OH}}^{-}}+{\mathrm{z}}^{+}{\mathrm{c}}_{\mathrm{t}}=0$$

(19)

in which ${\mathsf{\rho }}_{\mathrm{i}}$, $\mathrm{EW}$ and $\mathrm{IEC}$ represent the density of the ionomer, the equivalent weight of the polymer electrolyte, and the ion exchange capacity, respectively (Table 3). ${\mathrm{z}}^{+}$ is the positive charge. $\mathsf{\eta }$ is the overpotential in the Butler–Volmer equation, which is evaluated from the electronic and ionic potentials.

$$\begin{gathered} \mathsf{\eta }={\mathsf{\Phi }}_{\mathrm{e}}-{\mathsf{\Phi }}_{\mathrm{i}}-{\mathrm{E}}_{\mathrm{eq}} \\ {\mathrm{with} \mathrm{E}}_{\mathrm{eq}}=\{\begin{array}{cc}1.23-0.9\times {10}^{-3}\left(\mathrm{T}-298\right)& \mathrm{at\; the\; anode}\hfill \\ 0& \mathrm{at\; the\; cathode}\hfill \end{array} \end{gathered}$$

(20)

where ${\mathrm{E}}_{\mathrm{eq}}$ and $\mathrm{T}$ are the equilibrium potential and temperature. The volume fraction of the liquid phase is considered in the Butler–Volmer equation to consider the reduction in active surface area caused by gas production [22].

2.5. Water Transport through the Ionomer

Water transport in the ionomer is governed by electroosmosis, ion migration, and diffusion [20,23].

$${\mathrm{N}}_{{\mathrm{H}}_2\mathrm{O},\mathrm{i}}=-\frac{{\mathsf{\sigma }}_{\mathrm{i}}\mathsf{\xi }}{{\mathrm{z}}^{-}\mathrm{F}}\nabla {\mathsf{\Phi }}_{\mathrm{i}}-\left({\mathsf{\alpha }}_{\mathrm{l}}+\frac{{\mathsf{\sigma }}_{\mathrm{i}}{\mathsf{\xi }}^2}{{\mathrm{z}}^{-2}{\mathrm{F}}^2}\right)\nabla {\mathsf{\mu }}_{{\mathrm{H}}_2\mathrm{O},\mathrm{i}}$$

(21)

where ${\mathrm{N}}_{{\mathrm{H}}_2\mathrm{O},\mathrm{i}}$ represents the water flux through the ionomer and ${\mathsf{\mu }}_{{\mathrm{H}}_2\mathrm{O},\mathrm{i}}$ represents the chemical potential of water absorbed in the ion-conducting phase. Here, ${\mathsf{\alpha }}_{\mathrm{l}}$ corresponds to the liquid-equilibrated water transport coefficient [22].

$${\mathsf{\alpha }}_{\mathrm{l}}=\frac{{\mathrm{k}}_{\mathrm{sat}}}{{\mathsf{\mu }}_{\mathrm{l}}{\overline{\mathrm{V}}}_{\mathrm{l}}^2}{\mathsf{\epsilon }}_{\mathrm{i}}^{1.5}$$

(22)

in which ${\overline{\mathrm{V}}}_{\mathrm{l}}$ is the molar volume of liquid water. A full liquid equilibration is considered, which means that the volume fraction of absorbed water in the ionomer is equal to its maximum value, i.e., for water activity equal to unity. The resulting conservation equation for water transport through the ionomer is presented in Equation (23).

$$\nabla ·{\mathrm{N}}_{{\mathrm{H}}_2\mathrm{O},\mathrm{i}}={\mathrm{R}}_{{\mathrm{H}}_2\mathrm{O}}-{\mathrm{R}}_{\mathrm{ml}}$$

(23)

where ${\mathrm{R}}_{{\mathrm{H}}_2\mathrm{O}}$ represents the rate of water production at the anode and water consumption at the cathode.

$$\begin{gathered} {\mathrm{R}}_{{\mathrm{H}}_2\mathrm{O}}=+\frac{{\mathrm{i}}_{\mathrm{v},\mathrm{a}}}{2\mathrm{F}} \mathrm{at} \mathrm{the} \mathrm{anode} \\ {\mathrm{R}}_{{\mathrm{H}}_2\mathrm{O}}=+\frac{{\mathrm{i}}_{\mathrm{v},\mathrm{c}}}{\mathrm{F}} \mathrm{at} \mathrm{the} \mathrm{cathode} \end{gathered}$$

(24)

The rate of water absorption or desorption from the ionomer to the liquid phase, ${\mathrm{R}}_{\mathrm{ml}}$, is assumed to be proportional to the chemical potential difference:

$${\mathrm{R}}_{\mathrm{ml}}=\frac{{\mathrm{k}}_{\mathrm{ml}}}{{\mathrm{M}}_{{\mathrm{H}}_2\mathrm{O}}}\left({\mathsf{\mu }}_{{\mathrm{H}}_2\mathrm{O},\mathrm{i}}-{\mathsf{\mu }}_{\mathrm{liq}}\right)$$

(25)

in which ${\mathrm{M}}_{{\mathrm{H}}_2\mathrm{O}}$ and ${\mathrm{k}}_{\mathrm{ml}}$ are the molar mass of water and the rate constant for water transfer from the ionomer to the liquid phase. ${\mathsf{\mu }}_{\mathrm{liq}}$ corresponds to the chemical potential of liquid water [22]:

$${\mathsf{\mu }}_{\mathrm{liq}}={\mathrm{H}}_{\mathrm{l}}{\mathrm{M}}_{{\mathrm{H}}_2\mathrm{O}}\left(1-\frac{\mathrm{T}}{{\mathrm{T}}_{\mathrm{t}}}\right)+{{\mathrm{C}}_{\mathrm{p}}}_{\mathrm{l}}{\mathrm{M}}_{{\mathrm{H}}_2\mathrm{O}}\left(\mathrm{T}-{\mathrm{T}}_{\mathrm{t}}-\mathrm{T}\ln \left(\frac{\mathrm{T}}{{\mathrm{T}}_{\mathrm{t}}}\right)\right)+{\overline{\mathrm{V}}}_{\mathrm{l}}\left({\mathrm{p}}_{\mathrm{l}}-{\mathrm{p}}_{\mathrm{t}}\right)$$

(26)

where ${\mathrm{H}}_{\mathrm{l}}$, ${{\mathrm{C}}_{\mathrm{p}}}_{\mathrm{l}}$, ${\mathrm{T}}_{\mathrm{t}}$, and ${\mathrm{p}}_{\mathrm{t}}$ represent the enthalpy and heat capacity of liquid water and the temperature and pressure at the triple point, respectively (Table 3). Note that water vapor is not considered in the model, which means that the gases produced are pure and perfectly dry and that adsorption and desorption only occur between the liquid phase and the ionomer phase.

2.6. Boundary Conditions

Regarding laminar flows through the channels, a total KOH solution flow rate of 100 mL/min for the base case was applied to the inlet of the 24 channels at both the cathode and anode while a relative pressure equal to zero was set at the outlet of each channel (Figure 1). In addition, a free-porous interface boundary condition was used to couple the free and porous media flow between the channels/GDL interface at both sides with a ratio of free flow volume fraction equal to the unity between the phase saturation on the porous side and the phase volume fraction on the free-flow side. As mentioned in hypothesis 3, the no-slip condition is assumed at the free-porous interface on the free-flow side, which means that the tangential component of the mixture velocity is set to zero. No flux is considered on exterior boundaries where no mass flows in or out of the boundaries, which means that the total flux is zero. Finally, the electronic potential ${\mathsf{\Phi }}_{\mathrm{e}}$ was set to 0 V at the cathode land/GDL boundary, i.e., between the channels at the interface of the cathode GDL and bipolar plate (not represented in the model) and to the applied cell potential (from 1.4 to 2 V) at the anode land/GDL boundary, i.e., between the channels at the interface of the anode GDL and bipolar plate (not represented in the model).

2.7. Numerical Procedure

Figure 1 shows the computational domain of the AEMWE model, its components and the mesh configuration. The AEM cell of 1.2 × 1.5 cm² is composed of 24 straight channels on both sides for supplying KOH solution, two porous transport layers (PTLs) and catalyst layers (CLs), and one anion exchange membrane. This cell geometry is the one we are developing for our future experiments. The model dimensions are presented in

Table 5. Dimensions of the modeling domains.

Parameter	Values	Unit
Cell width	12	µm
Channel length	15	µm
Channel height	300	µm
Channel width	250	µm
Rib width	250	µm
CPTL thickness	190	µm
APTL thickness	270	µm
CL thickness	10	µm
Membrane thickness	50	µm

. All of the simulations were conducted at 60 °C and ambient pressure for different KOH concentrations from 0.01 to 1 M with a flow rate of 100 or 10 mL/min to both sides. All of the governing equations were solved by using version 6.2 of COMSOL Multiphysics^® software. A residual gas saturation of 1% was assumed in all layers for modeling convenience. The mesh is composed of 46,880 hexahedral elements. The calculation time is about 2 h 30 on Intel^® Core^TM i5-10505 CPU @ 3.20 GHz to obtain one polarization curve. A grid-independent validation was performed by using different level of finite element mesh refinement and similar results were obtained. The grid-independent validation can be found in Supplementary Material.

3. Results and Discussion

3.1. Model Calibration and Validation

The ability of the 3D AEMWE cell model to simulate the impact of KOH on cell performance is first validated through comparison of the simulation results with the experimental data published by Liu et al. [20]. The polarization curves computed with the 3D model (solid lines) and experimental data (symbols) obtained for various KOH concentrations ranging from 0.01 to 1 M at 60 °C and 1 atm are presented in Figure 2. The specific surface area was numerically calibrated to fit the experimental results obtained by Liu et al. [20] for a KOH concentration of 1 M (purple symbols in Figure 2) at high current (Table 3). The same specific surface area is considered for the other KOH concentrations as this parameter does not change with the KOH concentration. Overall, the influence of KOH concentration on the polarization curves is well captured by the simulations well. Nevertheless, if the simulated polarization curve obtained for a KOH concentration of 0.5 M (solid blue line in Figure 2) slightly deviates from the experimental data at low-current density, quite good agreement with the measurements is achieved at high current density for the other KOH concentrations. Indeed, for low current, our 3D model overestimates the performance obtained by Liu et al. [20] except for a KOH concentration of 0.01 M (solid orange line in Figure 2). The deviation between our 3D model and data from Liu et al.’s study can be attributed to differences in cell geometry, namely the triple serpentine flow-fields used by Liu et al. [20].

Figure 3 shows the total volumetric reaction current as well as the ionic and electronic potential distributions from the catalyst layer/membrane interface through the anode and cathode CL thickness at the middle of the cell length (green dashed line in Figure 1) at a cell current density of 0.5 A/cm² (dashed line in Figure 2) for different electrolyte concentrations. As the KOH concentration decreases, the local volumetric current density and ionic potential gradients increase while the electronic potential distribution remains almost uniform over the entire CL thickness. Volumetric current densities and ionic potential both exhibit higher gradients close to the catalyst layer/membrane interface. Indeed, ionic ohmic drop limitation becomes predominant when lowering the KOH concentration, reducing the catalyst utilization rate while the electronic conductors act as an equipotential. Under a KOH flow rate of 100 mL/min and at an applied current density of 0.5 A/cm², the shapes of the current density distribution are quite different from the ones predicted by Liu et al.’s 1D model [20] except at high KOH concentration. This difference is likely due to the fact that we assume a constant electrochemical surface area whatever the KOH concentration. It is also worth mentioning that the current density distributions do not change significantly along the channel length from the inlet to the outlet as the flow rate is quite high in this first case. Moreover, no rib/channel effects were noticed as the straight channel shape was considered.

Figure 4a shows the H₂ and O₂ gaseous volume fraction distribution at 2 V (i.e., at the maximum current density) for all of the KOH concentrations in porous media (i.e., CL and PTL) and gas channels. This latter is displayed along the cell length thanks to a through-plane cut in the middle of the cell (i.e., the red zone in Figure 1). The cuts respect an aspect ratio of one-quarter in length for better visibility with the inlet flow of the KOH solution at the bottom and the outlet at the top. The cathode compartment (channel-CPTL-CCL) is represented on the left of each cut, where HER takes place (negative electrode), and the anodic compartment (ACL-APTL-channel) on the right, where OER occurs (positive electrode). At low concentrations, the gas volume fraction is truly uniform along the porous layers and the gas channel because the current density is low, which consequently results in less gas production. It becomes heterogenous at higher current densities when the KOH concentration increases. Finally, Figure 4b shows the gaseous volume fraction at different cell voltages at 1 M KOH. As in the previous case, gas production becomes more heterogeneous when increasing the cell voltage (i.e., increasing the current density) because of the gaseous volume fraction accumulation along the cell length on the flow direction.

Liu et al.’s model was extended to be able to forecast the through-plane and in-plane gas distribution in the different layers of the cell. At a flow rate of 100 mL/min, the gas volume fraction distribution becomes strongly non-uniform at a current density above 1 A/cm² (corresponding to a cell voltage of 1.8 V at 1 M and 2 V at 0.1 M). The gas phases accumulate both in the PTL and the channel. The gas volume fraction is also more heterogeneous at the cathode (H₂) than at the anode (O₂). Regardless of the concentration, oxygen tends to stagnate more than hydrogen in its respective porous media. This can be explained by the higher density of oxygen compared to hydrogen, as well as an anode PTL that is one-third thicker (but with the same porosity) than the cathode PTL (Table 2 and Table 5). In the channel, the gas phases accumulate along the cell length from the inlet toward the outlet because of the electrolyte flows. Nevertheless, the impact of the two-phase flow is quite small when considering a high electrolyte flow rate. Therefore, a flow rate ten times lower is considered in the next section.

3.2. Flow Rate Study

Figure 5 shows the total volumetric reaction current through the anode catalyst layer and cathode catalyst layer at the middle of the cell length (the green dashed line in Figure 1) at a cell current density of 0.5 A/cm² for two different flow rates, namely 100 mL/min (solid curves) and 10 mL/min (dashed curves). Even if the flow rate of the KOH solution is ten times lower, the volumetric reaction current through the anode and cathode catalyst layers’ thickness at a cell current density of 0.5 A/cm² remains similar to the one obtained at a higher flow rate. Indeed, the bubble overpotential contribution is not so important when noble catalysts are used. The electrochemical reaction kinetics are fast enough to maintain performance. Liu et al.’s simulations confirm a quite small bubble overpotential lower than 50 mV even at high-current density [20].

To go further in our understanding of flow rate impact, Figure 6 shows the polarization curves for different KOH concentrations and the gaseous volume fraction for 1 M at 1.55 A/cm² for flow rates of 10 mL/min and 100 mL/min. Surprisingly, the performance of the AEMWE cell does not deteriorate so much when dividing the flow rate by ten (Figure 6a). Indeed, the cell voltage change does not exceed 100 mV in the worst situation (at 1.5 A/cm² which can be achieved only in 1 M KOH) despite the very inhomogeneous distribution of the gaseous volume fraction distribution along the channel length at a lower flow rate (Figure 6b). At a lower electrolyte flow rate, the gaseous volume fraction can reach high values at the channel outlet (0.7 at the cathode and 0.5 at the anode at a KOH concentration of 1 M) with a relative effect on cell performance. Moreover, the larger conversion ratio of water to H₂ and O₂ also induces a strongly non-uniform current density distribution through the catalyst layer thicknesses (Figure 7a,b). Nonetheless, these heterogeneities remain relatively small as the anode and cathode overpotential partly contribute to cell polarization when PGM catalysts are used. Indeed, the fast HER kinetics on the Pt catalyst are only slightly influenced by the decrease in the active surface area caused by the H₂ gaseous volume fraction. This non-uniform distribution may be more pronounced if PGM-free catalysts are used because of the necessity of a thicker catalyst layer in order to compensate for lower catalyst activity. Figure 7c,d displays the gaseous volume fraction for flow rates of 10 mL/min through the catalyst layer thickness. As the catalyst layer is thin with a PGM catalyst, no significant gradient can be observed even at a low flow rate.

4. Conclusions

AEMWEs have the potential to be an efficient and environmentally friendly method for producing hydrogen fuel. Nonetheless, they face certain constraints that could limit their practicality in certain scenarios. Hence, modeling proves to be a valuable approach for understanding and optimizing the performance of these systems.

In this work, a 3D two-phase flow multi-physics model of an anion exchange membrane water electrolyzer is developed to investigate the concentration influence of a potassium hydroxide electrolyte. The model was validated through comparison with Liu et al.’s experimental data obtained at different concentrations [20]. The simulated AEMWE polarization curves for various KOH concentrations are in good agreement with the experimental data. In the case of the 10 µm thick catalyst layer, the model predicts the appearance of strongly non-uniform H₂ and O₂ gas volume distributions for current densities above 1 A/cm². The impact of H₂ and O₂ gas production and transport was investigated. The performance of AEMWE cells varies locally along the channel due to changes in gas volume fraction distribution, namely at high current density and low electrolyte flow rate, reducing the catalyst utilization rate. According to the simulations performed in 1 M KOH flowing at 10 mL/min and at an applied current density of 1.55 A/cm², the gas volume fraction reaches 0.7 at the cathode (H₂) and 0.5 at the anode (O₂). Nevertheless, in the case of the PGM catalyst, the impact of the gas volume fraction is limited due to the fast electrochemical HER and OER reactions as well as the thin catalyst layer. The influence of the two-phase flow distribution will likely become more detrimental when PGM-free catalyst layers are used due to the use of thicker electrodes. Bubble overpotential may strongly increase, reducing catalyst layer utilization and cell performance. Therefore, the influence of bubble formation on the specific active surface area needs to be better understood. This topic is very tricky in alkaline membrane electrolyzers and requires intense research.