Three-Dimensional CFD Simulation of a Proton Exchange Membrane Electrolysis Cell

: The energy shift towards carbon-free solutions is creating an ever-growing engineering interest in electrolytic cells, i.e., devices to produce hydrogen from water-spli tt ing reactions. Among the available technologies, Proton Exchange Membrane (PEM) electrolysis is the most promising candidate for coping with the intermi tt ency of renewable energy sources, thanks to the short transient period granted by the solid thin electrolyte. The well-known principle of PEM electrolysers is still unsupported by advanced engineering practices, such as the use of multidimensional simulations able to elucidate the interacting ﬂ uid dynamics, electrochemistry, and heat transport. A methodology for PEM electrolysis simulation is therefore needed. In this study, a model for the multidimensional simulation of PEM electrolysers is presented and validated against a recent literature case. The study analyses the impact of temperature and gas phase distribution on the cell performance, providing valuable insights into the understanding of the physical phenomena occurring inside the cell at the basis of the formation rate of hydrogen and oxygen. The simulations regard two temperature levels (333 K and 353 K) and the complete polarization curve is numerically predicted, allowing the analysis of the overpotentials break-up and the multi-phase ﬂ ow in the PEM cell. An in-house developed model for macro-homogeneous catalyst layers is applied to PEM electrolysis, allowing independent analysis of overpotentials, investigation into their dependency on temperature and analysis of the cathodic gas–liquid strati ﬁ cation. The study validates a comprehensive multi-dimensional model for PEM electrolysis, relevantly proposing a methodology for the ever-growing urgency for engineering optimization of such devices


Introduction
The increasing focus on climate change and global warming has induced governments of industrialised countries to radically review future energy production and management strategies. Decarbonisation and environmental sustainability policies aim at promoting the research and development of alternative technologies, combining increasing energy demand with a net reduction of environmental impact. The use of renewable energy sources (RES) is a long-term strategy, requiring not only the capacity of RES to fulfil the energy production request but also the development of energy storage technologies, shaving the peaks of energy overproduction [1,2]. Energy chemical storage in the form of hydrogen production unveils several possibilities, being an energy buffer and carrier that can be used both as storage as well as in almost every application to replace fossil fuels. However, hydrogen must be produced in its molecular form (H2), which is extremely rare in nature, and hence it must be obtained by the separation of hydrogen-containing compounds (e.g., water or methane). The achievement of short-term decarbonisation requires not only the development of technologies to produce massive amounts of hydrogen, but also a shift in the leading technology away from the to-date ubiquitous use of Steam Methane Reforming (SMR) and coal gasification, which jeopardise both the decarbonisation and energy independence goals. The electrolysis of water powered by RES-derived electricity is the only way to produce carbon-free hydrogen, also named "Green Hydrogen". This process is characterized by sustainability and high efficiency and production volumes. Among the available technologies, the Proton Exchange Membrane Electrolyser Cell (PEMEC) is the most suitable to accommodate the highly unstable and fluctuating energy supply from RES overproduction to store it in the form of hydrogen, although technology advances are still needed to lower the specific cost of hydrogen production and to improve the profitability of PEM-based electrolysis plants [3].
A water electrolyser is an electrochemical device that splits water molecules into gaseous hydrogen and oxygen, thus directly converting electrical into chemical energy. Its operation is based on a complex interaction between fluid/thermal and electrochemical processes, such as the multi-phase flow field's evolution of the reactants supply and the removal of products from the reaction sites, with a direct effect on the electrochemical reaction rate and the required energy for the process. Also, the interaction of multiple parts requires a fundamental understanding of how the different materials interact with transport processes, such as heat rejection and the current of charges (both electric and ionic). A detailed knowledge of such aspects allows optimizing PEMEC design, thus reducing the amount of required catalyst or the formation of local stresses, for improved PEMEC economy with minimal degradation issues.
Numerical simulation is a powerful tool that allows the study and quantifying of all the mentioned phenomena [4]. This can be done on different levels, from zerodimensional models [5,6] to more complex three-dimensional ones. Multi-dimensional Computational Fluid Dynamic (CFD) models belong to the last class, allowing the simulation of the interaction of the heat/fluid/current transport and the distribution of each variable, thus complementing experiments for difficult to measure quantities. Given the analogy with the physics involved in fuel cells operation, CFD has been recently used to deepen the understanding of the operation of Proton Exchange Membrane Fuel Cells (PEMFC) [7]. However, the number of simulation studies on PEMEC is noticeably lower, hence leaving a critical knowledge gap on the topic. Among the available studies, Nie et al. [8] published a 3D-PEMEC model based on a mixture multi-phase approach [9]. The study refers to the anodic water flow and the two-phase flow field in the channel is analysed in terms of pressure and velocity distribution, to improve the design of bipolar plates. However, the electrochemical reaction is replaced by a constant and uniform oxygen production along the MEA/channel interface. They concluded that parallel configuration is not an optimal design for bipolar plates, because of the high nonuniformity of the oxygen concentration. Ruiz et al. [10] developed a high-temperature PEMEC model operating above 100°C, studying the impact of three different channel configurations on the cell performance (parallel, single-and multi-path serpentine). They analysed the thermal field, the distribution of pressure in the channels, the current density, and the molar fraction of hydrogen within the cell, asserting that the multi-path serpentine design shows the best performance. In this work, all the components of a PEMEC are considered, and mass and heat transfer within each media and the electrochemical kinetics at the electrodes are taken into account. However, due to the high working temperature, the presence of the liquid phase is completely neglected, so they adopt a single-phase approach. Similarly, Toghyani et al. [11][12][13] developed a single-phase 3D-CFD model finalised to study different channel configurations, obtaining results in line with [10] and confirming that the multi-serpentine design leads to the best performance. Porous media were modelled with a macro-homogeneous approach, i.e., replacing the complex porous micro-structure with statistical parameters (e.g., porosity, etc.). Despite the operating temperature being 353 K, the presence of liquid water is critically neglected, justified by considering the anodic two-phase flow as being more important at current densities higher than that considered. Ma et al. [14] proposed a 3D-CFD study of a triple serpentine PEMEC based on the mixture multiphase model, considering the development of flow, thermal, and electrochemical fields. Several processes are considered, such as electro-osmotic drag, membrane water diffusivity, and water evaporation. However, catalyst layers (CLs) are not directly modelled and they are replaced by a reacting interface. As commonly done in PEMFC models [15], this modelling simplification reduces the computational cost, but it also introduces a relevant simplification of the electrode physics [16]. Lafmejani et al. [17] used a Volume of Fluid (VOF) multi-phase model, numerically tracking the inter-phase surface and aiming at simulating the superficial motion. This was used to describe the unsteady gas/liquid anodic flow, including the gas bubbles morphology and the gas/liquid superficial interaction. Results provide a detailed description of the two-phase distribution in the channel, although modelling assumptions include isothermal processes, the neglect of capillary effects in porous media, phase transition, and electrochemistry. Another valuable VOF study of PEMEC is in [18] where an iterative procedure links the two-phase VOF method to the 3D full-cell model. The anodic gas/liquid phase evolution is studied utilizing the VOF method, obtaining the volume fraction distribution at the channel/electrode interface and using this information as a boundary condition for the full-cell model, resolving the oxygen pressure equation in the porous medium. Nevertheless, the computational cost is very high due to the reduced time step required.
The present work presents a modelling methodology for the three-dimensional fluiddynamic simulation of a PEMEC. The model is a non-isothermal and multi-physics approach to simulate the multi-phase flows (using the mixture multi-phase model) in the domain of a PEMEC. The study relevantly extends the use of a recently developed macrohomogeneous model of CLs [19], applied here to the reproduced PEMEC model and experimental data from [14] for several voltages and at two temperature levels, thus allowing us to exhaustively elucidate all the crucial phenomena at the basis of PEMEC operation. The 3D modelling of these devices allows for a comprehensive view of the concurring processes (e.g., distribution of the gas phase, temperature field, etc.), and their respective effects on cell performance, thus quantifying hard to measure quantities. This promotes the use of the presented model as a key step of the design and optimization phases of these devices.

Mathematical Modelling
The mixture multi-phase model (MMP) method is used in this study to consider the presence of multiple phases, where it is assumed that these are miscible and at equilibrium, and that their motion can be simulated as that of a unique continuum. Continuity, momentum, and energy equations are solved for the Eulerian mixture (e.g., the mixture density or mixture velocity vector are considered), and the phase subdivision is handled by a dedicated transport equation for the volume fraction [20][21][22]. The governing equations for the MMP model (Equations (1)-(6)) are synthesised in Table  1, together with the region-dedicated source term specification. In this study, the governing equations are simplified using the following modelling assumptions: 1. Laminar flow regime, because of the low fluid velocity both in channels and porous diffusion media. The calculated Reynolds number in the examined cases is approximately 1000 in the channels, confirming this assumption. 2. Ideal gases behaviour, given the relatively low pressure and temperature. 3. The gravitational force is not considered. 4. The gas diffusion layer (GDL) and CL are treated as isotropic and homogeneous porous media, characterized by effective permeability, uniform porosity, and tortuosity. 5. Butler-Volmer kinetics govern the electrochemical reaction at the anode and cathode. 6. The membrane is an impermeable solid medium, and the water flux is modelled by sorption reactions at the interfaces adjacent to CLs.
7. The simulations are steady state since the objective is to analyse the time-independent cell's performance at different voltages, aiming at understanding its stationary operation and numerically reproducing the polarization curve.

Governing Equation Source Term Specification
Continuity equation:

Porous Media Modelling
The chosen modelling approach involves neglecting the complex fibrous structure of the material, which is substituted by a macro-homogeneous approach. This allows consideration of the co-presence of both a liquid and a solid phase, and dedicated variables are introduced with this aim, such as porosity, tortuosity, permeability, etc. For the calculation of tortuosity, the equation proposed in [23] is used: As discussed in [19], the flow resistance exerted by the porous medium is expressed in Equation (8), where a viscous and an inertial contribution are present. Given the low velocities involved, the inertial term can be neglected, whereas the viscous term of the nth phase is calculated as in Equation (9) [19].
where ( ) , , and are the dynamic viscosity, the material permeability, and the relative permeability, respectively. Equation (10) with k = 4 is used to calculate the latter, in agreement with [24].
The contribution of capillary force is considered through a vector source term in the momentum equation, formulated as in Equation (11), including the gas-liquid surface tension () and the fluid/wall contact angle ( = 110°) and expressing the momentum exerted by the gradient of the capillary pressure. Here it is determined using a Leverett's function because of the lack of experimental data, although it has been derived for largely different porous media (homogeneous solid and sandy beds) as discussed in [25].

Membrane Modelling
The polymeric membrane is modelled as an impermeable solid electrolyte, preventing the crossover of fluids and allowing ion current. The ionic conductivity is calculated using the correlation proposed by Springer et al. [26] (Equation (12), justified by the analogy with PEMFC where the same perfluorosulfonic polymer (DuPont ® Nafion, USA) is used. Differently than in PEMFC technology, where reactant humidification techniques are mandatory for sufficient power generation [27], in the PEMEC simulation the membrane water content (λ) is imposed at λ = 22, i.e., a fully hydrated condition. This is justified by the massive presence of liquid water in contact with the membrane. This same assumption is made in [14].
The voltage to be applied to sustain the PEMEC reactions depends both on the opencircuit potential ( ) and on several potential losses as in Equation (15), i.e., , , , (cathodic and anodic activation overpotentials), (concentration overpotential), and (ohmic overpotential). Open-circuit voltage and activation overpotentials are calculated as in Equations (16)-(18) [28][29][30], where / are the charge transfer coefficients (symmetrical and equal to 1), is the Faraday constant, the ideal gas constant, the current density, and , / the anodic/cathodic exchange current density.
Exchange current densities are calculated as indicated in [14] to ensure validation (Equations (19) and (20)), whereas the ohmic overpotential is obtained by the thermal power dissipated by the Joule effect ( ) and the current ( ) in the different components of the cell (Equation (21) Finally, the cathodic/anodic volumetric current densities ( / ) are calculated using Butler-Volmer equations (Equations (22) and (23) The CL is modelled adopting a macro-homogeneous approach, which is implemented in the commercial CFD code SIMCENTER STAR-CCM+, licensed by SIEMENS DISW, via user-coding. This approach allows the three-dimensional modelling of the CL, which is represented as a finite-thickness part, and it ensures the development of electrochemical reactions on the unresolved internal surface. Despite not being geometrically detailed, it is estimated by means of average parameters, as reported in [19]. The solid phase of the CL consists of both an ionomeric solid phase, equivalent to that used for the membrane and deputed to ion transport, and an electronic solid phase, represented by carbon-supported platinum and used for the electron transport.

3D-CFD Model
The model of a PEMEC with a three-channel serpentine proposed by Ma et al. in [14], whose geometry and dimensions are shown in Figure 1a-c and Table 2, respectively, is reproduced in SIMCENTER STAR-CCM+ v2022.2 (SIEMENS DISW). This is a cell with a triple serpentine channel configuration measuring 5 cm × 5 cm. The number of computational cell layers for MEM, CL, and GDL are 4, 3, and 6, respectively. The finite volume mesh consists of approximately 1.35 M hexahedral cells, of which 120 k are for each CL. This was chosen from mesh independence tests until mesh-independent solutions were obtained. However, coarser grids (especially in the axial direction) will be evaluated in future studies based on the guidelines proposed in literature. Numerical simulations were run on a 32-core Linux cluster, with a computational time of approximately 20 h of physical time for each case. Material characterization and boundary conditions are reported in Tables 3 and 4, respectively. The thermal and electrical characterization of the boundary conditions at external walls is that of adiabatic ( ⁄ = 0, with being the wall-normal direction) and insulating ( ⁄ = 0) surfaces, except for the current collectors (top and bottom surfaces in Figure 1) treated as isothermal ( = 333 K/353 K, depending on the simulated operating condition) and isopotential surfaces, with imposed values at the cathode/anode of = 0/ , respectively, with  spanning the 1.4-2.0 V range with a step of 0.1 V, and a refined 0.05 V step for low-voltage operation. The temperature on the outer bipolar plate walls is set to the inlet temperature of water (i.e., 333 K and 353 K for the two tested temperatures), assuming an ideal cooling system behaviour. A face-normal velocity specification is imposed at anodic and cathodic inlets, with values based on the flowrate specification from [14] and reported in Table 4, with the single-phase (liquid) water temperature equal to 333 K or 353 K depending on the simulated condition, whereas ambient pressure is imposed at outlet sections for the formed multi-phase mixture (gaseous hydrogen/liquid water at the cathode and gaseous oxygen/liquid water at the anode, respectively).   Table 3. Physical and transport properties, reproduced from [14].

Results
The validation of the numerical model is achieved by comparing the simulated polarization curve with the experimental data reported in [14], for two temperature levels. As it can be seen in Figure 2, the presented results show a good agreement in terms of current/voltage for both the analysed operating temperatures. The experimentally observed voltage reduction with temperature is accurately reproduced by simulations and it is attributed to the increase in the exchange current density (Equations (19) and (20)) and in the membrane ionic conductivity (Equation (12)) with temperature. This results in a reduction of the activation and ohmic overpotentials, respectively, and the obtained results confirm a physics-based modelling of overpotentials contributions. To gain more insight into the performance, an analysis of the individual overpotentials is shown for each operating voltage in Figure 3. The largest contribution is the anodic activation overpotential, motivated by the reduced exchange current density value. The lowest loss is constituted by the mass transport overpotential, negligible for operations below 1.6 V and noticeable only around 2.0 V. This is due to the progressive filling of anodic CL with gaseous oxygen, obstructing the transport of reactants due to its high specific volume (gaseous) and the flow obstruction exerted by porous materials (CLs and GDLs). As might be expected, the ohmic overpotential tends to increase as the current supplied to the electrolyser increases. Figure 4 shows the volume fraction of gas at the interface between the GDL and anode-side channels. The horizontal bands with a reduced volume fraction of gas near the serpentine bends are due to the entry of liquid water into the gas diffusion layer, which tends to push the gas away and eject it into the adjacent channel, where bands with higher values can be seen. This phenomenon is caused by the pressure gradient between these adjacent channel pairs, drawing the liquid to bypass the channel path through the GDL. This phenomenon is explained in detail below.   Figure 5 shows the distribution of the volume fraction of gas in the sections within the GDL, highlighting the gas accumulation areas and the liquid by-pass flow. The x-direction component of the mixture velocity, whose negative values confirm the fluid bypass trend through the GDL under the BP in the indicated areas, is shown in Figure 6. The effect of this phenomenon decreases as the gas proceeds along the channel due to the reduction of the pressure gradient between adjacent channels, as will be explained in detail below.  The analysis of the gaseous volume fraction is extended to the gas channels, as the large specific volume of the gas-phase products alter the volumetric flow rate and hence the flow regime. Figure 7 shows the gas volume fraction distribution on the anodic and cathodic side of the cell, for two operating voltages (1.6 V and 2.0 V). As the voltage increases (i.e., for higher electric power), more oxygen and hydrogen are produced. As reactions develop on all the active surface, the gas phase tends to increase towards the outlet of the electrolysis cell, making gas removal from CLs and GDLs increasingly complex. The larger presence of gas in the cathode side is motivated by the lower inlet flowrate, as reported in Table 4, promoting the gas-phase build up. Figure 8 shows the gas volume fraction on transverse sections, highlighting the local stratification of the gas phase in the anodic CL and GDL and the consequential filling of the channels along the path.  From a comparison with the gas volume fraction fields reported in the reference article [14], a significant difference is noticed. The reference article shows higher gas volume fraction values compared to those obtained from the current model, which led to the decision to perform a quantitative analysis of the produced hydrogen flow rate based on Faraday's law to ensure the accurate prediction of the gas phase produced by this model. These results are presented in the final part of this section to validate the presented results. An interesting consequence of the multi-phase nature of flows in a PEMEC is related to the gas phase production and the consequent liquid fraction reduction, leading to a lower mixture density which induces a flow acceleration to conserve the mass flow rate. Figure 9a,b shows the velocity distribution for 1.6 V and 2.0 V, with the latter case showing higher flow velocities deriving from higher gas phase production. This causes an increased pressure drop (c,d) for high voltage operation, resulting in increased energy required for water pumps. Figure 10 shows the temperature ranges at the interface between the GDL and anode-side channels. The average temperature increases with the operating voltage due to the higher amount of power dissipated by the cell, and the areas of maximum temperature coincide with those of maximum gas volume fraction shown in Figure 4.  To complete the modelling methodology validation, the hydrogen production rate ( ) resulting from the 3D-CFD simulations is compared with the theoretical value predicted by Faraday's formula for electrolysis (Equation (24)) [31], which was not discussed or reported in [14]. The results are shown in Figure 11a, confirming that the correct amount of hydrogen is produced at all voltages and temperatures, thus reinforcing the CFD-based detailed analysis presented in the previous paragraphs. Finally, the PEMEC efficiency ( ) calculated as in Equation (25) [31] is reported in Figure 11b. The validation of this type of integral-scale variables (i.e., quantitative relationship between products' flow rate and input current) against theoretical predictions (Faraday's law) is fundamental to reinforce the detailed analysis presented in the previous sections.

Conclusions
In this study, a comprehensive multi-dimensional, multi-physics, and non-isothermal methodology for PEM electrolysis simulation is presented, investigating the interaction between fluid dynamics, electrochemistry, and heat transport. The 3D-CFD model for PEM electrolysis provides a detailed insight into gas-liquid stratification, overpotentials break-up, and the variation of the flow regime as a consequence of gas-phase generation. The model relevantly includes the use of a macro-homogeneous catalyst layer, developed by the authors for PEM fuel cells and here applied to PEM electrolysis.
The model is used to simulate a PEM electrolysis study from the open literature, providing valuable experimental and numerical data for model validation. The polarization curves obtained by the presented model are in excellent agreement with the experimental results at two temperature levels (333 K and 353 K), allowing the investigation of the complex multi-phase flow developing in the distributor channels and in the porous parts. Hydrogen production rates match the theoretically calculated values at all voltages and temperatures, providing a quantitative validation of the obtained model results.
The presented methodology offers an industry-relevant and ready-to-use tool for investigating the complex behaviour of PEM electrolysers. The study contributes to filling a critical engineering gap, supporting the transition towards carbon-free solutions in line with the global migration from fossil-fuel-based sources.