Development of a novel quasi-2D PEM Electrolyzer Model in Modelica

To increase the efficiency of PEM electrolysis, simulation models are required that accurately describe the system's electrochemical and thermal behavior in a computationally efficient manner and are thus suitable for developing control strategies. Therefore, a quasi-2D PEM electrolyzer model is presented in this paper, which is a compromise between the previously developed models regarding their model complexity. The electrochemical behavior is described with equations commonly used in the literature and the thermal behavior with correlations for gas-liquid heat transfer. Preliminary validation indicates that the model can describe the electrochemical behavior and thermal dynamics of a PEM electrolysis stack with good accuracy.


Introduction
Hydrogen will play a decisive role in the decarbonization of future energy systems.Consequently, the number of electrolyzers worldwide will have to increase significantly in the future to be able to produce the required quantities with low emissions.According to the National Hydrogen Strategy of the German Government, a hydrogen demand of approx.90 to 110 TWh a -1 is expected by the year 2030.Generation plants with a total installed capacity of 5 GW are planned to meet this demand (BMWK 2020).The European Commission (2020) even expects a total installed capacity of 40 GW in the EU.
Proton exchange membrane (PEM) electrolysis is of particular importance here because of its suitability for coupling with volatile sources of electrical energy due to its rapid start-up and shutdown behavior and its partial and overload capability.But according to the current state of technology, only stack efficiencies between 56 and 74 % based on the lower heating value of hydrogen are achieved (Tjarks 2017).Approximately one-third of the electrical energy supplied is thus converted into heat.This heat is currently mostly dissipated directly to the environment via heat exchangers.The overall efficiency and costeffectiveness of PEM electrolysis could be significantly increased by using this waste heat.However, to exploit this unused potential, models are needed that accurately describe the heat transfer processes in the electrolyzer.In particular, for the development of control strategies, models are needed that can realistically represent the thermal dynamics of electrolysis stacks.
While their electrochemical behavior has been extensively studied, their dynamic thermal behavior has been mostly simplified.A large number of models use the so-called lumped parameter approach, where incoming and outgoing heat fluxes are calculated based on the assumption that the electrolysis stack has a uniform temperature at each time step (Crespi et al. 2023;Espinosa-López et al. 2018;García-Valverde, Espinosa and Urbina 2012;Sood et al. 2020).This type of model has already been implemented in Modelica by Webster and Bode (2019).They are computationally efficient, but cannot represent the heat transfer within the stack and require intensive experimental studies for parameterization.On the other hand, there are multiple models describing and investigating the heat transfer in the electrolysis cell using complex 3D finite volume approaches (Ma et al. 2021;Toghyani, Afshari and Baniasadi 2019;Zhang and Xing 2020).Although these models represent heat transfer very accurately, they are not suitable for dynamic simulation over longer periods due to their complexity.
In the field of fuel cell modeling, some authors discretize the cell components in only one dimension and then couple the discretized cell components with each other.These approaches are called 1D+1D or quasi-2D models (Gong et al. 2022;Tang et al. 2017).Some similar approaches have also been developed in the field of PEM electrolysis but without heat transfer description in the flow channels (Kim, Park and Lee 2013;Lin and Zausch 2022).Therefore, a quasi-2D model of a PEM electrolyzer is presented in this work, which can describe the heat transfer processes in the individual cells and thus the thermal dynamics of the entire stack.

General Structure
Figure 1 shows the basic structure of a PEM electrolysis cell.In an electrolysis stack, several cells are connected in series, with the bipolar plate of the cathode serving as the bipolar plate of the next cell's anode.In the present model, a parallel flow field design is assumed and the cells are discretized in flow direction only.Temperatures, flow velocities, etc. perpendicular to the channel orientation are thus assumed to be uniform in the respective volumes.Furthermore, the porous transport and catalytic layers as well as the PEM are combined to a uniform thermal mass and together represent the Membrane Electrode Assembly (MEA), where the conversion of water to hydrogen and oxygen takes place.Figure 2 shows the discretized electrolysis cell's structure in the Modelica development environment.The TILSuite package by TLK Thermo GmbH serves as the model basis for the fluid data calculation and the creation of boundary conditions.To be precise, the Connector, Boundary, Splitter, Joiner, and TILMedia substance data models were used.
The anode and cathode channel volumes each have a gas and a water inlet and outlet at the top and bottom, allowing the individual volumes to be interconnected.Except for the anode's water inlet, all inlets of the first cell volumes are provided with boundaries whose mass flow is equal to zero, since process water is usually the only incoming mass flow.The MEA is connected to the flow channel volumes via three connectors to describe the gas, water, and heat exchange between them.The heat ports on the left and right connect the bipolar plate and cathode flow channel volumes to the neighboring cells.The heat ports at the top and bottom connect the bipolar plate to the ambient.Based on the single-cell model, the stack is now discretized in cell direction.For this purpose, the left and right heat ports of the individual cells are connected.The flow inputs and outputs are connected using so-called joiners and splitters (Figure 3).Joiners add the individual cells' mass flows and form the arithmetic average of their temperatures.Splitters distribute the incoming mass flow evenly over the cells.The heat ports of the first and last cells are connected to the end plates.These have a significantly larger volume than the bipolar plates and are therefore considered separately.All remaining heat ports are connected to heat boundaries representing the ambient temperature.Due to the asymmetric structure of the single cell, a discretization according to Figure 3 would model a redundant bipolar plate.This was solved in the original model by adding a single cell without a bipolar plate but is not shown in Figure 3 for reasons of clarity.
The following sections explain how the individual components can be described mathematically.

Bipolar / End Plate
The bipolar and end plates of the electrolyzer are treated as lumped capacitance masses C th with a uniform temperature T at every timestep.It is assumed that they are made of titanium.The temperature can be calculated through the following ordinary differential equation: The heat flows to the flow channels Q ̇BP,j are calculated in the respective cells.Q ̇amb,i describes the heat flows to the surrounding ambient.They can be defined as: For the heat transfer coefficients, values are taken from the work of Tjarks (2017).There, α bp = 2.5 W m -2 K -1 for the edges of the bipolar plates and α ep = 3.6 W m -2 K -1 for the end plates were calculated, which is in good accordance with the values determined for alkaline electrolyzers by Diéguez et al. (2008).

Anode / Cathode Flow Channel
Because of the production of hydrogen and oxygen in the MEA, gas-liquid flow occurs in the flow channels.The following section describes the mass balance, energy balance, and the calculation of heat transfer in them.
Pressure loss, on the other hand, is ignored and all fluids are assumed to be incompressible.Therefore, the momentum balance is not presented.
In the anode flow channel, incoming process water gets mixed with oxygen from the MEA.In addition, water flows into the MEA due to its electrochemical conversion and electro-osmosis.In the cathode flow channel, water and hydrogen enter from the MEA.Therefore, the mass balances of the anode and the cathode flow channel volumes can be described with equations 3 and 4.
In reality, hydrogen and oxygen cross the MEA as well due to diffusion and pressure differences.In this work, however, these mass flows are neglected in the calculation of flow conditions and heat transfer in the channels.They will be considered later when describing the gas flows into the flow channels, to be able to represent the effective hydrogen and oxygen production correctly.
The steady-state energy balance can be formulated via equation 5 with the sum of enthalpy flows into and out of the flow channel and the heat flows to the MEA and bipolar plate.A dynamic energy balance was not introduced since the storage capacity of the flow channel volumes is assumed to be negligible.Furthermore, it is assumed that both fluids in the respective flow channels have the same temperature when leaving the channel volume (T gas,out = T H2O,out ).
The heat flows Q ̇MEA and Q ̇BP can be described in analogy to equation 2. The fluid temperatures in the volumes T fluid are defined as the arithmetic average between the inlet and outlet of the flow channel volume (T fluid = 0.5T fluid,in + 0.5T fluid,out ).The overall flow channel volume temperature T fc again is defined as the arithmetic average between the gas and water temperature (T fc = 0.5T H2O + 0.5T gas ).
Because the flow channels are not in contact with the MEA and bipolar plate over the complete cell area A cell , the correction factor n A is introduced, which is assumed to be n A = 0.6 (Figure 4).Consequently, the contact area with the MEA becomes A MEA = A cell • n A .Since the flow channels are in contact with the bipolar plate on three sides, the contact area with the bipolar plate becomes A bp = 3A cell • n A if the flow channels are assumed to be quadratic.To reduce the model complexity, the heat transfer between the MEA and bipolar plate is neglected.The flow channels are assumed to have a thickness of w fc = 1 mm.phase were not present.Then, the heat transfer coefficient of the gas-liquid mixture α TP is calculated using the velocity ratio between the pure gas and pure liquid phase u r .There are three different formulations for α TP depending on the pure liquid's Reynolds nu er R LS : For 15 < Re LS < 175:   =   (1 +   ) 0.25 (6) For Re LS ≤ 15: = 0.75  (1 +   ) 0.25 (7) For Re LS > 175: Factor E is calculated using the Froude number of the pure liquid phase Fr LS (equation 9).For Fr LS > 10, it becomes E = 1.
The velocities in the flow channels are determined by their average single-phase volume flows and the total flow channel cross-sectional area.However, this velocity calculation is only valid if there is a parallel flow field design.Since the flow channels are assumed to be quadratic, the hydraulic diameter is d h = w fc .
The heat transfer coefficient for the liquid phase is calculated by the correlation of Sieder and Tate (1936) for the laminar and by the correlation of Dittus and Boelter (1985) for the turbulent regime, where L is the total channel length, λ LS is the thermal conductivity and Pr LS is the Prandtl number of the pure liquid (equations 10 and 11).Because the transition from laminar to turbulent occurs at significantly lower Reynolds numbers in gasliquid flow than in single-phase flow, the correlation for laminar flow applies only up to a Reynolds number of Re LS < 170.

Membrane Electrode Assembly
In the MEA, the conversion of water to hydrogen and oxygen takes place.The efficiency of this process depends on the cell voltage V cell , which can be described as the sum of the open-circuit voltage V ocv , the activation overvoltage V act , and the ohmic overvoltage V ohm (equation 12).The concentration overvoltage is neglected in this work because of its minimal effects at typical operating densities (Espinosa-López et al. 2018).
The open-circuit voltage describes the electromotive force that is required to start gas production (equation 13).It depends on the reversible cell voltage V rev , the partial pressures of hydrogen, oxygen, and water vapor in the MEA (p H2,MEA , p O2,MEA , and p H2O ), and its temperature T MEA .R = 8.3145 J mol -1 K -1 and F = 96 485 C mol -1 represent the universal gas constant and Faraday's constant.The reversible cell voltage can be described as a function of MEA temperature using the standard temperature T std = 298.15K and the reversible cell voltage at standard conditions V std = 1.23 V (equation 14).
To determine the partial pressures of hydrogen and oxygen in the MEA, their partial pressures in the flow channels p H2,cat and p O2,an must be quantified first.They can be calculated via Dalton's Law using the absolute pressures in the flow channels p an and p cat and the water vapor partial pressure (equations 15-17).The formulation for the water vapor partial pressure is taken from the work of Biaku et al. (2008)  When the electrolyzer is in operation, a pressure difference is established, so the partial pressures of hydrogen and oxygen in the MEA are higher than in the flow channels.The factors A cat and A an describe their linear dependency on the current density i (equations 18 and 19).In the work of Schalenbach et al. (2013), they are specified as A cat = 2.4 bar cm 2 A -1 and A an = 2.8 bar cm 2 A -1 .
The activation overvoltage describes the energy that is required to start the electrochemical reaction at the anode and cathode.According to Espinosa-López et al. ( 2018), it can be calculated only considering the activation overvoltage at the anode since it is significantly larger than at the cathode (equation 20).It depends on the charge transfer coefficient α an , and the exchange current density i 0,an .The latter is temperature-dependent and can be defined using an Arrhenius expression, with i 0,an,std being the exchange current density at standard conditions and )) The ohmic overvoltage describes the voltage loss due to the electrolyzer components' resistance to electric flow.Followin Oh 's law, the overvolta e is defined as the product of current density and the sum of electrical resistances.According to Olivier, Bourasseau and Bouamama (2017), it is valid to consider the e rane's electrical resistance R mem as the only resistance since it is the dominant factor (equation 22).The e rane's electrical resistance can be expressed in terms of the membrane thickness δ mem and its protonic conductivity σ mem .The membrane thickness is assumed to be δ mem = 183 µm, which is equivalent to the thickness of a Nafion TM 117 membrane (Chemours 2023).The protonic conductivity can be described with an Arrhenius expression as a function of membrane temperature, with σ mem,std being the protonic conductivity at standard conditions and E pro the activation energy required for the electron transport in the membrane (equation 23).
To define the E 's mass balance, the produced and permeated fluid flows have to be determined.The produced oxygen and hydrogen flow ṁ H2,prod and ṁ O2,prod and the consumed water flow ṁ H2O,cons can be calculated through the electrical current density and the respective molar masses M (equations 24-26).The Faraday efficiency is not introduced, since it depends mainly on the permeated mass flows, which are calculated separately.
According to Fick's law, the permeated oxygen and hydrogen flows can be described using the e rane's permeability to hydro en and oxy en ε H2 and ε O2 , its thickness, and the oxygen and hydrogen partial pressures in the MEA, assuming that the partial pressures of the permeated gases are small in comparison to the product gas partial pressures (equations 27 and 28).Schalenbach et al. (2013) determined ε H2 = 4.65 • 10 -11 mol cm -1 s -1 bar -1 and ε O2 = 2 • 10 -11 mol cm -1 s -1 bar -1 for the permeability of a Nafion TM 117 membrane at T MEA = 80 °C.In reality, these values are temperature-dependent, however, they are assumed to be constant in this work since electrolyzers are operated to a large extent at membrane temperatures close to 80 °C.
In addition, the water mass flow ṁ H2O,ed is transported from the anode to the cathode through the MEA due to electro-osmosis (equation 29).The factor n ed describes the percentage of proton transport through the membrane that involves water molecules.In the work of Santarelli, Torchio and Cochis (2006) It is assumed that the product gases exit the MEA completely saturated with water.According to Dalton's law, the water vapor mass flows ṁ vap,an and ṁ vap,cat can be calculated via the mass flows of the product gases and the pressure ratio between the water vapor partial pressure and the pressure in the anode and cathode, leading to the following equations for the water vapor mass flows: Consequently, the mass flows into and out of the MEA can be described by the following balance equations: To determine the E 's temperature and thus the product gas temperatures and the heat flows into the flow channels, the heat flow rate generated by the electrolysis reaction due to overvoltages Q ̇ely must be known.It depends on the current density and the difference between cell and thermoneutral voltage V tn = 1.48 V.
In addition, a latent heat flow Q ̇out,lat is removed from the MEA due to water saturation of the product gases.It can be determined by multiplying the vapor mass flows and the water enthalpy of vaporization ΔH vap = 40.65 kJ mol -1 .
The MEA itself is described as a lumped capacitance mass analogous to the bipolar / end plates.It is assumed that the porous transport layers represent its only relevant thermal mass.The material is assumed to be titanium.The porous transport layers are filled with a certain percentage of water.The porosity value of Φ = 0.37 is taken from the work of Grigoriev et al. (2009) and the total thickness is assumed to be w MEA = 1 mm.
The energy balance can be calculated from the sum of the incoming and outgoing enthalpy flows, the heat flow generated by the electrolysis reaction, the total heat flow into the flow channels Q ̇MEA,tot , and the latent heat flow removed by the vapor mass flows (equation 39).It is assumed that the temperature of the outgoing mass flows is equal to the E 's operating temperature (T fluid,out = T MEA ).

Model Validation
To validate the simulation model, experimental data from a 1 kW PEM electrolysis test stand at the City University of Applied Sciences Bremen was used.Although the measured data are not ideal for validating the present model due to a lack of large load steps, they can be used to demonstrate the general functionality of the model.The electrolyzer's technical specifications and simulation parameters required for the validation are listed in Table 1.Based on the stack dimensions it is assumed that the bipolar plates have a thickness of w bp = 3 mm and the end plates of w ep = 60 mm.The test series used for the validation primarily served to determine the currentvoltage characteristic of the electrolyzer.For this purpose, current densities from i = 0,1 A cm -2 to i = 2.5 A cm -2 were set and the cell voltages and water temperatures were measured.
Figure 5 shows the stack's current-voltage characteristic at T MEA ≈ 65 °C and the pressures listed in Table 1.To derive the missing parameters α an , i 0,an,std , E exc , E pro , and σ mem,std from this curve for the simulative mapping of the current-voltage relationship as described in section 2.4, the SciPy Python package was used.Its curve_fit function performs a non-linear least squares analysis to fit a set of m observations with a model that is non-linear in n unknown parameters.The fitted curve and the missing parameters are presented in Figure 5.They show a high agreement with the corresponding values from the literature review conducted by Espinosa-López et al. (2018).After finalizing the model parameterization with the calculated parameters, a simulation was performed using the input current profile with which the current-voltage characteristic was determined.The current profile and the resulting measured and simulated stack voltages are shown in Figure 6.It can be seen that in both simulation and measured data, the stack voltage drops and rises with decreasing and increasing current following the currentvoltage characteristic.The simulation result shows high accuracy with a mean absolute error of ΔV MAE = 0.088 V and a maximum deviation of ΔV max = 0.40 V. To verify the thermal modeling, the measured anode water outlet temperature was compared with the simulated one.
For the simulation, the measured anode inlet temperature was used.The result is shown in Figure 7.It can be seen that the temperature difference between the inlet and outlet in both simulation and measured data drops and rises as the current decreases and increases due to varying stack heat production.The simulated and measured temperature curves agree to a large extent, the mean absolute error is ΔT MAE = 0.239 K, and the maximum deviation ΔT max = 0.630 K.

Discussion
Using the obtained experimental data, the general model functionality was successfully demonstrated.However, minor differences between simulation results and experimental data were observed during the simulation of the anode water mass flow's outlet temperature.As the specifications of the employed temperature sensors were unknown, the deviations could be within their measurement inaccuracies.Furthermore, the water mass flow was not measured during the experiment but determined afterward by metering.Therefore, an incorrect mass flow rate may have been used for the simulation.This is also indicated by the fact that the stack's thermal energy balances calculated from the experimental data, once using mass flow and temperature difference and once using equation 38, show significant differences, which cannot be justified by additional heat losses or sources.
Having said this, it is important to note that the presented model is a work in progress.During the development and parameterization, a multitude of assumptions and simplifications were made, e.g. the cell element dimensions, the parallel flow field design, and the steady-state energy balance in the flow channels.Also, as mentioned before, it should be highlighted that the experimental data utilized for validation was not ideal for assessing the dynamic thermal behavior, as it lacked significant load variations.
Furthermore, the model incorporates calculations for substance transport through the membrane to realistically capture the quantities of produced oxygen and hydrogen.However, product mass flows, impurity gas fractions and cathode water mass flow were not measured during the experimental investigation.Therefore, it is necessary to conduct more comprehensive validation studies in the future, aiming to verify the accurate representation of all physical phenomena within the model.

Conclusion
In this paper, a quasi-2D model of a PEM electrolyzer was presented.The individual cells were discretized in flow direction of the flow channels, and the electrochemical and thermal behavior was described using analytical equations.Furthermore, the cells in the stack were thermally coupled to each other.
To validate the developed model, experimental data from a 1 kW PEM electrolyzer stack at the City University of Applied Sciences Bremen were utilized.Comparison of the measured data with simulation results demonstrated high accuracy in capturing the electrochemical behavior of the stack.Smaller deviations between thermal simulation results and measurement data can most likely be justified by imprecise data acquisition.
However, not all modeled physical phenomena could be validated using experimental data.In addition, a large number of assumptions were made regarding the design of the stack, which could not be substantiated.Further experimental investigation will be necessary in the future to comprehensively validate the model.

Figure 1 .
Figure 1.Basic structure of a PEM electrolysis cell.

Figure 2 .
Figure 2. Structure of a single cell.

Figure 3 .
Figure 3. Structure of the stack.

Figure 6 .
Figure 6.Input current and measured vs. simulated stack voltage.

Figure 7 .
Figure 7. Anode inlet temperature and measured vs. simulated anode outlet temperature.
E exc the activation energy required for the electron transport in the anode electrode (equation 21).
it was given as n ed = 0.27.