Dynamic modelling of a direct internal reforming solid oxide fuel cell stack based on single cell experiments

with direct predicted by the two the need to identify rate limiting steps for the reforming and hydrogen oxidation reactions on anodes of functional SOFC assemblies. The modelling approach can be used to study off-design conditions, transient operation and system integration, as well as to develop adequate energy management and control strategies.


Introduction
Global agreements to eliminate the net emissions of greenhouse gases and hazardous air pollutants call for the development and implementation of more efficient and cleaner energy conversion technologies [1][2][3]. Solid oxide fuel cell (SOFC) systems enable fuel to electricity conversion with high efficiencies and virtually no hazardous air pollutant emissions, and even higher efficiencies are achieved when they are combined with thermal cycles [4,5]. Moreover, SOFCs can be operated with non-hydrogen fuels with minimal pre-processing requirements and have higher resilience to fuel impurities due to their relatively high operating temperature compared to other fuel cell types [6,7].
Another advantages of SOFCs over other fuel cell technologies is their capability of direct internal reforming (DIR) of energy dense hydrocarbons to a hydrogen rich gas mixture on the fuel electrode [8]. An example of this is direct internal methane steam reforming (MSR) in SOFC systems fuelled with natural gas using the heat and steam produced by the electrochemical hydrogen oxidation reaction. Moreover, DIR reduces the need for external reforming as well as cooling with cathode air [9]. However, DIR imposes the risk of solid carbon formation on the anode, and can be expected to increase temperature gradients within the stack [10]. The resulting thermal stresses may inflict damage to the brittle ceramic cells [11].
The thermal stresses in SOFCs are primarily affected by the stack design, in-and outlet temperatures of fuel and air, fuel composition and electric load [12,13]. The spatial distributions of the electrochemical and reforming reactions play an important role, as these determine the local heat consumption and production and, as such, the temperature profile within the stack [11,14]. Therefore, detailed information on the kinetics of those reactions in the SOFC anode is required to accurately model their spatial distributions, which is useful to develop appropriate control strategies [15]. Kinetic models for the MSR reaction on SOFC anodes vary from simple first order dependencies on the methane partial pressure to multi-step reaction mechanisms including up to 42 individual reactions on the catalyst surface [16,17]. Multi-step mechanisms capture the elementary reactions on the catalyst surface and can thus be applied to any anode if its morphology and all temperature dependent rate constants are known [18]. However, this requires a high level of modelling detail and extensive experimental effort, while first order models can be used in basically any model with determination of only one temperature dependent rate constant. However, multi-step methods are expected to give more accurate results.
The kinetics of both the electrochemical and reforming reaction depend on the local temperature, partial pressures of reactants and products and the catalytic properties of the anode. However, as both reactions in turn affect the spatial distributions of temperature and partial pressures, they are strongly coupled. Although the individual reactions have been studied extensively on substrate materials, deviating operating conditions and cell materials make it difficult to directly implement these kinetic mechanisms in SOFC stack models [19][20][21].
A number of studies have numerically modelled DIR SOFCs, usually to evaluate their electrochemical performance or develop appropriate control strategies. The MSR reaction is most commonly assumed to be in chemical equilibrium and thus considered to proceed infinitely fast. Alternatively, authors have implemented kinetics which assume a first order dependency on the methane partial pressure only [22][23][24]. The multi-step reaction mechanism derived by Hecht et al. [17] for nickelyttrium stabilised zirconia (Ni-YSZ) cermet anodes has been implemented in more comprehensive computational fluid dynamics (CFD) models [25,9].
In control-oriented dynamic SOFC models the level of detail is typically limited to reduce the computational demand. Multi-step reaction mechanisms are, therefore, less suitable for such models, while simple first order kinetics may yield inaccurate results [18]. Alternatively, Langmuir-Hinshelwood, Hougen-Watson (HW) or Eley-Rideal kinetics can be used if the reaction can be described by a single rate determining step on the catalyst surface. Power law (PL) kinetics may allow the inclusion of the effects of other reactants and products on the reaction rate without knowledge on the reaction mechanism [26][27][28]. However, it is unknown if any of those kinetic mechanisms can be used to accurately predict the spatial distribution of the reforming reaction within the stack.
Acknowledging the need for reforming kinetics derived at conditions relevant for stack operation, MSR experiments were carried out on a functional single cells with a nickel-gadolinium doped cerium oxide (Ni-GDC) cermet anode in a previous study [21]. The data was then used in a follow-up study to derive PL and HW kinetics [29]. In addition, a CFD model of the test setup was developed and both kinetic models where implemented. Although both models accurately predicted the overall methane conversion, different spatial distributions of the MSR rates were obtained. However, whether this is also the case for stack operating, where the conditions differ substantially from the controlled environment in single cell experiments, is not clear.
Since the DIR experiments were carried out on the same commercially available single cells applied in stacks, direct transfer of the kinetics to stack models is in principle possible. Two 1D dynamic SOFC models are developed in this study, one for the single cell test station and the other representing a stack based on the same type of cells. The structure and equations are the same for the two models, but the boundary conditions are adjusted to account for the physical differences between single cell setup and stack operation.
The stack model developed in this study not only includes the active area of the stack, but accounts for heat transfer in the inactive in-and outflow sections and heat losses to the surrounding as well. Therefore, a commercially available integrated stack module (ISM) is modelled, which contains two 30 cell stacks in series with internal fuel manifolding, placed in a thermally insulated box and external air manifolding [30]. The ISM is equipped with connections for the fuel and air in-and outlets, electric power cables and temperature monitoring [31]. However, the ISM does not include balance of plant components such as a (pre-) reformer, afterburner and heat exchanger, as is the case in other ISMs [32,33].
The ISM model is validated with operating parameter and power curves published by the manufacturer for hydrogen-nitrogen, catalytic partial oxidation (CPOX) and steam reformate (SR) fuels. The single cell model is validated using experimental data and the CFD model developed in previous work [21,29]. Finally, the two DIR models obtained in the single cell experiments are implemented in the stack model to compare the predicted spatial distributions of species concentrations,  temperatures and reaction rates. A schematic overview of this approach is shown in Fig. 1.

Model description
This section describes the 1D modelling framework developed to simulate single cells, stacks and integrated stack modules based on a single set of equations. Separate control volumes are defined for air, fuel, interconnect and the positive electrode-electrolyte-negative electrode (PEN) assembly, which are then discretised in the flow direction. The model formulation is dynamic to enable application in future transient simulation and development of control strategies.
It is commonly assumed that the solid temperatures dominate the transient behaviour of SOFCs and should thus be solved dynamically, while the remaining equations are can be formulated quasi-static [23,24]. However, the accuracy of the time dependent solid temperature depends on the time interval at which the quasi-static properties are updated in that case. Moreover, quasi-static mass balances require iterative solving, while stiff solvers can deal with the different time scales encountered, yielding acceptable runtimes. Therefore, all mass and energy balances are implemented dynamically.
Hydrogen is electrochemically oxidised at the anode of the SOFC by oxygen ions diffusing from the cathode through the dense oxide electrolyte. Two electrons are released by this reaction (1) subsequently traveling to the cathode through the external circuit, thus delivering a useful current. The electrons are consumed in the reduction of oxygen in the cathode: The overall reaction, thus results in the electrochemical oxidation of hydrogen in the anode compartment by oxygen from the cathode side, delivering a net electronical current. Light hydrocarbons can be converted to a hydrogen rich mixture via the steam reforming reaction. For methane this gives: Although carbon monoxide can be directly electrochemically oxidised on SOFC anodes, the hydrogen oxidation reaction is generally reported to dominate the electrochemical oxidation practice, while carbon monoxide reacts via the water gas shift reaction [34]:

Mass and energy balances
Dynamic mass and energy balances are implemented in the model. The gas channels of the SOFC are modelled as a series of continuously stirred-tank reactor control volumes, in analogy to Hosseini et al. [35]. The local time derivative of the molar concentration a species follows from a molar balance, divided by molar capacity of the control volume, determined by its size and ideal gas law: The outlet flow of species i is calculated from its concentration multiplied by the total molar outflow, which follows from the total molar inflow and the sum of moles produced and consumed by reactions. It is assumed that the changes in the total molar flow settle infinitely fast and can thus be assumed quasi-static.
It is assumed that all chemical reactions take place on the solid anode and cathode catalyst interface of the SOFC, and the heat from chemical reactions is, therefore, assigned to the PEN control volume and not to the gases. The dynamic energy balance of a gaseous control volume thus only depends on the change in enthalpy of the inlet gas flow due to the temperature difference with the previous control volume and the convective heat transfer to or from the solid parts, divided by the heat capacity of the control volume while the enthalpy change from inlet to outlet due to chemical reactions is thus accounted for in the energy balance of the PEN. The thermodynamic properties of the gases, such as heat capacities, enthalpies and entropies are calculated using the Shomate equation with coefficients provided by the National Institute of Standards and Technology Chemistry WebBook [36,37]. The dynamic energy balance of the solid control volumes has the following form: Here the first two terms on the right hand side represent the convective heat transfer to or from the gases and heat conduction, and apply to both the interconnect and PEN, including its inactive area. The sum of the heat of reactions and electric power drawn from a control volume only apply to the active parts of the PEN. The heat loss to the surroundings Q loss is applied to the boundaries and zero everywhere else, as is discussed in Section 2.6.

Chemical reactions
The chemical reactions considered to take place on the anode of the SOFC are the WGS and MSR reaction. Since the WGS is assumed to proceed infinitely fast, its reaction quotient Q WGS is assumed to be equal to the chemical equilibrium constant K WGS along the active area of the SOFC which is achieved by selecting an arbitrary high value for the frequency factor of the WGS reaction. Assuming ideal gas behaviour for the reactants and products the fugacity constant 1 and the activities are calculated from: where p 0 is the standard pressure. A similar approximation may be used for the MSR reaction if the kinetics are not known, yielding: The equilibrium constants in Eqs. (9) and (11) are obtained from the Gibbs free energy change of the reaction at standard state: However, a kinetic model is required if the reaction does not proceed infinitely fast, which is commonly reported to be the case for methane. Therefore, two kinetic MSR models were derived in a previous study using data obtained from experiments on an electrolyte supported cell with a Ni-GDC anode [21,29]. One is a global kinetic model of the PL type and depends on reaction orders and for the reactants methane and steam respectively and the activation energy E a : The second kinetic model is of the HW type, assuming that the dehydrogenation of the CHO radical on the catalyst surface is rate determining. This mechanism was proposed by Xu et al. [38] for steam reforming reactors with nickel catalysts and yields the following expression for the MSR rate: where the denominator is a Langmuir adsorption isotherm accounting for the surface coverage by adsorbed oxygen from a steady state of steam adsorption and hydrogen desorption, leading to: In this isotherm, A 0 is pre-exponential factor of the temperature dependent adsorption equilibrium constant with associated enthalpy H O . Four parameters should be obtained from experimental data in both kinetic models. While these are fully independent in the PL expression, two temperature dependent constants are required for the HW kinetics as the dependence on the reactant and product partial pressures depends on the rate limiting reaction step.

Electrochemical reactions
The spatial distribution of the electrochemical reaction rate, more commonly referred to as the current density, is calculated using the equipotential assumption, i.e. the voltage is uniform on the cell plane. The bisection algorithm is used to determine the resulting cell voltage for the total current drawn. This means that the current density distribution is calculated for which the sum of the overpotentials equals the difference between cell voltage and the electrochemical equilibrium potential, better known as the Nernst voltage: Assuming that the kinetics of hydrogen oxidation dominate the electrochemical reaction, the Nernst voltage is given by: The ohmic resistance of the PEN structure depends on the individual resistances of the anode, cathode and electrolyte, and is proportional to their thickness and disproportional to their electronic or ionic conductivity. In addition, a contact resistance R contact is included to account for non-ideal electrical contacts in the single cell test station and stack assembly: It is assumed that the electrical conductivities of the anode and cathode can be estimated with constant values for the operating condition of the stack, since their temperature dependence is limited and their conductivity is high for SOFC operating temperatures. The ionic conductivity of the electrolyte, on the other hand, has a strong temperature dependency and is thus calculated for the local PEN temperature. A reasonable value for the contact resistance of a well-designed stack is obtained from Liu et al. [39], while the contact resistance for the single cell test setup was estimated based on an experimental IV-curve.
Concentration losses typically start to dominate the electrochemical losses only for higher current densities. Even if the overall current density is limited, local current densities can be high, especially toward the outlet of a stack where the temperatures are usually highest. The departure of the cell voltage from Nernst as a result of the concentration gradient in the electrode can be shown to follow from [22].  (19) in which p i tpb , is the partial pressure of species i in the triple phase boundary, calculated from the bulk species partial pressure, corrected using the local current density, electrode thickness and an effective diffusion coefficient D eff : The activation losses in the electrodes are calculated using the Butler-Volmer equation. Assuming symmetry between the anodic and cathodic reaction, the overpotentials follow from [40]: = Although mathematically convenient, it should be noted that the anode exchange current density j 0 in the Butler-Volmer equation depends on the concentrations of the reactants and products of the electrochemical reaction. This is partly due to the equilibrium potential effect discussed by Bessler et al. [41], since the forward and backward reactions on the electrodes are affected by the electrical potential difference, which depends on the reactant and product concentrations. The global influence of the reactants and products concentrations on j 0 can be shown to depend on these concentrations themselves and the rate limiting charge transfer mechanism [41].
It should be noted that Eq. (23) is only valid if a single charge transfer reaction is rate limiting and the reactant and product concentrations are constant, which is generally not the case for cermet electrodes. In order to deal with the limited understanding of the fundamental electrochemical oxidation kinetics, it has become customary to introduce global reaction orders for the reactant and product activities in the exchange current densities [42]: The values of the reaction orders , and depend on the rate limiting kinetics of the electrochemical reactions, temperature, reactant and product concentrations and absolute electric potential difference. A value of 1/4 is often used for , but for and values from 0 to 1 and −1/2 to 1 respectively can be found in literature [41,43]. Since the partial pressures of hydrogen and steam vary substantially along the anode, especially under internal reforming conditions, values for and are fitted using power curves provided by the manufacturer.

Effective diffusion coefficients
Effective diffusion coefficients are required to calculate the mass transfer losses in the anode and cathode. Diffusion in the porous SOFC electrodes has been extensively studied and various models have been presented, usually accounting for multicomponent diffusion and wallsurface interaction in the channels of the electrode material. In this model, effective diffusion coefficients for hydrogen and steam in a multicomponent gas mixture and oxygen in air are calculated within the porous electrodes, since they are only used to calculate concentration overpotentials for limiting current densities.
The approach in this study is similar to the one proposed by Chan et al. [44]. Effective diffusion coefficients for hydrogen and steam in a mixture are calculated from: Since the pore diameters in the porous electrodes of SOFCs are typically comparable to the mean free path of the gas molecules, Kundsen diffusion should be accounted for. The Bosanquet formula is used to calculate an effective diffusion coefficient for species i, which is then corrected for the tortuous path of the molecule and the porosity of the electrode: In analogy with Yakabet et al. [45] the Knudsen diffusion coefficient for species i is calculated from kinetic theory where r is the mean pore radius of the electrode and M i the molecular mass of species i. The molecular diffusion coefficient of species i in a multicomponent mixture is calculated from using Blanc's law, in which the binary diffusion coefficients D ij are calculated using Fuller's method [46,47]:

Heat transfer coefficients
Heat transfer in SOFCs proceeds via conduction in the solid parts, convection between the gases and the solids and radiation. It has been shown that radiative heat transfer can be ignored in most cases for planar stack designs, since the temperature gradients perpendicular to the cells are usually small [48]. The convective heat transfer coefficient between the solid parts and the gases in the SOFC can be calculated from the Nusselt number which is assumed to be independent of the Reynolds number due to the laminar flow conditions, and has a constant value of 3.09 [22]. The hydraulic diameter of the square gas channels follows from their width w ch and height ch where the width of the channels is calculated from the width of the cell, number of channels and interconnect thickness: The heat conductivity of the gas mixture is calculated using the Wassiljewa equation [49] = y y , where i is the thermal conductivity of the individual gas species, calculated from an estimation method recommended by Todd et al. [47] = C T 0.01 1000 , i n n n (35) and ij a function depending on the thermal conductivity and molecular mass of the species involved, proposed by Mason et al. [50]:

Model parameters and boundary conditions
In principle, all equations discussed so far apply to both the model of the single cell test setup and the stack. The difference between the two models is in the geometrical parameters and boundary conditions assumed. The cell parameters are based on ESC2 cells obtained from Kerafol/H.C. Starck, for which MSR kinetics were derived in previous work [21,29]. These cells have a Ni-GDC anode and an 8YSZ/LSM-LSM double layer cathode supported on a dense 3YSZ electrolyte [51]. The stack model is based on Mk200 stacks and the ISM V3.3 produced by Sunfire/Staxera, which relies on the same cells [52]. The parameters assumed in the model are summarised in Table 1. Fig. 2 illustrates the differences between the ISM and single cell configuration. The Sunfire/Staxera ISM V3.3 contains a single tower of two Mk200 stacks consisting of 30 cells each with internal fuel manifolding. The stack is placed in a thermally insulated box which provides external air manifolding and has connections for the fuel and air in-and outlets, electric power cables and instrumentation. The ISM offers a flexible solution for system integrators, since other balance of plant components, such as heat exchangers, blowers or (pre-) reformers, are not included. An inactive in-and outflow manifold is present in the stack where only heat transfer is expected to occur, since the cells are contained in metal cassettes in the stack. The temperature of the ISM is sustained by the heat produced by the electrochemical reaction. Overheating of the stack is prevented by control of the cathode air flow, although heat will be lost through the insulation as well. The active area of the stack is discretised into 50 control volumes, and the inactive in-and outlet areas are subsequently divided into 21 control volumes each.
In the single cell test setup 10 × 10 cm cells are placed in a ceramic holder with fuel and air manifolding. Compression seals are used to seal the anode and cathode compartments, and nickel meshes are used as current collector. The ceramic holder is placed in an isolated furnace equipped with electric heaters and temperature control. This is used to maintain a constant cell temperature in the ceramic block during the MSR experiments. The active area of the single cell had to be discretised into 250 control volumes to accurately capture the sharp temperature gradients encountered at the inlet, a result of the prescribed wall temperature boundary condition.
Spatial variations perpendicular to the flow direction are ignored, since the model is formulated in 1D. Therefore, only boundary conditions on the inlets, outlets and perpendicular to the cell assembly need to be defined. A periodic boundary condition is assumed on the interconnect in the stack model, thus assuming an infinitely repeated stack assembly. Furthermore, it is assumed that the temperature gradient in the boundaries is negligible due to the insulation material applied: The heat loss to the environment is calculated from Newton's cooling law, and thus proportional to the temperature difference with the environment multiplied by an effective heat transfer coefficient ins . This heat loss term is subtracted in the boundaries of the interconnect of each cell: An environment temperature of 25°C is assumed. The heat transfer coefficient of the ISM isolation is estimated based on a reference operating point specified by the manufacturer of the ISM, which results in a total heat loss of 250 W for the 60 cell ISM. This is in good agreement with the value reported by a system integrator [52].
A prescribed temperature boundary condition is applied to the cell boundaries and the walls representing the ceramic block in the single cell model, since the temperature of the ceramic block was controlled in the experiments with the furnace heating: The system of equations described in this section is implemented in Matlab/Simulink using S-Functions. The system is solved dynamically using implicit Runge-Kutta with a trapezoidal rule as a first and backward differentiation as a second stage (ode23tb).

Stack model validation and evaluation
The ISM model is validated with power curves published by Sunfire/ Staxera for their ISM V3.3 [30,31]. These specify the stack power for different stack currents at reference conditions for three gas compositions: a hydrogen-nitrogen mixture, a CPOX reformate and a SR. An overview of these reference operating conditions and fuel compositions is given in Table 2. The manufacturer advices temperature control trough manipulation of the cathode airflow or inlet temperature. Therefore, a control loop is implemented in the model which adjusts the airflow such that a maximum PEN temperature of 850°C is maintained.
The hydrogen content of the fuel mixture specified varies from 31% in the CPOX reformate to 53% in the SR. In addition, the SR contains 24% steam, while the hydrogen-nitrogen mixture is dry. The hydrogen and steam partial pressure affect the cell voltages and consequently the power curves due to changes in the Nernst potential and the anode exchange current density [41]. However, the dependency of the anode exchange current density is a subject of debate and thus unknown. Therefore, appropriate values are determined for the investigated stack using the power curves reported for three different gas compositions.
The activation polarization model proposed by Costamagna et al. [43] is implemented, which relates exchange current densities to global dependencies on the reactants and products and an Arrhenius temperature dependency. In addition, all parameters for the cathode are adopted, i.e. the cathode pre-exponential factor k ca 0, , reaction order for oxygen , as well as the activation energies for the exchange current densities of both electrodes. However, reaction orders for hydrogen and steam , ranging from 0 to 1 and −1/2 to 1 respectively, are evaluated. The value of the pre-exponential factor is determined for hydrogen-nitrogen operation at a stack current of 27 A for every combination of and . Table 3 presents a selection of evaluated values for Eqs. (24) and (25). Fig. 3 shows the power curves simulated with the ISM model for four exchange current density models compared to the values published by the manufacturer for the three fuel compositions. The model predicts the power curve for hydrogen nitrogen operation with reasonable accuracy, regardless of the formulation of the anode exchange current density. Most formulations perform reasonable well for CPOX reformate as well, but clear differences are observed for SR. The model with proportional dependence on the steam and hydrogen partial pressure predicts higher cell voltages for SR operation than hydrogen nitrogen operation, while the opposite is reported by the ISM manufacturer. The best results are obtained for a square root dependence on only the hydrogen partial pressure (R 2 = 0.99). Therefore, this anode exchange current density formulation is selected. Table 4 compares the simulation results to reference operating parameters specified by the manufacturer. The simulated cathode outlet temperature is only 3.2°C lower than the manufacturer reference, probably because a 0.7 Nl min −1 higher cathode airflow is required to limit the PEN temperature to 850°C in the model. The stack voltage and power are in good agreement as well, demonstrating the ability of the model to simulate the behaviour of a commercial stack with high accuracy. The results confirm that the assumed heat transfer coefficient results in the expected heat loss of 250 W at reference conditions [52].
The advantage of 1D models over the more commonly used lumped parameter models is in the ability to resolve spatial distributions of temperatures, concentrations and reaction rates along the flow direction. Fig. 4 shows the temperature profiles of air, fuel, interconnect and PEN structure along the inactive and active area of stack for the operating conditions shown in Table 4. The temperature difference between the four layers is relatively small along the flow direction, with the exception a small inlet section. The temperatures are substantially higher at the beginning and end of the active area than in the inactive boundaries due to heat losses to the surroundings. The PEN temperature initially drops while heat is transferred from the hot fuel to the cold cathode air, then increases along the active area as waste heat is produced by the electrochemical reaction, and finally drops again due to heat loss to the surroundings. Fig. 5a compares the PEN temperature profile for hydrogen-nitrogen, CPOX and SR operation. The current density distributions and temperature gradients predicted for the active area are shown in Fig. 5b and c respectively. The average PEN temperature is highest for hydrogen-nitrogen operation, closely followed by CPOX. Endothermic cooling from the reforming reaction, assumed to react to chemical equilibrium instantly, causes the PEN temperature to drop sharply at the beginning of the active area for SR.
The current density distribution is relatively homogeneous for the hydrogen-nitrogen fuel mixture, since the effect of a decreasing hydrogen and increasing steam concentration on the Nernst voltage along the flow direction is compensated by the increasing PEN temperature. The presence of steam in fuel and the higher air flow required to cool the stack result in a larger current density variation for CPOX. Similarly, the low inlet temperatures caused by the endothermic DIR reaction induce large current density variations for SR. Therefore, the temperature gradients are substantially higher compared to hydrogen-nitrogen and CPOX fuel mixtures.

Transient simulations
The stack model developed is dynamic and can thus be used to simulate transient operation of the ISM. Although the validation of the stack model is limited to steady-state performance in this study due to the lack of reliable data, the transient capabilities of the model are demonstrated by simulating the start-up behaviour of the ISM for completeness. According to the manufacturer, the stack should be preheated to a temperature 650°C by a start-up burner after which current is drawn to support further heating to the operating temperature. The current should not be increased faster than 2 A min −1 and a stack voltage below 36 V has to be avoided. Fig. 6 shows the simulated ISM start-up behaviour assuming hydrogen-nitrogen operation at reference conditions. Two scenarios are simulated: in the first one the current is ramped to the nominal value of 27 A with the maximum rate specified by the manufacturer, while the current is stepped to the final value in the second scenario. Fig. 6b shows that application of the specified current ramp results in a gradual change of the stack voltage, while a step change of the current causes an instant drop below the minimum specified for safe operation. The relatively high voltage drop for a step change is a result from the high electrolyte resistance at a lower average PEN temperature.
Safe heating of the stack to the desired temperature with the specified current ramp can take up to half an hour, as can be seen in Fig. 6c. A higher current ramp may bring this down to 15 min, but will cause the stack voltage to drop below acceptable values. This may cause large cell-to-cell variations in the stack, resulting in potentially deteriorating local hot spots. In addition, rapid ramping may induce relatively large local temperature gradients, as is evident from the maximum PEN temperature gradient in Fig. 6d. The resulting thermal strain may damage the stack.

Air and fuel compositions
Air

Internal reforming in single cell experiments
Direct internal MSR was experimentally studied by Fan et al. [21] on 10 × 10 cm single ESC2 cells obtained from Kerafol (former H.C. Starck). The reforming data was used to parameterise two different kinetic models, one of the PL and the other of the HW type [29]. Both rate equations were then implemented in a CFD model. Both kinetic models predicted comparable overall MSR rates, but higher temperature gradients were predicted by the HW kinetics due to their nonmonotonic dependency on the steam-to-hydrogen ratio.
The objective of this study is to transfer the kinetic models derived in previous work to the 1D ISM model. A general problem when transferring kinetic mechanisms derived on substrate materials to stack models is the selection of an appropriate value for the frequency factor k 0 in Eqs. (13) and (14). Data fitting typically yields an arbitrary rate constant in mol s −1 , which can be normalised by, for example, the anode volume, active area or the weight or surface area of the catalyst. However, it is difficult to account for geometrical, structural and material variations between different anodes.
It would be more appropriate to obtain an appropriate value for k 0 from stack operation data, but this is difficult as all methane is typically converted within the stack and hence no reforming rates can be deduced. However, in this case the reforming parameters are derived using data obtained from same cells used in the Sunfire ISM V3.3. In theory, this enables direct transfer of the MSR kinetics to the stack model. To demonstrate the validity of this approach, the original experimental conditions are simulated using a single cell version of the 1D model. Table 5 presents an overview of the simulated experimental conditions. The values of k 0 are chosen such that the simulated methane conversion matches the experimentally observed value for gas composition (GC) 1 at 725°C and open circuit conditions. However, the values are found to be within 10% of those obtained in the fitting procedure, in which a simplified isothermal ideal plug flow reactor model was assumed. An overview of the parameters used in the PL and HW rate equation is given in Table 6. Fig. 7a shows temperature profiles along the flow direction calculated with 1D single cell model and those obtained with the 3D CFD model published earlier for GC 1 at 725°C and open circuit conditions [29]. The temperature profiles predicted with the 1D model are qualitatively in good agreement with the CFD model, especially considering that it lacks most geometrical information, such as the shifted position of the fuel in-and outlet tubes. Both models predict higher temperature gradients for the HW kinetics.

Table 3
Parameters used in Eqs. (24) and (25) to calculate the exchange current densities for different global dependencies on the oxygen, hydrogen and steam partial pressures.    Fig. 4. Temperature profiles of air, fuel, PEN structure and interconnect along the flow direction, both in the inactive and active area of the stack, for the conditions presented in Table 4. Fig. 7b shows a comparison of the current-voltage characteristics predicted with the PL and HW rate equations to experimental data. These are obtained with the electrochemical model derived for the ISM model, although the contact resistance is adjusted to account for nonideal contacting in the experimental setup. This uncertainty is difficult to eliminate since the contact resistance depends on the contact between the anode and the current collector, which varies from experiment to experiment, for example due to changes in the compression seals. Nonetheless, it can be seen that the predicted open circuit voltages agree well with the experimental value. Higher voltages are predicted using HW kinetics, probably since the reforming reaction proceeds faster at the entrance according to this reaction model. Fig. 8 shows a comparison between the experimental MSR rates and simulation results for some of the conditions specified in Table 5. Both the PL and HW model show good agreement with the experimental data, although the PL rate equation performs slightly better, in contrast to finding using the CFD model [29]. In general, the validation results give confidence for implementation of the MSR kinetics in the ISM model.

Reforming kinetics in stack modelling
Fig. 9b compares the MSR rates predicted along the active area of the stack with the PL and HW rate expression to the original results assuming chemical equilibrium, for ISM operation at the reference conditions specified in Table 2 for SR fuelling. As expected, the MSR rates are considerably lower for the kinetic models, revealing that the MSR reaction is kinetically limited for typical stack operating conditions.
The PL kinetics predict higher reaction rates than the HW rate equation at the entrance of the active area, in contrast to the findings for the experimental conditions. This is attributed to the relatively low steam-to-hydrogen ratio in the partly pre-reformed fuel compared to the experiment, where the fresh fuel consisted primarily of methane and steam only. Still, even with the slower HW kinetics all methane is reformed within the stack length. Fig. 9a compares the PEN temperature profiles in the active and inactive area of the stack predicted by the PL and HW rate expression to assuming chemical equilibrium. The wider distribution of the endothermic MSR reaction results in a more gradual increase of the PEN temperature. In addition, the somewhat unlikely cold spot predicted assuming chemical equilibrium disappears. As a result, the PEN temperature gradients are smoothened as well, as is shown in Fig. 9c. The maximum PEN temperature gradients are 45.7, 31.9 and 30.2°C cm −1 for chemical equilibrium, PL and HW kinetics respectively. Fig. 10 shows there is little influence of the internal MSR kinetics on the overall power production by the stack. With the HW kinetics the predicted power curve is closest to the one reported by the manufacturer, followed by the PL kinetics and assuming chemical equilibrium. However, the overall maximum deviation between the predicted stack powers is only 5 W. This emphasises that realistic predictions of the DIR rate is primarily important for accurate thermal stress predictions.  (Fig. 4), as well as the current density distribution (Fig. 5b) and PEN temperature gradient in the active area (Fig. 5c) for a reference operating conditions and fuel compositions specified in Table 2 and a stack current of 27 A.

Discussion
The 1D dynamic modelling platform presented in this study enables the simulation of both a single cell experimental setup and a commercial ISM with a single set of equations, changing only geometrical parameters and boundary conditions. Therefore, it enables validation of kinetics derived from single cell experiments and allows for direct transfer to models for stack simulation. Single cell experiments are expected to yield more reliable reforming kinetics than data collection from substrate reactors, while being substantially less complicated and easier to instrument and control than experiments on complete stack assemblies.
Dynamic models of the same stack design have been developed by Kupecki et al. [53,54], Sorce et al. [24] and Greco et al. [55]. Kupecki et al. [54] validated the voltages predicted by their quisi-1D model dynamically with current ramps for two methane containing gas compositions. Greco et al. [55] used a full 1D dynamic model to study faulty states, for example due to reformer malfunctioning. However, it appears that heat transfer in the inactive area of the stack was not included in these models. Moreover, chemical equilibrium was assumed for the internal reforming reaction.
It was shown in this work that assuming chemical equilibrium for the reforming reaction results in unrealistic temperature profile predictions. Thermal stresses can be more accurately estimated using appropriate kinetic models for the MSR reaction on the anode. However, the PL and HW kinetics derived in previous work yielded different reaction rates temperature profiles in the 1D stack model, even though  (Figs. 6b-6d) for a simulated system start-up, either by a step change to the stack current or with the maximum allowable current ramp specified by the ISM manufacturer (Fig. 6a).

Table 5
Overview of the conditions simulated for the experiments on ESC2 cells obtained from Kerafol/H.C. Starck. [21]. The volume flows are specified for atmoshperic pressure and a temperature of 120°C, such that all steam is evaporated.  both were derived from the same experimental data and shown capable to reproduce the experimentally observed conversions in the 1D single cell model. The total number of parameters is higher for the HW than the PL kinetics. However, both the PL and HW kinetics have four free-fitted parameters, since the reactant and product partial pressure dependency is not fitted but based on an intrinsic rate limiting mechanism. The parameters are entirely independent in the PL kinetics, while two preexponential factors and their respective energies are determined for the HW kinetics. Statistically, the PL kinetics showed slightly better agreement with the experimental data, but the HW kinetics ideally contain information on the intrinsic rate determining mechanism and are, therefore, expected to predict the internal reforming kinetics more accurately.
Although PL kinetics do not contain any mechanistic information, they may still yield more realistic temperature profiles than assuming the reforming reaction to be in equilibrium. Assuming proportionality to the methane partial pressure might be acceptable as well if some inaccuracy in the predicted thermal stresses is allowed. However, it is not clear if this holds for more substantially deviating operating conditions, such as changes in the extent of pre-reforming, oxygen-to-carbon ratio in the fuel, inlet temperatures of air and fuel, fuel utilisation and anodic off-gas recirculation. Further study is required to discriminate between different kinetic models and determine the rate limiting step(s) of the reforming reaction on different SOFC anode materials.
The electrochemical reactions in the anode were accounted for using a symmetric Butler-Volmer equation, as it is numerically convenient and Noren et al. [40] showed that errors due to this assumption are small for typical SOFC operation. This result was confirmed in our simulations as well. Global reaction orders for hydrogen and steam were determined in the anode exchange current density to account for the rate limiting hydrogen oxidation kinetics, using power curves from the stack manufacturer for three different fuel compositions.
Bessler et al. [41] pointed out the limitations of using global dependencies, since the hydrogen oxidation reaction is a complicated multi-step process and there is no agreement on exact mechanism. Depending on the rate determining step, the global reaction orders may depend on the hydrogen and steam partial pressures themselves, temperature and anode materials. However, implementing a multi-step heterogeneous reaction mechanism for the hydrogen oxidation reaction may be overly complicated for control-oriented computationally inexpensive dynamic models.
Good results were obtained for a square root dependency of the anode exchange current density on the hydrogen partial pressure, while no evidence of the influence of the steam partial pressure was found. Ignoring these global dependencies or using values reported in literature yielded inaccurate predictions of the ISM power when the fuel composition was changed. Although the global reaction order for hydrogen supports a hydrogen spill-over charge transfer mechanism, a positive steam partial pressure dependency is reported for most charge Fig. 7. Comparison of the temperature profiles predicted with the PL and HW kinetics to those obtained with CFD modelling in previous work [29], and predicted IV curves compared to the experimental data. transfer rate limiting reactions. While hydrogen oxidation kinetics have been studied extensively on Ni-YSZ anodes, this result highlights the necessity to study them on other anodes such as Ni-GDC and wider steam-to-hydrogen ratios as well.
The approach presented was shown to enable simulation of the performance of both single cell experiments as well as a commercial ISM with high accuracies. The model formulation is dynamic, and load transient could be simulated with good accuracy and numerical stability as well as limited computational time. In addition, thermal stresses induced by DIR in a commercial ISM could be predicted based on reforming data obtained in single cell experiments. Therefore, the modelling platform may be a powerful tool for future studies on design, operation and control of SOFC stacks and systems.
The model could only be validated with the available data specified by a single manufacturer in this work. However, the approach may be extended to different stack designs and validated for more operating conditions, such as variations in the extent of pre-reforming, oxygen-tocarbon ratio in the fuel, inlet temperature of air and fuel, fuel utilisation and anode off-gas recirculation. In addition, further validation of the transient predictions is required, as well as further study on the intrinsic kinetics of both the reforming and electrochemical reactions.
The dynamic modelling approach developed in this study can be used to study off-design operation of commercial ISMs and integration with balance of plant components. The models developed can be used to predict the electrochemical performance of the stack, air flow required to maintain a constant temperature as well as thermal stresses induced. In addition, the model may be used to simulate transient operation and develop adequate control logic for integrated systems.

Conclusions
A 1D dynamic SOFC modelling platform was developed in this study to simulate a single cell experimental setup as well as a commercial ISM, changing only geometrical parameters and boundary conditions.  9. PEN temperature profiles in the active and inactive area of the stack (Fig. 9a), and methane steam reforming rates (Fig. 9b) and PEN temperature gradients (Fig. 9c) along the active area for SR operation with different reforming models. The models were used to validate PL and HW MSR kinetics derived in single cell experiments and use these to predict the spatial distribution of the MSR reaction in a commercial ISM equipped with the same electrolyte supported SOFCs with Ni-GDC anodes.
While the HW kinetics predicted higher reforming rates for the operating conditions in the single cell experiments compared to PL kinetics, this trend reverses for the investigated stack operating conditions due to the lower steam-to-hydrogen ratio in the pre-reformed fuel used in stack operation. As a result, the predicted temperature gradients are slightly smaller for the HW kinetics. However, both the PL and HW kinetics predicted more realistic temperature profiles in the ISM model than assuming chemical equilibrium, thus indicating that DIR is kinetically limited for the investigated conditions and giving a more reasonable estimation of the local temperature gradients.
Power curves published by the ISM manufacturer for hydrogen-nitrogen, CPOX reformate and SR fuels were used to determine the global reaction orders for the hydrogen and steam in the anode exchange current density. Good results were obtained with a square root dependency of the anode exchange current density on the hydrogen pressure, while no evidence for an effect of the steam partial pressure was found. The 1D dynamic model was shown to simulate SOFCs in both single cell and stack configuration with good accuracy. The modelling approach can be used to simulate off-design conditions and transient operation of commercial ISMs, study system integration and develop adequate energy management and control strategies of SOFC systems. However, further study on the rate limiting steps in the reforming and electrochemical oxidation reaction is required. In addition, the transient predictions of the model need further validation.