Experimental validation of a dynamic modelling of a Reversible Solid Oxide Cells (rSOCs)

. The integration of Hydrogen technologies in different end-uses such as transport, electric microgrids, residential and industrial applications, will increase exponentially soon. Hydrogen as energy carrier allows more favourable energy conversion than other conventional systems and is crucial in worldwide decarbonize end uses. The production of green hydrogen, using RES, is a key area for the evolution of this technology. In this context, SWITCH is a Horizon 2020 European Project that aims to design, build and test an in-situ fully integrated and continuous multisource hydrogen production system, based on solid oxide cell technology. Reversible Solid Oxide Cell (rSOCs) technologies allow to convert renewable energy as hydrogen in the power-to-gas application (P2G) and in reversible mode is able to produce electricity from hydrogen stored, power-to-power application (P2P). rSOCs are really interesting to stabilize the random nature of RES because a combined electrolysis and fuel cell system should be able to switch between the two modes as quickly as possible in order to optimize the integration and the use of RES. However, rSOCs need a complex BoP from the thermal point of view, able to guarantee high efficiency even at partial load mode as well as easy start-up and shutdown procedures. In this work, a Stack Box Module dynamic model was developed in Modelica environment as a dynamic tool for the definition and optimization of BoP requirements. Stack model was validated in SOFC (Solid Oxide Fuel Cell) and SOE (Solid Oxide Electrolyser). The results of the simulation provide verification of the technical/thermodynamic behaviour and flexibility of a stack box of 70 cells. Dynamic modelling allows to evaluate the effect of the reagent inlet temperatures on the operation and hydrogen production/consumption in terms of yield as well as the transients between the different operative modes. Model has been validated by experimental measurements performed in the laboratory. In particular, the kinetics of the reactions governing steam methane reforming (SMR) was considered from data found in the literature, while the ASR (Area Specific Resistance) value was calibrated according to experimental data. The results of the dynamic model show as model can be a useful design and optimization tool for the SOCs technology.


Introduction
According to the European Hydrogen Strategy, mass production of electrolysers is expected to be deployed on the market with capacities of 6 GW by 2024 and 40 GW by 2030 [1].Reversible Solid oxide cells (rSOCs) occurs in electrochemical devices that operate at high temperature (800-1000 °C).They are particularly attractive as electrochemical devices since they can operate both as fuel cells and as electrolysers depending on the application and needs.The behaviour of both Solid oxide fuel cells (SOFCs) and Solid oxide electrolyser (SOE) has been studied and reported in many articles [2][3][4] [5] and the reversibility has been investigated since the 1980s when both steam [6] [7].The reversible working mode is gaining interest due to its high flexibility of operation than which supports the integration of fluctuating energy produced from renewable energy sources (RES).In fact, rSOCs can work as electrolysers when there is excess production from RES.The hydrogen produced during periods of RES abundance can be stored in dedicated systems until RES can meet the energy demand, at which point the stored hydrogen is used for energy production [8].However, a limited amount of studies has been carried out on the combined design of a reversible fuel cell characterized by SOFC and SOE behaviour like the one that the SWITCH project aims to implement.[9][10] [11].In this study, values of area specific resistance at temperature T0 (T0 =1073K) ASR0 and activation energy (Ea) valid for the calculation of area specific resistance (ASR) for the SOE and SOFC cases respectively were extrapolated through laboratory tests and used to validate the model and the simulation results of rSOFC behaviour.

Methodology
The dynamic modelling of the Reversible Solid Oxide Cells (rSOCs) was developed by using Modelica language [12] with the Dymola software [13].The one-dimensional models and sub-models that characterizes the system is discretized over space with a number of n intervals.The choice of this discretization number is set at the top level of the model's system and then it results equal for each subsystem and sub model.In addition, a kinetic model for the reactions describing the steam methane reforming was considered within the anode channel of the stack and within an external reformer while a reversible equivalent electrical model for the water electrolysis and fuel cell operation was considered within the elementary cell model.The general layout used in the model is presented in Figure 1.The general layout consists in different blocks as the stack, the external reformer, the mixer, the hydraulic losses blocks, the temperature sensors, the sinks and the fluid sources.In the hydraulic losses blocks the pressure difference between the inlet and the outlet port is caused by friction losses.In the friction model the mass flow rate (mflow) is in first approximation function of nominal operating point specified by mass flow rate (mflow0), density(ρ0) and pressure drop (dp0 ) : The diagram and the breakdown of the stack is shown in Figure 2. The Stack model (1) consists of other sub-models such as one relating to thermal losses (4) with the environment and one relating to the sub-stack model( 2).The sub-stack is itself composed of a model relating to the anode and cathode channels (3) and one relating to the elementary cell model ( 4), as is possible to see in Figure 2. In the model that describes the stack it is possible to link several sub-stacks in parallel.The mass ports (orange circles) used for the flow of fuel and air to the anode and cathode, respectively, are connected to particular manifolds that consider the number of sub-stacks chosen.
According to SOE/SOFC modes, gases can exit or enter the sub-stack and then are collected in the exhaust manifolds, at the exit of the stack.The heat exchange between manifolds and the sub-stack is accounted in the manifold model.Moreover, a medium (NASA reformate long) from the Modelon fuel cell library [14][15] was used for the fuel that feeds the anode in both SOE and SOFC phases, while a medium with a fixed composition was used for the air.In the SOFC phase a mixer is used to mix steam and methane in order to have a Steam Carbon (S/C) ratio of about 2.3, while in the SOE phase water is used with a mass fraction of about 12% to avoid undesirable phenomena at the cathode that can cause irreversible damage [16].The compositions and mass fractions of inlet fuel and air for the SOE and SOFC case are expressed in the following Table 1 and Table 2.A cell area Acell of 80 cm 2 I provided by the LSM during SOFC operation ISOFC and that absorbed by the SOE mode ISOE is an input data of the model, taken by experimental data.The current sign indicates the outgoing (positive) or incoming (negative) direction as per the standard.The number of H2 moles, ݊ ு మ , is calculated considering the Faraday's law as follows: Where ‫ܨ‬ ௨ is the fuel utilization of 0.75, F is the Faraday constant ‫ܨ(‬ = 96485.3‫,)݈݉/ܥ‬ ncell is the number of cells (70), and nsub is the number of sub-stacks.The input for anode channel mass flowrates of the fuel are defined as follows: Where ‫ܯܯ‬ is the molar mass of the different elements.

Cell model
The single elementary cell model is the lowest modelling level developed for this study.The single cell can be traced back to an electrical circuit with voltage and resistance generator.The voltage delivered from the cell is closely related to the electrical resistance of the cell, the current density, and the molar flows coming from the upper level of the sub-stack.In the elementary cell model, the open circuit potential at various temperatures can be evaluated with Nerst's equation, considering the evolution of the oxidation-reduction reaction: With the Nerst's equation: where F is the Faraday constant ‫ܨ(‬ = 96485.3‫,)݈݉/ܥ‬ R is the universal constant of gas ( ܴ = 8.314 ‫ܬ‬ /K mol) [17]while ‫ܩ߂‬ is the Gibbs free energy for the complete reaction (1) calculated by using the difference between the enthalpy formation ( ‫ܪ߂‬ ) and the product of cell temperature and formation entropy (ܶ ߂ܵ ) of both anode and cathode semi-reactions [18].In addition, for the calculation of the enthalpy and entropy of formation, polynomials which are functions of the cell temperature were used [19].The partial pressures of both reactants and products are evaluated for each segment of the onedimensional discretization with a penalty to avoid negative pressure values : In this way, by using Nernst equation, the cell equilibrium potential is calculated as a function of temperature and partial pressures.In real conditions, irreversible voltage losses occur when an electrical load is connected to the cell and a current flow through the cell.The overall voltage losses can be divided in three main categories: ohmic losses related to the resistivity of the solid oxide cell materials, electrode activation losses, and concentration losses.In this model, only an ohmic loss term (ߟ ) is considered and represents all the voltage losses occurring in the SOE/SOFC model.The voltage at the cell external connectors is calculated by Kirchhoff voltage law (KVL) on the circuit [20], shown in Figure 3, as follows: The second term is positive for SOE and negative for SOFC since the resistances must be overcome to allow the reactions to occur.Each cell is discretized by dividing the spatial domain in N equivalent elements cells (i.e.݊ = 70).Therefore, the single-cell resistance is given by the sum of the resistances over the N cell elements.As a consequence, the computation of the overall stack resistance (ܴ ௧ ) is given by : Where ASR is the area-specific resistance of each cell and Acell is the equivalent active area of the stack.An empirical law [4] is used to describes the ASR as a function of the cell temperature (Tcell): where ‫ܴܵܣ‬ is the area specific resistance at temperature T0 (T0 =1073 K) [9].The data extrapolated from the laboratory experiments are expressed in Table 3 Once the electrical model is defined, the current density (jionic) is defined as follows: Where Icell is the current applied to the whole cell.Depending on the current direction, the ASR0 and Ea values associated with the electrolyser or fuel cell mode are considered within the elementary cell model.Therefore, the equivalent electrical circuit with the voltage generator represented by the potential of Nerst and an electrical resistance that depends on the current (pin_n.i)and the system temperature through the ASR, is represented in the block diagram in Figure 3.The current generator (source.V) is connected to the electric connectors (pin_p, pin_n), passing through the internal resistance block (internalRes).Tcell is defined as the temperature of the cell that is directly connected to the thermal port wall (wall.T = Tcell).A thermal port allows temperature and heat flow continuity.Moreover, two mass ports that allow only the stoichiometric mass flow to diffuse from the cell to the channels and vice-versa according to hydrogen consumption or production.The temperature of the cell is calculated by the following expression for each element of discretization : where Mcp is the thermal capacity (400 J/K), Qan is the heat flux between the cell and the anode, Qcat is the heat flux between the cell and the cathode, Qcell is the heat produced in the i-th cell, and Qwall is the heat exchanged between cell and external environment.The thermal ports referred to heat transfer between the cell and the anode or cathode (wall_an and wall_cath, respectively) are characterized by temperatures (Twall_an and Twall_cat).The heat fluxes are defined: where Kc = 250 W/(m 2 K) is defined as the heat transfer coefficient between fluid and substrate and Acell = 80e-4 m 2 .
ௗ ௗ௧ represents the speed with which reagent species are consumed and produced, i.e. the kinetic speed of reaction (mol/s).The global reaction introduced before can be split in the following two semi-reactions : The mass flow rate of H2O consumed can be evaluated according to Faraday's law, which allows at the same time to evaluate the O2 and H2 produced: The sign is positive (+) for substances entering into a block and negative (-) for those exiting the block, i.e. positive for H2O and negative for H2 and O2.Therefore, it is possible to determine the flow through the mass ports and the bulk enthalpy of the mixture ‫ܪ(‬ ௧ ̇ and ‫ܪ‬ ௗ ̇), since specific enthalpy of each component is known.
Finally, Qcell is defined through the following energy balance: Where ܲ = ܸ ‫ܫ‬ is the cell power while Vcell is calculated as potential difference between the two electric connectors (pin_p, pin_n).

Reaction channel model
In the SOFC mode, the reactions describing the steam methane reforming will take place in the anode channel and in the external reformer.Subsequently, the hydrogen produced by these reactions will be converted into electricity produced from the fuel cell by using the equivalent electrical approach described in the elementary cell model.The stoichiometric coefficients of the following reactions are defined at first in the reaction channel: The reaction rates are expressed as R1, R2 and R3, respectively and defined as follows [4]: Where ܵ ௧ is a coefficient that represents the effect of a possible catalyzer, and in this model, it has been fixed with a different constant value in external reformer and in internal anode channel of the stack (4e-5 and 1e-5 respectively).Ω is a parameter expressed by the following equation: Moreover, kj is the Arrhenius reaction constant and is defined as follow: The values of kj are calculated by using the parameters listed in Table 4 [9] ‫ܭ‬ ூ , ‫ܭ‬ ூூ and ‫ܭ‬ ூூூ are the equilibrium constant related to the corresponding reactions and are expressed as a function of the cell temperature.On the other hand, the Van't Hoff species adsorption constants ( ‫ܭ‬ ) introduced for the definition of Ω are calculated with a similar Arrhenius expression by using the parameters listed in Table 5 [11].The reaction molar rate (rZ) are expressed as follows: 0.07, v ଶ = 0.06 and v ଶ = 0.7.These reaction molar rates are considered in the consumption or production rates of each substance by the following equation: Where deplZ is a protection factor used against depleted species.It is calculated by using the smoothing splice function present in the math library of Modelon library while ‫ݔܼ‬ ௪ is the molar fraction and is defined as: In which MMX is the molecular weight of the medium while ݉ܺ ௪ is the mass flow rate that is defined in the reaction channel as the product of mass flow by the mass fraction.݉ܺ ௪ = ݉ ௪ ܺ ( 22)

Results and Discussion
In this section, simulation results are presented in comparison with experimental results.The results referring to the behaviour of the rSOCs refer to SOFC behaviour with a current ramp delivered from 1 to 9 A and SOE behaviour with a current ramp from 0.5 to 3.5 A. In the first case, an increase and decrease of the supplied current within the declared limits was considered, with a trapezoidal trend, while in the second case only a linear increase of the absorbed current was considered.Simulation times in accordance with those of the experiments have a duration of 3450 and 22900 s for SOFC and SOE, respectively.The voltage trend as a function of time is presented in Figure 4 for the SOFC case (top) and for the SOE case (bottom).It is possible to notice that the trend of the simulated voltage (Voltage Sim) in both cases is very consistent to that of the experimental voltage (Voltage Exp).In the SOFC case there is an increase and subsequent decrease in voltage of approx.8.5 V, while in the SOE case there is an increase of approx.7 V.The Figure 5 shows the relative error (gap) between the simulated and experimental voltage in the SOFC case (top) and in the SOE case (bottom).It can be observed that in both cases, the relative error is less than 3%.In the SOFC case when the current reaches plateaux and is stationary with a value of 9 A, the relative error increases from a value close to zero to a maximum value slightly higher than 2.5%.In the SOE case instead, the relative error varies from a minimum value of about 0.25% to a maximum value of about 2.25%.In Figure 6, the evolution of voltage versus current (V-I trend) in the SOFC case (top) and in the SOE case (bottom) is presented.The trend of the polarization curve allows appreciating the nature of the ohmic losses.These in fact present a typically linear trend both in the experimental case (linear lines) and in the simulated case (dotted lines).Particularly in the SOFC case, it is possible to appreciate a difference in the V-I trend in the ramp-up and ramp down phases.For the first phase, the lines have a blue and red (dotted) colour, while for the second phase they have light blue and orange (dotted) color.It can be noticed that the curves have a parallel trend in the current range from 3 to about 7 A, while outside this range the curves have a slight deviation from the experimental data.For the SOE case, on the other hand, the curves show a linear trend for the complete current range.Finally, with this model it was possible to demonstrate that the system can work reversibly as a fuel cell or electrolyser.The calculation and the relative error estimation allow the model to be considered robust as the simulation results converge with the experimental ones with a maximum margin of error about 0.028.

Conclusions
In this paper a dynamic model of a reversible Solid Oxide Cells (rSOCs) has been developed in Dymola [13].The dynamic model results were validated through experimental data with a relative error of less than 3 % The model constructed in this way is enough precise to simulate rSOCs for different applications as in power to power system, or to study the integration of rSOCs systems in electrical grid to manage in smart way the use of renewable energy.The reversible behaviour between electrolyzer (SOE) and fuel cell (SOFC) has been implemented in the elementary cell sub-model using experimental data.This would allow to simulate the storage of a surplus of energy produced by renewable sources [21] in the form of hydrogen and vice versa the conversion of hydrogen into electrical energy The use of SOCs systems at high operating temperatures shows a low impact of activation overvoltages compared to low temperature systems [9] and these systems show negligible mass transport effects at high currents and low fuel utilization (steam or natural gas).Moreover, the experimental data used for the validation of the model did not cover this type of case and for these reasons, the use of a simple ohmic model does not introduce significant errors.However, since the ASR used is determined by experimental data and is a function of temperature, inaccuracies in the model may arise if current and fuel utilization values are used outside the range of the experimental values.To achieve greater flexibility in the behaviour of the modelled current-voltage curve, it would be appropriate to implement more complex functions to describe the overvoltage losses due to activation and mass transfer [22][2].However, in reality, the operating condition of an LSM is usually handled within a certain current range for which the ASR function used can be adjusted although there are obviously deviations that are considered negligible.

E3SFigure 1 :
Figure 1: Layout of the general model

Figure 2 :
Figure 2: Diagram and breakdown of the model that describes the stack

Figure 3 :
Figure 3: Block diagram of the i-th elementary cell element

Figure 4 :
Figure 4: Voltage trend as a function of time for SOFC (top) and SOE (bottom) respectively

Figure 5 :
Figure 5: Relative error trend as function of time for SOFC (top) and SOE (bottom) respectevely

Table 1 :
Mass Fraction of the fuel in SOE and SOFC mode

Table 2 :
Mass Fraction of the air 3. Concentrated parameters model, physical phenomena occur only in the respective component ; 4. Diffusion phenomena are not considered.

Table 3 :
Experimental value of ASR0 and Ea for SOFC and SOE case