Next Article in Journal
Application of Digital Image Correlation to Evaluate Strain, Stiffness and Ductility of Full-Scale LVL Beams Strengthened by CFRP
Next Article in Special Issue
Elucidating the Mass Transportation Behavior of Gas Diffusion Layers via a H2 Limiting Current Test
Previous Article in Journal
Acoustic Insulation Characteristics and Optimal Design of Membrane-Type Metamaterials Loaded with Asymmetric Mass Blocks
Previous Article in Special Issue
Improved Lead Sensing Using a Solid-Contact Ion-Selective Electrode with Polymeric Membrane Modified with Carbon Nanofibers and Ionic Liquid Nanocomposite
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Three-Dimensional Numerical Simulation of the Performance and Transport Phenomena of Oxygen Evolution Reactions in a Proton Exchange Membrane Water Electrolyzer

1
Department of Aeronautics and Astronautics, Fudan University, Shanghai 200433, China
2
School of Chemical Engineering and Technology, Hainan University, Haikou 570228, China
3
Department of Aeronautics and Astronautics, Zhejiang University, Hangzhou 310027, China
*
Author to whom correspondence should be addressed.
Materials 2023, 16(3), 1310; https://doi.org/10.3390/ma16031310
Submission received: 20 December 2022 / Revised: 26 January 2023 / Accepted: 1 February 2023 / Published: 3 February 2023
(This article belongs to the Special Issue Analysis of Electrode Materials)

Abstract

:
Proton exchange membrane (PEM) water electrolysis, which is one of methods of hydrogen production with the most potential, has attracted more attention due to its energy conversion and storage potential. In this paper, a steady state, three-dimensional mathematical model coupled with the electrochemical and mass transfer physical fields for a PEM water electrolyzer was established. The influence of the different operation parameters on the cell performance was discussed. Moreover, the different patterns of the flow field, such as parallel, serpentine, multi-serpentine, and interdigitate flow fields, were simulated to reveal their influence on the mass transfer and current distribution and how they consequently affected the cell performance. The results of the numerical modeling were in good agreement with the experimental data. The results demonstrated that a higher temperature led to a better mass transfer, current distribution, and cell performance. By comparing the polarization curve, current, velocity, and pressure distribution, the results also indicated that the PEM water electrolyzer with a parallel flow field had the best performance. The results in this study can help in optimizing the design of PEM water electrolyzers.

1. Introduction

Hydrogen is considered to be a clean energy source and a promising replacement for traditional fossil fuel energy [1,2,3]. Proton exchange membrane (PEM) water electrolyzers are a type of device used in the production of hydrogen with the most potential [4,5,6]. Compared with other hydrogen production devices, PEM water electrolyzers have several advantages including a high purity of hydrogen production, high efficiency, compact structure, etc. [7,8,9].
The performance of PEM water electrolyzers is affected by many factors, including operation parameters such as temperature, pressure, reactant flow rate, internal structure parameters such as the thickness and the pore size of the liquid/gas diffusion layer, different flow fields, catalyst layer properties, and morphology [10,11,12]. In order to better understand the mechanism of electrochemical and mass transfer in PEM water electrolyzers and improve their performance, mathematical modeling and simulation can be used as a promising way to theoretically conduct investigations with less time and cost requirements, and this can then be used to optimize the design and improve the performance and efficiency.
Over the past few decades, most researchers have mainly focused on thermodynamic methods. Choi et al. [13] established a 1D steady-state model of a PEM water electrolyzer. The cell voltage was composed of Nernst potential, activation overpotential, ohmic overpotential and concentration overpotential. An equivalent electrical circuit analogy was provided for the sequential kinetic and transport resistances. They found that the reduction kinetics at the cathode were relatively fast, while the anodic overpotential was mainly responsible for the voltage drop. Onda et al. [14,15] conducted a 2D model simulation of a single cell and predicted the stack performance. The convection heat, conduction heat, heat transfer, latent heat of water evaporation, the absorbed heat by the entropy change, and the generated heat by the overpotential for the water electrolyzing reaction were discussed in detail. The shunt resistance was considered in the ohmic loss for the first time. The influence of different temperatures and pressures on the cell performance was compared and measured by experimentation. Gorgun et al. [16] first established a dynamic model under one-dimensional adiabatic conditions. They calculated the partial pressure of the substances at both ends of the anode and cathode through the molar conservation and ideal gas equation and then discussed the phenomena of the electro-osmotic drag force and the diffusion of water in a membrane electrode. The model was established by Simulink, and the polarization curve, efficiency, partial pressure, and current transient change were discussed. Marangio et al. [17] established a one-dimensional non-adiabatic steady-state model. From the perspective of the polarization curve, they studied the activation overpotential, concentration overpotential, and ohmic overpotential of the cathode and anode, respectively. The influence of mass transfer phenomenon on the activation/concentration overpotential and the influence of electrode and flow field plate shunt resistance on the ohmic overpotential were taken into account.
In previous studies, most models have used CFD software with the finite volume method to analyze the flow phenomena in the flow field and porous media. Nie et al. [18] established a 3D dynamic model of the flow field of a PEM water electrolyzer. Through the fluid control equation, the generation rate and the distribution of hydrogen/oxygen in the flow field were calculated by FLUENT at a certain water flow rate. Lobato et al. [19] studied the influence of different flow fields and established a 3D model, including three kinds of flow fields. The pressure drop, velocity distribution, and oxygen mole fraction distribution of the different flow fields were calculated.
Recent years, in order to establish a more accurate model, more researchers have tried to solve the questions of multi-physical field coupling. Kaya et al. [20] developed a 2D isothermal model to investigate the performance of a PEM water electrolyzer, and the effects of Pt and Pt–Ir anode catalysts were compared at the same conditions under different temperatures, membrane thicknesses, current collector lengths, and molar fraction distributions in the channel. Sezgin et al. [21] developed a 3D steady-state isothermal model to observe the effects of various parameters, inlet velocities of the gases, and conductivities of the polymer membrane, and the simulation results of the concentration, current density, and polarization curve were compared with experimental results. Haghayegh et al. [22] established a 3D isothermal steady-state model, coupled the mass transfer phenomenon in the flow field and porous media through the fluid control equation, and calculated the reaction kinetics in a membrane electrode, and they discussed the concentration distribution, pressure distribution, and the influence of temperature and pressure on the current density. Zhang et al. [23] established a 3D steady-state model. The mass transfer, heat transfer, and current distribution were coupled to figure out the influence of co-flow and counter-flow regimes and the influence of different widths of the flow field on the cell performance.
In this study, a steady state, three-dimensional model of a PEM water electrolyzer was established by coupling the electrochemical model, the mass conservation model and the hydrodynamics model. Hence, the current distribution, pressure distribution, transport of species and fluid dynamics were coupled by the electrochemical and mass transfer physical field. The performance of the PEM water electrolyzer under different working conditions were discussed, including the current density, temperature, and pressure. Different patterns of the flow field were simulated, and the performances of the cells using these flow fields were compared. The results can help in optimizing the design of PEM water electrolyzers.

2. Materials and Methods

The structure of a typical PEM water electrolyzer is shown in Figure 1. The polymer electrolyte membrane is in the middle of cell, which plays the role of conducting protons, and it can also isolate the connection between the reactants and products on both sides, ensuring the purity of the products. The catalyst layer is the reaction site in the PEM water electrolyzer, at which water is decomposed into hydrogen and oxygen. The catalyst layer is placed at both ends of the proton exchange membrane by ultrasonic spraying or other technologies, which are the anode catalyst layer and cathode catalyst layer, respectively. This combined structure is called the catalyst coating membrane (CCM). In order to transfer the reactants and products, a layer with a porous structure, entitled the diffusion layer, needs to be added outside the anode and cathode, respectively. The combined structure of the above CCM and the diffusion layer is collectively referred to as the membrane electrode assembly (MEA). Next to the MEA, from the inside to the outside, there are flow field plates, electrified plates, and end plates. The flow field plate is mainly used to control the internal flow state, and the electrified plate is used to connect with an external power supply to apply current to the PEM water electrolyzer. The outermost end plate plays a role of fixing the components.
The MEA is the core component of a PEM water electrolyzer. In addition to the role of the diffusion layer of transferring reactants, the CCM plays a catalytic role in the water electrolysis. The concept of a three-phase boundary (TPB) exists in the water electrolysis reaction [24,25,26]. It is believed that in the water electrolysis reaction, the oxygen evolution reaction (OER) and hydrogen evolution reaction (HER) only occur in a limited-space area, which is at the junction of the PEM, the catalyst, and the pores of the diffusion layer. Among them, the reaction kinetics of the OER have a great impact on the performance of water electrolysis.
The reaction in a PEM water electrolyzer can be expressed as follows:
Anode :   H 2 O 2 H + + 2 e + 1 2 O 2
Cathode :   2 H + + 2 e H 2
Total :   H 2 O 1 2 O 2 + H 2
A flow field plate is also very important for the water electrolysis reaction [27]. In the process of water electrolysis, the reactant (water) and the product (hydrogen) exist in the anode simultaneously. The flow characteristics, diffusion phenomena, mass transfer, and heat transfer characteristics of the mixture play a decisive role in the concentration polarization, affecting the performance of the PEM water electrolyzer. At the same time, the flow field plate and the electrified plate also act the part of conducting electricity, ensuring the continuous progress of the electrochemical reaction. The grooves in the flow field plate are called flow channels. The main types of flow fields are shown in Figure 2, including serpentine flow fields, parallel flow fields, multi-channel serpentine flow fields, interdigital flow fields, grid flow fields, spiral flow fields, etc. The reactant (water) produces an oxygen and hydrogen flow in the diffusion layer and the flow field. According to the principle of PEM water electrolyzers, the supply and flow of water will affect its electrolysis efficiency after reaching the CCM, and different types of CCM will also affect the electrolysis efficiency, affecting the degree of electrochemical polarization. At the same time, the amount of products obtained by electrolysis at the cathode and anode and the difficulty of discharging the products from the diffusion layer and the flow field will cause a change in the substance concentration in the PEM water electrolyzer, which will have a certain impact on the concentration polarization. Under the combined action of electrochemistry and concentration polarization, PEM water electrolyzers with different flow fields and MEAs will have different performances, so it is necessary to study the types and properties of different structures.

2.1. Electrochemical Analysis

The PEM water electrolysis reaction is affected by the electrochemical reaction kinetics and the transport of protons and electrons. Among them, the mass transfer resistance in the flow field and the diffusion layer of the PEM hydrolysis cell have a very important influence on the above physical and chemical processes. The mass transfer resistance is mainly divided into the diffusion resistance of the reactant products, electro-osmotic drag resistance, and cross resistance [21,23,28]. Diffusion resistance is the overvoltage caused by the change in the equivalent resistance caused by the change in the concentration of the cathode and anode, and electro-osmotic drag resistance is the resistance caused by protons passing through the membrane, which is mainly related to the humidity of the PEM. Cross resistance is caused by the gas passing through the PEM, which has little impact. The mass transfer resistance will greatly affect the electrochemical reaction process and produce a certain degree of overpotential, affecting the performance of the PEM water electrolyzer.
By establishing a numerical model of the electrochemical reactions and mass transfer processes in a PEM water electrolyzer, the performance of the PEM water electrolyzer and the distribution of the physical quantities such as speed, pressure, material concentration, and current density distribution inside the PEM water electrolyzer can be simulated, which can give a more intuitive and clear understanding of the internal situation of the water electrolytic cell and also ensure the correctness of the performance comparison under different working conditions [29]. In this section, some assumptions about the mathematical model in this study are made, and then a complete water electrolysis model is established in three parts, namely an electrochemical model, a mass conservation model, and a momentum conservation model. Finally, the boundary conditions are set.
An electrochemical model was developed for a PEM water electrolyzer. A PEM water electrolyzer needs to consume external electric energy to split water. The polarization curve, which is the relationship between the current density and cell voltage, depicts the performance of a water electrolyzer. The cell voltage mainly consist of four parts, which can be expressed as follow [30]:
E c e l l = E O C V + η a c t + η o h m + η d i f f
where E O C V is the open circuit voltage, η a c t is the activation overpotential, η o h m is the ohmic overpotential, and η d i f f is the diffusion overpotential. The theoretical energy is the required energy for water splitting without any losses, which is represented by the open circuit voltage. The open circuit voltage is the reversible potential, which can be determined by the Nernst equation and can be calculated as follows [31]:
E O C V = Δ G n F
where Δ G is the Gibbs free energy, n is the number of electrons transferred in the semi-reaction, and F is the Faraday constant.
During the electrochemical reaction in water splitting, the Gibbs free energy is affected by the temperature and pressure conditions. The difference in the behavior between a real gas and an ideal gas is dependent only on the pressure and the temperature and not on the presence of any other gases. At a given temperature, the “effective” pressure of a gas is given by its fugacity. By historical convention, fugacities have the dimension of pressure, so the dimensionless activity is given by:
α i = f i p * = φ i y i p p *
where φ i is the dimensionless fugacity coefficient of the species, y i is its mole fraction in the gaseous mixture ( y = 1 for a pure gas), and p is the total pressure. The value p * is the standard pressure; it may be equal to 1 atm (101.325 kPa) or 1 bar (100 kPa). The Gibbs free energy at a constant temperature and pressure can be calculated by the Nernst equation:
Δ G = Δ G * + R T ln p H 2 p O 2 1 / 2 α H 2 O
where Δ G * is the reference Gibbs free energy at standard pressure and different cell temperatures. The species enthalpies and entropies are used to calculate the reference Gibbs free energy as follows [32]:
Δ G * = Δ H * T Δ S * = ν i H i T T ν i S i T , P
where ν i is the stoichiometric coefficients.
The steady form of the charge balance can be expressed as:
j l + j s = 0
where j l is the ion current in the membrane, and j s is the electron current in the electrodes, both of which can be calculated by Ohm’s law:
j l = σ l eff ϕ l
j s = σ s eff ϕ s
where is the Laplace operator, σ l eff is the effective ion conductivity, and σ s eff is the effective conductivity of the electrodes, determined by the Bruggeman correction:
σ l eff = σ l ε l 3 / 2
σ s eff = σ s ε s 3 / 2
where σ l is the membrane ion conductivity coefficient, and ε l is the porosity. This can be calculated as follows [33]:
σ l = 0.005193 λ 0.00326 exp 1268 1 303 1 T
where λ is the humidification degree.
The activation and diffusion overpotentials at the anode and cathode catalyst layer depend on the electrode potential, electrolyte potential, and open circuit voltage, which can be expressed as follows:
η a c t + η d i f f = ϕ s ϕ l E O C V = E c e l l η o h m E O C V
The activation overpotential is the energy loss caused by the electrochemical reaction, which can be calculated by the Butler-Volmer equation as follows [34]:
j v = a v j 0 exp α a F R T η a c t exp α c F R T η a c t
where a v is the electrode specific surface area, α a and α c are the charge transfer coefficients on the anode and cathode, respectively, R is the gas constant, and T is the cell temperature. j 0 is the exchange current density, which influences the diffusion overpotential, and it can be calculated as follows [35]:
j 0 = j 0 , r e f i c i c i , r e f α a v i n i c i c i , r e f α c v i n
where the c i is the species concentration, and c i , r e f is the reference concentration.

2.2. Momentum Conservation

The momentum equation in a porous domain is described by the Brinkman equation, in which the gas velocity is approximated by Darcy’s law and the continuity equation, which can be calculated as follows [36]:
ρ ε u u ε = p + μ ε u 2 μ 3 ε u μ k + Q u
ρ u = Q
where ρ is the density of the gas mixture, ε is the electrode porosity, u is the mass velocity, μ is the viscosity, k is the electrode permeability, and Q represents the mass source, which can be calculated as:
Q = m i R i , m M i
The Navier–Stokes and continuity equations without mass creation depict the flow in the gas channels, which can be expressed as follows [37]:
Ρ u u = p + μ u 2 3 μ u
ρ u = 0

2.3. Mass Conservation

The flow distribution in the gas channel can be described by the Maxwell–Stefan equation because of the multicomponent diffusion and convection conditions [20]:
t ρ ω i + ρ ω i j = 1 N D i j M M j ω j + ω j M M + x j ω j P P + ω i ρ u = q i .  
where ω i is the mass fraction, P is the partial pressure, M is the molar mass, and q i is the species flux, which is a function of the local current density.
q H 2 = j c 2 F M H 2
q O 2 = j a 4 F M O 2
q H 2 O = j a 2 F M H 2 O
The diffusivity in the free space can be expressed as follows [38]:
D i j = a T T c i T c j b P c i P c j 1 / 3 T c i T c j 5 / 12 1 M i + 1 M j 1 / 2 P
The effective diffusion coefficient in the electrodes, which is approximated by the Bruggeman correction, can be expressed as:
D i , j eff = D i , j ε 1.5
The density of the gas mixture can be calculated as follows:
ρ = p + p ref M RT
M = i w i M i 1
x k is the molar fraction, which can be calculated as:
x k = w k M k M n
In this research, the water inlet was related to the stoichiometry, and the mass flow rate can be expressed as follow:
M f l u x i n = λ M H 2 O j A 2 F
The electrochemical model was established by using the Nernst equation, the Butler–Volmer formula, Ohm’s law, and some semi-empirical formulas, and the mass conservation model was established by using the mass conservation equation, the Stefan–Maxwell formula, and some basic thermodynamic formulas. The hydrodynamic model was established by using the kinetic conservation and the Brinkman equation. By coupling the above three models, a complete mathematical model of a PEM water electrolyzer was established. The parameters and boundary conditions of the model were set on the basis of actual working conditions.

2.4. Physical modeling

The geometrical model used in this paper is shown in Figure 2. The structure from top to bottom was mainly composed of four parts: the anode flow field, which was composed of 1 mm × 1 mm microchannels; the anode electrode, including the anode gas diffusion layer and the anode catalyst layer; the PEM; and the cathode electrode, including the cathode gas diffusion layer and the cathode catalyst layer.
In addition to the parallel flow field, the other three common flow fields were also numerically investigated in this paper. The different flow fields are as show in Figure 3.
By using COMSOL Multiphysics, the three-dimensional geometric model of a PEM water electrolyzer was established. FEM was used to mesh the geometric model and verify the grid quality. The mathematical model was substituted by region. The parameters were changed according to different experimental conditions, and the accuracy of the model was verified. The results showed that the mathematical model was in good agreement with the experimental results and had a high accuracy.
Geometric models of PEM water electrolyzers with different flow patterns were established and substituted into the same mathematical model.

3. Results and Discussions

3.1. Model Validation

The mathematical model in this study was established and calculated by using COMSOL Multiphysics. To validate the present model, the numerical results were compared with the previous experimental results reported by Hansen et al. [39] at the same working conditions and with a similar geometric form. The working conditions included a cell temperature of 403   K and a pressure of 1   atm , and the geometric parameters and other working conditions are shown in Table 1.
Figure 4 illustrates the comparison of the present results with the experimental data, and the results showed a good agreement with the experimental data. The coefficient of determination of fitting was 0.97. After the validation, the effective model was used to calculate and analyze the impact of different working conditions and geometric forms on the cell performance.

3.2. Effect of Different Current Density on Mass Transport

The hydrogen and oxygen production rates were closely related to the current density. According to Faraday’s law, the molar production rate n ˙ i corresponded to the current density, which can be expressed as follows [30]:
n ˙ i = j A n F
The mole concentration was proportional to the molar production rate, which can be calculated by Equation (32). Figure 5 illustrates the oxygen mole concentration distribution at the anode at different current densities. It can be seen that when the current density was increased, the oxygen mole concentration increased in the anode, while most of the increase was at the edge of electrode. It is worth mentioning that when the current density was too high, the oxygen generated at the electrode was hard to remove, which slowed the continuous progress of the reaction.
According to Equation (30), the current density was related to the inlet mass flow rate, which influenced the velocity in the flow field. Figure 6 illustrates the velocity field at the same current density as above. At the inlet of the flow field, the velocity at 1.5   A   cm 2 was 2.1   m   s 1 , which was three times larger than the velocity at 0.5   A   cm 2 and twice the velocity at 1.0   A   cm 2 . This can also explain the increase in the molar concentration.

3.3. Effect of Cell Temperature

From the present model, it was found that the open circuit voltage, activation overpotential, mass concentration, diffusion coefficient, and membrane conductivity were closely related to the cell temperature. Figure 7 illustrates the oxygen mole concentration distribution at different temperatures. When the temperature increased, the molar concentration dropped evidently at the edge of electrode, and concentration drop was small at the boundary of the electrode and the flow field.

3.4. Effect of Different Flow Fields

Figure 8 shows the velocity distribution in different flow fields when the current density was 1.0   A   cm 2 . It can be seen that at the channel turning point in the flow field, the velocity decreased along the channel, which was caused by the higher pressure difference between the adjacent channels compared to other areas. Comparing the velocity distribution of each flow field, the maximum velocity in the parallel flow field was 1.32   m   s 1 , the maximum velocity in the single-channel serpentine flow field was 1.62   m   s 1 , the maximum velocity in the multi-channel serpentine flow field was 1.56   m   s 1 , and the maximum velocity in the interdigital flow field was 1.36   m   s 1 . The velocity distribution in the two serpentine flow fields was relatively uniform. The peak velocity of the parallel flow field and interdigital flow field was mainly at the inlet and outlet, and the velocity in the flow channel was relatively small, which led to a slower migration of the reactants through the porous media. In this way, although the reaction could be carried out more fully, the products were discharged slowly.
The anode pressure distribution of hydro digesters in different flow fields under a current density of 1.0   A   cm 2 is shown in Figure 9. The pressure distribution map shows which flow field had a low pressure drop. A PEM water electrolyzer with a low-pressure-drop flow field can consume less energy when distributing the reactants, thus improving the performance to a certain extent. In the figure, the pressure drop of the parallel flow field was 3.46   Pa , and the pressure dropped uniformly along the axis connecting the inlet and outlet. The pressure drop of the single serpentine flow field was 38.9   Pa , and that of the multi serpentine flow field was 11.6   Pa . The pressure of both serpentine flow fields decreased gradually along the direction of the flow passage. The pressure drop of the interdigital flow field was 7.08   Pa , and the pressure of the inlet side and outlet side of the flow channel decreased slowly along their respective flow channels, while the pressure of the flow channel from the inlet to the outlet dropped rapidly. It can be clearly seen that the pressure drop of the parallel flow field was lower than that of the other flow fields. By comparing the two serpentine flow fields, it can be seen that there were fewer bends in the channel of the parallel flow field. Based on the analysis of the velocity distribution in the previous section, the flow velocity in the parallel flow field was also relatively slow, so the flow was less obstructed. Compared with the interdigital flow field, the main reason for this is that the mixed gas can be directly discharged to the outlet, while the interdigital flow field must pass through the porous media, which increases the pressure drop.
The oxygen concentration distribution is the result of the interaction between the oxygen evolution rate and the mixture flow. On the one hand, a higher oxygen concentration can reflect a faster overall reaction rate, and the same current density indicates more reaction points, and the electrochemical polarization is weak, so it has a better electrolytic performance. On the other hand, this may indicate that the slower flow rate affects the discharge of oxygen, and a stronger concentration polarization, a worse electrolytic performance, and an uneven distribution of the oxygen region also verifies this point. By comparing the oxygen concentration with the current density distribution, the point which dominates the reaction can be obtained.
Figure 10 shows the anodic oxygen concentration distribution in anode electrode of the PEM water electrolyzer with different flow fields when the current density was 1.0   A   cm 2 . The oxygen concentration in each flow field increased from the inlet to the outlet. Among them, the maximum oxygen concentration of the anode electrode of the PEMEC with a parallel flow field reached 13.7   mol   m 3 , which was located in the middle of the side near the outlet, indicating that the water entered the porous medium for a full reaction and had the highest oxygen generation rate on the whole. The maximum oxygen concentration of the interdigital flow field reached 13.5   mol   m 3 , which was mainly located at the end of interdigital flow field on the inlet side and at the edge of electrode. At the same time, it can be seen that the parallel and interdigital flow fields had a high oxygen concentration due to the slow flow rate and the accumulation of products. The maximum oxygen concentration of the two serpentine flow fields reached 11.0   mol   m 3 and increased uniformly from the inlet to the outlet.
As the reactants entered the flow field and penetrated into the electrode, they were continuously consumed, resulting in a concentration difference. In the flow process, pressure gradients appeared in the middle and edge of the flow channel. Figure 11 shows the corresponding current density distribution, with an average current density of 1.0   A   cm 2 . It can be seen from the figure that the water electrolysis reaction was mainly distributed at the junction of the flow field and the electrode. The maximum current density of the anode in the parallel flow field was 1.87   A   cm 2 , and the maximum current density was mainly at the junction of the inlet and outlet section and the flow channel. The maximum current density of the single-channel serpentine flow field was 1.76   A   cm 2 , and that of the multi-channel serpentine flow field was 1.8   A   cm 2 . The maximum current density was mainly at the position where the flow channel bends. The current density distribution of the above three flow fields was relatively uniform. The highest current density of the interdigital flow field anode was 2.1   A   cm 2 . The highest current density was mainly at the junction of the inlet and outlet sections and the flow channel. At the same time, the reaction points were concentrated at the junction of the flow field and the electrode, which seriously affected the reaction efficiency. In combination with the oxygen concentration distribution, it can be considered that the high oxygen concentration position of the interdigital flow field was due to the accumulation of oxygen that was not easily discharged. To sum up, the reasons for the performance differences between the different flow fields could be analyzed from the perspective of the oxygen concentration distribution and current density distribution.

4. Conclusions

In this research, a steady state, three-dimensional mathematical model coupled with electrochemical and mass transfer physical fields for a PEM water electrolyzer was established and validated. Under the same working condition, the flow velocity, oxygen and water concentration distribution, current density distribution, and polarization curve of the different models were compared. The results showed that the single-channel and multi-channel serpentine flow fields had faster flow velocities than the parallel and interdigital flow fields. Although the distribution of the reactants and products was more uniform, they caused a greater pressure drop, so the diffusion overpotential increased. The oxygen concentration distribution of the two serpentine flow fields was similar, and the parallel flow field and interdigital flow field had slightly higher values. However, the interdigital flow field was mainly caused by the difficulty of the mixture entering and exiting the porous media, the low oxygen generation rate, and the difficulty in discharging. Taking the polarization curve as the main factor and considering many aspects comprehensively, the parallel flow field had a good performance in terms of electrolyzing water.

Author Contributions

Conceptualization, J.M. and J.Z.; methodology, J.Z.; software, J.Z.; validation, J.Z.; formal analysis, J.Z. and B.H.; investigation, J.Z. and J.M.; resources, J.M.; data curation, J.Z.; writing—original draft preparation, J.Z.; writing—review and editing, J.M., J.Z., B.H. and Z.K.; visualization, J.Z.; supervision, J.M.; project administration, J.M.; funding acquisition, J.M. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by National Natural Science Foundation of China (11902080) and the start-up fund from Fudan University (JIH2126004Y).

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Not applicable.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Abe, J.O.; Popoola, A.P.I.; Ajenifuja, E.; Popoola, O.M. Hydrogen energy, economy and storage: Review and recommendation. Int. J. Hydrogen Energy 2019, 44, 15072–15086. [Google Scholar] [CrossRef]
  2. Turner, J.A. A Realizable Renewable Energy Future. Science 1999, 285, 687–689. [Google Scholar] [CrossRef] [PubMed]
  3. Edwards, P.P.; Kuznetsov, V.L.; David, W.I.F.; Brandon, N.P. Hydrogen and fuel cells: Towards a sustainable energy future. Energy Policy 2008, 36, 4356–4362. [Google Scholar] [CrossRef]
  4. Turner, J.A. Sustainable hydrogen production. Science 2004, 305, 972–974. [Google Scholar] [CrossRef]
  5. Carmo, M.; Fritz, D.L.; Mergel, J.; Stolten, D. A comprehensive review on PEM water electrolysis. Int. J. Hydrogen Energy 2013, 38, 4901–4934. [Google Scholar] [CrossRef]
  6. Barbir, F. PEM electrolysis for production of hydrogen from renewable energy sources. Sol. Energy 2005, 78, 661–669. [Google Scholar] [CrossRef]
  7. Wang, Y.; Chen, K.S.; Mishler, J.; Cho, S.C.; Adroher, X.C. A review of polymer electrolyte membrane fuel cells: Technology, applications, and needs on fundamental research. Appl. Energy 2011, 88, 981–1007. [Google Scholar] [CrossRef]
  8. Chi, J.; Yu, H. Water electrolysis based on renewable energy for hydrogen production. Chin. J. Catal. 2018, 39, 390–394. [Google Scholar] [CrossRef]
  9. Kumar, S.S.; Himabindu, V. Hydrogen production by PEM water electrolysis—A review. Mater. Sci. Energy Technol. 2019, 2, 442–454. [Google Scholar] [CrossRef]
  10. Din, M.A.U.; Irfan, S.; Sharif, H.M.A.; Jamil, S.; Idrees, M.; Khan, Q.U.; Nazir, G.; Cheng, N. Nanoporous CoCu layered double hydroxide nanoplates as bifunctional electrocatalyst for high perfor-mance overall water splitting. Int. J. Hydrogen Energy 2022, 48, 5755–5763. [Google Scholar] [CrossRef]
  11. Din, M.A.U.; Idrees, M.; Jamil, S.; Irfan, S.; Nazir, G.; Mudassir, M.A.; Saleem, M.S.; Batool, S.; Cheng, N.; Saidur, R. Advances and challenges of methanol-tolerant oxygen reduction reaction electrocatalysts for the direct methanol fuel cell. J. Energy Chem. 2022, 77, 499–513. [Google Scholar] [CrossRef]
  12. Rashid, M.; Al Mesfer, M.K.; Naseem, H.; Danish, M. Hydrogen production by water electrolysis: A review of alkaline water electrolysis, PEM water electrolysis and high temperature water electrolysis. Int. J. Eng. Adv. Technol. 2015, 4, 3. [Google Scholar]
  13. Choi, P.; Bessarabov, D.G.; Datta, R. A simple model for solid polymer electrolyte (SPE) water electrolysis. Solid State Ion. 2004, 175, 535–539. [Google Scholar] [CrossRef]
  14. Onda, K.; Murakami, T.; Hikosaka, T.; Kobayashi, M.; Notu, R.; Ito, K. Performance Analysis of Polymer-Electrolyte Water Electrolysis Cell at a Small-Unit Test Cell and Performance Prediction of Large Stacked Cell. J. Electrochem. Soc. 2002, 149, A1069–A1078. [Google Scholar] [CrossRef]
  15. Onda, K.; Kyakuno, T.; Hattori, K.; Ito, K. Prediction of production power for high-pressure hydrogen by high-pressure water electrolysis. J. Power Sources 2004, 132, 64–70. [Google Scholar] [CrossRef]
  16. Görgün, H. Dynamic modelling of a proton exchange membrane (PEM) electrolyzer. Int. J. Hydrogen Energy 2006, 31, 29–38. [Google Scholar] [CrossRef]
  17. Marangio, F.; Santarelli, M.; Calì, M. Theoretical model and experimental analysis of a high pressure PEM water electrolyser for hydrogen production. Int. J. Hydrogen Energy 2009, 34, 1143–1158. [Google Scholar] [CrossRef]
  18. Nie, J.; Chen, Y. Numerical modeling of three-dimensional two-phase gas–liquid flow in the flow field plate of a PEM electrolysis cell. Int. J. Hydrogen Energy 2010, 35, 3183–3197. [Google Scholar] [CrossRef]
  19. Lobato, J.; Cañizares, P.; Rodrigo, M.A.; Pinar, F.J.; Mena, E.; Úbeda, D. Three-dimensional model of a 50 cm2 high temperature PEM fuel cell. Study of the flow channel geometry influence. Int. J. Hydrogen Energy 2010, 35, 5510–5520. [Google Scholar] [CrossRef]
  20. Kaya, M.F.; Demir, N. Numerical Investigation of PEM Water Electrolysis Performance for Different Oxygen Evolution Electrocatalysts. Fuel Cells 2017, 17, 37–47. [Google Scholar] [CrossRef]
  21. Sezgin, B.; Caglayan, D.G.; Devrim, Y.; Steenberg, T.; Eroglu, I. Modeling and sensitivity analysis of high temperature PEM fuel cells by using Comsol Multiphysics. Int. J. Hydrogen Energy 2016, 41, 10001–10009. [Google Scholar] [CrossRef]
  22. Haghayegh, M.; Eikani, M.H.; Rowshanzamir, S. Modeling and simulation of a proton exchange membrane fuel cell using computational fluid dynamics. Int. J. Hydrogen Energy 2017, 42, 21944–21954. [Google Scholar] [CrossRef]
  23. Zhang, Z.; Xing, X. Simulation and experiment of heat and mass transfer in a proton exchange membrane electrolysis cell. Int. J. Hydrogen Energy 2020, 45, 20184–20193. [Google Scholar] [CrossRef]
  24. Kong, W.; Zhang, M.; Han, Z.; Zhang, Q. A Theoretical Model for the Triple Phase Boundary of Solid Oxide Fuel Cell Electrospun Electrodes. Appl. Sci. 2019, 9, 493. [Google Scholar] [CrossRef]
  25. Vijay, P.; Tadé, M.O.; Shao, Z.; Ni, M. Modelling the triple phase boundary length in infiltrated SOFC electrodes. Int. J. Hydrogen Energy 2017, 42, 28836–28851. [Google Scholar] [CrossRef]
  26. O’Hayre, R.; Prinz, F.B. The Air/Platinum/Nafion Triple-Phase Boundary: Characteristics, Scaling, and Implications for Fuel Cells. J. Electrochem. Soc. 2004, 151, A756–A762. [Google Scholar] [CrossRef]
  27. Arvay, A.; French, J.; Wang, J.-C.; Peng, X.-H.; Kannan, A. Nature inspired flow field designs for proton exchange membrane fuel cell. Int. J. Hydrogen Energy 2013, 38, 3717–3726. [Google Scholar] [CrossRef]
  28. Lee, C.H.; Banerjee, R.; Arbabi, F.; Hinebaugh, J.; Bazylak, A. Porous Transport Layer Related Mass Transport Losses in Polymer Electrolyte Membrane Electrolysis—A Review. In Proceedings of the ASME 14th International Conference on Nanochannels, Microchannels, and Minichannels, Washington, DC, USA, 10–14 July 2016. [Google Scholar] [CrossRef]
  29. Mo, J.; Kang, Z.; Retterer, S.T.; Cullen, D.A.; Toops, T.J.; Green, J.B.; Mench, M.M.; Zhang, F.-Y. Discovery of true electrochemical reactions for ultrahigh catalyst mass activity in water splitting. Sci. Adv. 2016, 2, e1600690. [Google Scholar] [CrossRef] [PubMed]
  30. Mench, M.M. Fuel Cell Engines; John Wiley & Sons: Hoboken, NJ, USA, 2008. [Google Scholar]
  31. Pilatowsky, I.; Romero, R.; Isaza, C.; Gamboa, S.; Sebastian, P.; Rivera, W. Thermodynamics of fuel cells. In Cogeneration Fuel Cell-Sorption Air Conditioning Systems; Springer: Cham, Switzerland, 2011; pp. 25–36. [Google Scholar] [CrossRef]
  32. Gibbs, J. A Method of Geometrical Representation of the Thermodynamic Properties by Means of Surfaces; Yale University Press: New Haven, CN, USA, 1957; pp. 33–54. [Google Scholar]
  33. Kang, Z.; Mo, J.; Yang, G.; Li, Y.; Talley, D.A.; Han, B.; Zhang, F.-Y. Performance Modeling and Current Mapping of Proton Exchange Membrane Electrolyzer Cells with Novel Thin/Tunable Liquid/Gas Diffusion Layers. Electrochim. Acta 2017, 255, 405–416. [Google Scholar] [CrossRef]
  34. Dickinson, E.J.; Wain, A.J. The Butler-Volmer equation in electrochemical theory: Origins, value, and practical application. J. Electroanal. Chem. 2020, 872, 114145. [Google Scholar] [CrossRef]
  35. Plieth, W. Electrochemistry for Materials Science; Elsevier: Amsterdam, The Netherlands, 2008. [Google Scholar] [CrossRef]
  36. Whitaker, S. Flow in porous media I: A theoretical derivation of Darcy’s law. Transp. Porous Media 1986, 1, 3–25. [Google Scholar] [CrossRef]
  37. Raja, L.L.; Kee, R.J.; Deutschmann, O.; Warnatz, J.; Schmidt, L.D. A critical evaluation of Navier–Stokes, boundary-layer, and plug-flow models of the flow and chemistry in a catalytic-combustion monolith. Catal. Today 2000, 59, 47–60. [Google Scholar] [CrossRef]
  38. Han, B.; Steen, S.M., III; Mo, J.; Zhang, F.Y. Electrochemical performance modeling of a proton exchange membrane electrolyzer cell for hydrogen energy. Int. J. Hydrogen Energy 2015, 40, 7006–7016. [Google Scholar] [CrossRef]
  39. Hansen, M.K.; Aili, D.; Christensen, E.; Pan, C.; Eriksen, S.; Jensen, J.O.; von Barner, J.H.; Li, Q.; Bjerrum, N.J. PEM steam electrolysis at 130 °C using a phosphoric acid doped short side chain PFSA membrane. Int. J. Hydrogen Energy 2012, 37, 10992–11000. [Google Scholar] [CrossRef]
Figure 1. Schematic diagram of a standard PEM water electrolyzer.
Figure 1. Schematic diagram of a standard PEM water electrolyzer.
Materials 16 01310 g001
Figure 2. The computational domain of the model.
Figure 2. The computational domain of the model.
Materials 16 01310 g002
Figure 3. The different flow fields of PEMEC. (A) Parallel, (B) serpentine, (C) multi-serpentine, and (D) interdigitate.
Figure 3. The different flow fields of PEMEC. (A) Parallel, (B) serpentine, (C) multi-serpentine, and (D) interdigitate.
Materials 16 01310 g003
Figure 4. Comparison and validation of the present model with experimental data.
Figure 4. Comparison and validation of the present model with experimental data.
Materials 16 01310 g004
Figure 5. Oxygen mole fraction distribution in the anode under different current densities: (A)   0.5   A   cm 2 , (B) 1.0   A   cm 2 , and (C) 1.5   A   cm 2 .
Figure 5. Oxygen mole fraction distribution in the anode under different current densities: (A)   0.5   A   cm 2 , (B) 1.0   A   cm 2 , and (C) 1.5   A   cm 2 .
Materials 16 01310 g005
Figure 6. Velocity distribution in the flow field under different current densities: (A)   0.5   A   cm 2 , (B) 1.0   A   cm 2 , and (C) 1.5   A   cm 2 .
Figure 6. Velocity distribution in the flow field under different current densities: (A)   0.5   A   cm 2 , (B) 1.0   A   cm 2 , and (C) 1.5   A   cm 2 .
Materials 16 01310 g006
Figure 7. Oxygen mole fraction distribution under different temperatures: (A)   40   ° C , (B)   80   ° C , (C)   120   ° C , and (D)   160   ° C .
Figure 7. Oxygen mole fraction distribution under different temperatures: (A)   40   ° C , (B)   80   ° C , (C)   120   ° C , and (D)   160   ° C .
Materials 16 01310 g007
Figure 8. Velocity distribution in microchannels with different flow fields. (A) Parallel, (B) serpentine, (C) multi-serpentine, and (D) interdigitate.
Figure 8. Velocity distribution in microchannels with different flow fields. (A) Parallel, (B) serpentine, (C) multi-serpentine, and (D) interdigitate.
Materials 16 01310 g008
Figure 9. Pressure distribution in microchannels with different flow fields. (A) Parallel, (B) serpentine, (C) multi-serpentine, and (D) interdigitate.
Figure 9. Pressure distribution in microchannels with different flow fields. (A) Parallel, (B) serpentine, (C) multi-serpentine, and (D) interdigitate.
Materials 16 01310 g009
Figure 10. Oxygen mole concentration in anode electrode with different flow fields. (A) Parallel, (B) serpentine, (C) multi-serpentine, and (D) interdigitate.
Figure 10. Oxygen mole concentration in anode electrode with different flow fields. (A) Parallel, (B) serpentine, (C) multi-serpentine, and (D) interdigitate.
Materials 16 01310 g010
Figure 11. The current density distribution on anode electrode under different flow fields. (A) Parallel, (B) serpentine, (C) multi-serpentine, (D) interdigitate.
Figure 11. The current density distribution on anode electrode under different flow fields. (A) Parallel, (B) serpentine, (C) multi-serpentine, (D) interdigitate.
Materials 16 01310 g011
Table 1. Basic parameters used in the PEMEC modeling.
Table 1. Basic parameters used in the PEMEC modeling.
Description, SymbolValue, Unit
MEA active area, A 4   cm 2
Operating pressure, P 1 ,   2 ,   4   a t m
Operating temperature, T 40 ,   80 ,   120 ,   160   ° C
Electrode porosity, ε 0.4
Electrode specific surface area, a v 10 9   cm 2
Anode charge transfer coefficients, α a 2
Cathode charge transfer coefficients, α c 0.5
Anode exchange current density, α 0 a 6 × 10 10   A   cm 2
Cathode exchange current density, α 0 a 0.34   A   cm 2
PEM N a f i o n @   115
Electrode thickness, δ 200   mm
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Share and Cite

MDPI and ACS Style

Zheng, J.; Kang, Z.; Han, B.; Mo, J. Three-Dimensional Numerical Simulation of the Performance and Transport Phenomena of Oxygen Evolution Reactions in a Proton Exchange Membrane Water Electrolyzer. Materials 2023, 16, 1310. https://doi.org/10.3390/ma16031310

AMA Style

Zheng J, Kang Z, Han B, Mo J. Three-Dimensional Numerical Simulation of the Performance and Transport Phenomena of Oxygen Evolution Reactions in a Proton Exchange Membrane Water Electrolyzer. Materials. 2023; 16(3):1310. https://doi.org/10.3390/ma16031310

Chicago/Turabian Style

Zheng, Jinsong, Zhenye Kang, Bo Han, and Jingke Mo. 2023. "Three-Dimensional Numerical Simulation of the Performance and Transport Phenomena of Oxygen Evolution Reactions in a Proton Exchange Membrane Water Electrolyzer" Materials 16, no. 3: 1310. https://doi.org/10.3390/ma16031310

APA Style

Zheng, J., Kang, Z., Han, B., & Mo, J. (2023). Three-Dimensional Numerical Simulation of the Performance and Transport Phenomena of Oxygen Evolution Reactions in a Proton Exchange Membrane Water Electrolyzer. Materials, 16(3), 1310. https://doi.org/10.3390/ma16031310

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop