Coupled Transport Effects in Solid Oxide Fuel Cell Modeling

With its outstanding performance characteristics, the SOFC represents a promising technology for integration into the current energy supply system. For cell development and optimization, a reliable quantitative description of the transport mechanisms and the resulting losses are relevant. The local transport processes are calculated by a 1D model based on the non-equilibrium thermodynamics (NET). The focus of this study is the mass transport in the gas diffusion layers (GDL), which was described as simplified by Fick’s law in a previously developed model. This is first replaced by the Dusty-Gas model (DGM) and then by the thermal diffusion (Soret effect) approach. The validation of the model was performed by measuring U,j-characteristics resulting in a maximum deviation of experimental to simulated cell voltage to up to 0.93%. It is shown that, under the prevailing temperature, gradients the Soret effect can be neglected, but the extension to the DGM has to be considered. The temperature and heat flow curves illustrate the relevance of the Peltier effects. At T=1123.15 K and j=8000 A/m2, 64.44% of the total losses occur in the electrolyte. The exergetic efficiency for this operating point is 0.42. Since lower entropy production rates can be assumed in the GDL, the primary need is to investigate alternative electrolyte materials.


Introduction
The demand for electrical energy has been steadily increasing in recent years, and, along with it, the importance of environmentally friendly energy converters. Electrochemical systems, such as fuel cells, are such alternative technologies for directly converting the chemical internal energy of a fuel into electrical energy. Hydrogen as an energy carrier becomes inexhaustible due to its simple production using solar-based carbon free energy sources and enables its use without direct CO 2 emissions, making it environmentally friendly, also with regard to the Paris Climate Agreement.
The solid oxide fuel cell (SOFC) is a promising technology due to its superior performance indicators. These include the integration into the current energy supply system through the use of carbon-containing fuels, such as natural gas, the variability of the purity of the hydrogen required, and high efficiency due to the high operating temperatures. According to the current state of development, the electrical efficiency of SOFCs reaches values of up to 60-65% [1]. For the development of improved SOFCs, it is of great importance to be able to simulate the processes and the operating behavior taking place inside of the cells with the help of models. This is especially true as an experimental assessment of the ongoing processes inside the cell is prohibitive because of the temperature level around 1000 K. Many models already exist in the literature. These range from modeling individual layers within one cell to modeling the SOFC as a stack system. Three-dimensional models were set up in the work of Anderson et al. [2], Mauro et al. [3], Yakabe et al. [4], and Peksen et al. [5], in which the SOFCs are modeled based on the 3D finite element method. The momentum transport is described in these works with the help of the Navier-Stokes equations. In the work of Peksen et al., the SOFC is modeled, including peripheral components, such as the heat exchangers and the recirculation loops. Individual layers were investigated, for example, by Xu and Dang [6], Niu et al. [7], and Joos [8]. They model the porous electrodes of the SOFC in terms of their exact microstructure and investigate the relationship between the microstructure and the electrochemical reactions at the three-phase boundaries. In all these models, the transport processes are described using the classical empirical equations for the transport of heat, charges, molecules, and momentum. Why these empirical equations may not be sufficient for an exact description of the transport mechanisms in an SOFC is addressed in this work.
The performance and lifetime in an SOFC are determined by large temporal and spatial temperature gradients. Zeng et al. [9] shows an extensive literature review on heat transfer in SOFC stacks and derived thermal management methods. The study shows that excessive temperature gradients in SOFCs can lead to delamination and cracks in the electrolyte and electrodes. Depending on the operating point, temperature gradients in a planar co-flow SOFC along the electrolyte in the main flow stream of up to 30 K/cm can be reached. To avoid damage in the SOFC, a temperature gradient of 10 K/cm should not be exceeded [10]. Such temperature gradients can be reduced in the system by means of an effective heat management.
Temperature gradients are also caused by local heat generation and by flow field arrangement. Two types of heat generation are relevant in an SOFC. Reversible heat describes the reaction entropy of the system as a function of temperature, while irreversible heat is caused by dissipated energy. Dissipated energy arises due to ohmic resistance and overvoltages at the electrodes. The heat generation due to an electric current is described by means of the first Joule's law, which is why the irreversible heat is also called Joule's heat. In contrast, the electrochemical reaction at the three-phase interface leads to a thermoelectric effect in which reversible heat is released or absorbed at each electrode, also called Peltier heat.
In addition to a temperature gradient or a gradient in the electric potential, a gradient in the chemical potential of a species can also lead to a heat flow, called the Dufour heat. The influence of this coupling of multiple transport mechanisms caused by gradients of different intensive variables on the heat flux in a PEMFC is investigated in Reference [11], where individual effects on heat generation (Fourier, Peltier, and Dufour) are examined. As expected, at typical current densities, there is a net heat flux pointing out of the anode and cathode. An increasing current density leads to a higher potential difference and a higher electro-osmotic effect, which increases the Peltier and the Dufour effect. However, above a current density of 9700 A/m 2 , Fourier heat dominates due to the temperature rise of the cathode surface caused by the irreversible activation overvoltage. Valadez Huerta [12] investigates, among other things, the contribution of the Peltier effect to the heat fluxes in an SOFC. A significant contribution is shown to come from the electrolyte, where the heat flux can flow from a lower to a higher temperature within the electrolyte due to the coupling effects mentioned above.
At the same time, the temperature gradients which develop affect the charge flow (Seebeck effect) and the material flow (Soret effect). The Seebeck effect describes a charge transport caused by a temperature gradient. Thus, due to this thermoelectric effect, a voltage can be measured between two different conductors when a temperature gradient is imposed on them. The magnitude of this voltage is described by the Seebeck coefficient, which depends on the material and the temperature. In the works of Kjelstrup et al. [13,14], a temperature gradient is applied to an electrochemical cell in order to measure the Seebeck coefficient. From this, the entropy of the oxygen ions transported in the electrolyte can be determined, which can also be used to calculate the Peltier coefficient.
Fick's law, the Stefan-Maxwell diffusion, or the Dusty-Gas model (DGM) are widely used to describe mass transport within a cell. In Reference [12], the coupling of the mass flow with other flows is neglected. Using Fick's law, the concentrations of the individual components are calculated using multi-component diffusion coefficients from the work of Costamagna et al. [15]. In the work of Suwanwarangkul et al., the approaches for describing the mass transport by means of Fick's law, Stefan-Maxwell diffusion, and the DGM for determining the concentration overvoltage at the anode are examined. The modeling with the DGM achieves the best results, whereas the approaches using Fick's law and Stefan-Maxwell diffusion are good approximations, especially for small electrical current densities and large pore diameters [16].
In Reference [11], the modeling of a PEMFC using non-equilibrium thermodynamics (NET) takes into account the Soret effect, whereby a material flow is driven by a temperature gradient. In operating ranges of high electric current densities, a significant influence of the Soret effect is observed in the membrane of the PEMFC. This is due to the strongly increasing temperature gradients in the membrane with increasing electric current density.
Considering these results, the coupling effects for a reliable description of the transport mechanisms in fuel cells cannot be neglected. It is not sufficient to address these coupling effects individually and simply add them up. It is very important to take an integrated approach as will be shown below. An integrated coupling of several transport processes can be based on theory of non-equilibrium thermodynamics (NET), which has only been considered in a few models so far. The first known model based on NET is the 1D model by Kjelstrup and Bedeaux [17], which can be used to calculate the potential field and the temperature curve of the cell. Within the NET theory, a generalized flux J i results from the linear combination of all occurring forces X j [18]: where the phenomenological coefficients L ij , also called conductances, are the connection between forces and fluxes. Forces are gradients of intensive variables. From the phenomenological equations, a matrix with the N 2 conductances L ij necessary for modeling is obtained. By Onsager's reciprocity relation L ij = L ji , the problem reduces to the upper and lower triangular matrices of this L ij matrix, including the coefficients L ii of the main diagonal [19].
For the optimization of a technical system, it is of great importance to understand these local transport mechanisms and their driving forces, as they all dissipate energy, i.e., create entropy. With the help of the local entropy production rate, loss mechanisms can be identified and quantified as exergy losses. In terms of NET, the local volume specific entropy production rateσ i due to the transport process is generally calculated by multiplying a flux J i by the respective corresponding force X i [20]: thus allowing the second law of thermodynamics to be integrated into the model approach. However, if several fluxes are present simultaneously in a process, the total local entropy production is calculated by summing all products of occurring fluxes J i with the respective corresponding force X i [20]:σ Sauermoser et al. concluded that, in a PEMFC, the total local entropy production rate is highest at the cathode, which is explained by the potential profile. In contrast, Valadez et al. break down the total local entropy production rate of an SOFC into its individual contributions. In the gas diffusion layers (GDL), the diffusion process provides the largest contribution to the local entropy production rate, with losses at the cathode GDL being higher than losses at the anode GDL. In the electrolyte, the local entropy production rate is dominated by ionic conduction. If the Peltier effect is neglected, the GDLs show reduced values for the heat flux and potential contribution. As a result, a change of direction of the heat flux in the electrolyte happens. In both cases, the highest local entropy production rate results in the electrolyte, with up to 73%.
Our work aims at providing a consistent model approach for a single solid oxide fuel cell based on the NET theory. The 1D NET model by Valadez Huerta [12] provides the starting model for our approach. However, the mass transport described by Fick's law is extended by transport equations that consider the coupling of several transport mechanisms. In addition, the phenomenological equations of the electrolyte are transformed to model the electrolyte with measurable coefficients. With the localization of the individual loss mechanisms, the system is then evaluated exegetically.

The Thermodynamic System-Solid Oxide Fuel Cell (SOFC)
For modeling purposes, the SOFC is divided into five different layers; see Figure 1. The membrane electrode assembly comprises a homogeneous bulk phase, which is separated by two surfaces. The boundaries to these two surfaces are two GDLs. Via the flow field plates, which are not considered in this model, the reactants are homogeneously distributed to the anode and cathode GDL. The surfaces also form interfaces between the three homogeneous bulk phases. The three homogeneous bulk phases, i.e., electrolyte membrane (e), anode GDL (a), and cathode GDL (c), each have a thickness of ∆y i . The oxygen reduction reaction (ORR) takes place on the surface of the cathode catalyst layer, and, at the surface of the anode catalyst layer, the negative oxygen ions react, together, with hydrogen in the hydrogen oxidation reaction (HOR): Atomic vacancies in the electrolyte allow diffusion of the negative oxygen ions. In this study, 8 mol% yttria-stabilized zirconia (8YSZ) is investigated as the electrolyte material. The electron transport is ensured via wires that are connected to the catalytic layers. Nickel (Ni) is used at the anode, and platinum (Pt) at the cathode. The cell can be described with the following standard electrochemical notation: The upper script describes the anode (a) or the cathode (c). We define y as the spatial coordinate perpendicular to the active surface of the SOFC. The relevant molar fluxes J i are defined according to their flow direction, and the heat fluxes J q are defined as positive in the direction of the y axis. The flux J a q (y) describes the heat through the homogeneous anode side GDL, while the heat flux occurring at the reaction layer is defined as J a q (∆y a ) = J d,a q .
The heat flux J e q (∆y a ) = J e,a q , which has been embossed by the half-cell reaction, flows from the anode reaction layer into the electrolyte. The definition for the cathode side is analogous. The direction of the electric current density j is assumed to be the technical direction of the current against the electron flow. Throughout the entire simulations, an equimolar mixture of hydrogen and water on the anode side and air on the cathode side is supplied. Ideal gas mixtures are assumed. The thermodynamic reference state point is set to be at T = 1123.15 K and p = 1 bar.

Mass and Energy Balance Equations for the Gas Diffusion Layers (GDL)
In the steady state, some component flow densities and the electric current density are in a fixed relationship to each other. The following relationships result from a steady-state component and charge balance: J i k describes the flux of each component, j is the electric current density, and F is the Faraday constant. The steady-state energy balance for the GDL of the anode and cathode side of an in y-direction infinitesimal volume element reads for the positive y direction: Heat J i q , electric power j · dφ, and the gaseous components i with a molar enthalpy H m,i are supplied to or removed from the GDL. As ideal gas behavior is assumed in our approach, and no enthalpy of mixing has to be considered. If the relationship from the component balance is used, the energy balances result in: Cathode: 0 = dJ c q dy + j dφ dy The molar isobaric heat capacities C iG m,p,i (T) are calculated using a power series approach according to Kabelac et al. [21].

Transport Equations for the Homogeneous Phase
Based on Equation (1), the general transport equations for the homogeneous phases are [20]: with the local thermodynamic temperature T = T(y), the chemical potential µ k of component k, and the chemical potential µ j of component j = k. Unfortunately, the phenomenological coefficients L ij are not well known at this point, and only some measurements of these special coefficients have been performed. They differ from the well-known coefficients from the empirical transport equations, but some helpful relations can be established using the following definition for an SOFC [20]: Diffusion coefficient : Electrical resistance : Peltier coefficient : Thermal diffusion coefficient : Very much attention has to be paid to the state variables which have to be constant for these relations. For an SOFC, simplified transport equations result if ideal gas mixture assumptions are made: The electrical resistance r i is calculated via an Arrhenius approach according to Hajimolana et al. using the pre-exponential factor r i,0 [22]: To calculate the Peltier coefficient, the entropy flow within the homogeneous phase i is derived according to NET, as follows, which, of course, also includes the local entropy production [20,23]: With the correlation from Equation (15), obeying the conditions dT = 0, dµ i,T = 0 and the assumption that no electroosmosis takes place, the following equation holds true for the Peltier coefficient: with S m,k as the transported molar entropy of all n − 1 components. The n-th component is chosen as the reference of frame. In this case, this is the positive ion lattice of the electrolyte as an unmoved reference wall, so that, for the anode and cathode reaction layer, the following equations result: The molar entropy S m,k of component k at temperature T is calculated by the entropic equation of state of ideal gases with the approach according to Kabelac et al. [21] for the molar isobaric heat capacities. The entropies of electrons in platinum and nickel are taken from Reference [24,25]. Approaches for determining the diffusion coefficient D k and the thermal diffusion coefficient D T are presented in the following sections. The local entropy production rate of a homogeneous phase can now be estimated from Equation (2) as follows:σ

Modeling of the Gas Diffusion Layers (GDL)
Fick's law can be derived by simplifying assumptions from the more general Stefan-Maxwell diffusion and, thus, represents the simplest of the diffusion modeling. For the description and rating of a more complete and real mass transport modeling in the GDLs, Stefan-Maxwell diffusion (1) is described first and then overlaid with Knudsen diffusion (2), convection (3), and thermal diffusion (4). For a comparison of all these approaches, the mass transport is additionally modeled with Fick's law. The Fick's diffusion coefficients are calculated from the Knudsen diffusion coefficients, and the binary diffusion coefficients using the Bosanquet formula [16,26]: Krishna et al. [26] state that, for molecules, such as H 2 , the Bosanquet formula is applicable to a wide range of pore diameters with reasonable accuracy. Since porous structures are assumed, some parameters, such as the permeability B i 0 and the diffusion coefficients D, must be corrected by the porosity i , the tortuosity τ i , and the mean pore diameter d i p . The NET approach according to Equation (10), respectively, Equation (18), represents another possibility for modeling the GDLs, and it is subsequently compared directly with the Dusty-Gas model (DGM) (1)-(4) described previously. (1) If a mixture is considered, a diffusive movement relative of the individual components can occur. This relative movement leads to collisions between the different molecules. The driving force, caused by a gradient in the partial pressure of a component k, is in equilibrium with the frictional force, which is caused by the intermolecular collisions [27]. The equilibrium between the driving force and the frictional force in an ideal gas mixture of N components in a porous structure can be described by the following equation:

Mass Transport by Stefan-Maxwell Diffusion
For the binary diffusion coefficients, we have Ð eff kj = Ð eff jk . For the calculation of the effective binary diffusion coefficients, the approach of Fuller et al. is used [28]: Here, v are the dimensionless diffusion volumes of the two components in a binary mixture, which are tabulated, e.g., in the VDI Heat Atlas [29]. (2) In narrow pores, in which the mean free path length of the molecules is considerably greater than the pore diameter, there are increased collisions of the molecules with the pore wall [30]. Because of these collisions, the following gradient in partial pressure result for the individual component k according to Reference [27]:

Extension by Knudsen Diffusion
For the mass flow density of nitrogen, J N 2 = 0 applies in a steady state operation, as this inert gas does not react at the reaction layers. The effective Knudsen diffusion coefficients D Kn,eff k of the individual components are calculated according to Krishna and Wesselingh [27] via the kinetic gas theory: To superimpose the Stefan-Maxwell diffusion and the Knudsen diffusion, the partial pressure gradients of the two transport mechanisms may be added [27,30,31]. The idea follows exclusively from the additivity of the momentum transfer, without any theoretical proof [31]. The following transport equation results in a steady state operation: 2.2.6. Extension by Convection (3) In convection, the gas mixture moves as a continuum and is driven by a gradient in the total pressure. Typically the adhesion condition applies to the walls of the flow channel from a macroscopic view, which means that the velocity of the viscous flow is zero directly at the wall interface. The convective mass flow density of the respective component can be calculated by Darcy's law, with the permeability B i 0 and the dynamic viscosity of the mixture η [27]: The permeability B i 0 depends on the geometry of the pores and must also be corrected due to the porous structure. For the effective permeability within circular open pores, Equation (33) applies [30]: The total mass flow density, which occurs in a steady state operation, results from the summation of the diffusive and convective component [27,30,31]: This results in the following transport equation, which, in addition to Knudsen and Stefan-Maxwell diffusion, also takes convection into account: dp k dy The dynamic viscosities of the individual gases η iG k are calculated using Lucas et al.'s correlation equation recommended in the VDI Heat Atlas, which is permissible due to the low operating pressure in the SOFC [29]. All quantities required for the calculation of the dynamic viscosities of the pure substances can be taken from Reference [29]. Subsequently, the dynamic viscosities of the gas mixtures η are calculated for both the anode and the cathode on the basis of the mass fractions and the dynamic viscosities of the pure substances via the mixing rule according to Wilke [32]: with (4) Thermal diffusion, also called the Soret effect, is a coupled process, in which a flow of substances is caused by a temperature gradient. If the mixture is subject to a temperature gradient, gradients are established in the concentration of the individual components [33]. Heavy molecules arrange themselves more in the colder regions, and light molecules in the warmer regions. Thermal diffusion belongs to multi-component diffusion. Therefore, a version of the Stefan-Maxwell diffusion according to Krishna and Wesselingh, extended by the thermal diffusion rate, is used to consider thermal diffusion [27]:

Extension through Thermal Diffusion
with The extended velocity w T k is obtained by the summation of the diffusion velocity w diff k with the velocity due to thermal diffusion, which can be calculated by the thermal diffusion coefficient D T kj . A positive thermal diffusion coefficient means that the component moves toward colder regions, and a negative coefficient means that the component moves toward warmer regions [27]. For the thermodiffusion coefficient, we have [31]: The thermal diffusion coefficients are often expressed by the thermal diffusion rate k T kj or the thermal diffusion factor α T kj [34]: Taking thermal diffusion into account, the general transport equation within the GDL results: dp k dy From the equation, it can be seen that thermodiffusion can both improve or dampen mass transfer. This depends on the signs of the thermal diffusion factor and the temperature gradient. If a component moves against the expected transport direction as a result of the thermal diffusion, the thermal diffusion acts as an additional transport resistance that must be overcome by a higher partial pressure gradient.
The thermal diffusion factor α T kj can be calculated by means of the kinetic theory of gases, which relates the macroscopic transport coefficients to the intermolecular interactions of gases [35]. A model that describes the intermolecular potential as a function of the distance between the molecular centers r * is the Lennard-Jones (12-6) potential φ LJ . Data for the maximum attraction between molecules k and for the characteristic Lennard-Jones length σ k are taken from Reference [36]. Widely used calculation approaches for the thermal diffusion factor α T kj , which are based on the kinetic gas theory, are the first two approximations according to Chapman and Cowling and the first approximation according to Kihara. The first approximation according to Kihara provides better results than the approaches of Chapman and Cowling [37], which is why it is used in this work [38]: Here, S k and Q k are functions of the molecular masses, the Lennard-Jones parameters, and the collision integrals which account for the different interactions between the individual molecules in a binary mixture. Additional collision integrals are summarized in C * 12 . These multi-dimensional collision integrals were calculated and tabulated in Reference [35] for different potential functions, depending on the temperature. By convention, the index 1 stands for the molecule with the greater mass.

Mass Transfer Approach According to NET
Through comparison of the coefficients of Equations (10) and (18), it follows that, for the phenomenological coefficients: and The diffusion coefficient D k can be determined using Equation (26). The thermal diffusion coefficient D T is calculated via the approach of Equation (40).

Modeling of the Reaction Layers
In this model, reactions are assumed to take place on the surfaces of the electrodes only. A homogeneous temperature profile is additionally assumed at the electrode surfaces. However, jumps in the y-direction occur in the course of the electrical potential and the heat flux density. The heat flux density is calculated via an entropy balance, including the reaction entropy. The material current densities are expressed by the electric current density. For the oxygen ion flux, with The entropy of the oxygen ion S O 2− is calculated using the approach of Kjelstrup and Tomii [14]. The entropy production ratesṡ i irr result from the activation overvoltages at the respective electrode. Butler-Volmer kinetics is used to calculate the activation overvoltages: A steady-state charge balance yields j = 2 · F · ξ a HOR for the anode and j = −2 · F · ξ c ORR for the cathode. The exchange current densities are determined by an Arrhenius approach according to Yonekura et al. [39]: where the pre-exponential factors γ a and γ c . The activation energies E i A,j 0 have been determined empirically by electrochemical impedance spectroscopy in Reference [12]. Using the activation overvoltages, the entropy production rates are finally calculated as: The potential curve at the reaction boundary layer is calculated as: Thus, the potential immediately behind the reaction boundary layer is calculated from the potential immediately before the reaction boundary layer, the electrode potentials in no-load operation, and the activation overvoltages. As the equilibrium electrode potentials are calculated with the reactant concentrations in the reaction zones, a determination of the concentration overvoltage η i,con is not necessary. The electrode potentials in an open circuit situation at the individual electrodes are obtained by applying the equilibrium condition of an electrochemical reaction and correspond to ∆φ i 0 = (μ O 2− − 2μ e − )/2 · F [12]. The electrode potentials are given by the following equations:

Modeling of the Electrolyte
In an SOFC electrolyte, the diffusion of the oxygen ions leads not only to a mass transport but also to a proportional charge transport. This relationship is explained by the definition of the electrochemical Guggenheim [40], where µ O 2− is the chemical potential of the anion, z * O 2− = 2 is the charge number of the anion, and ψ is the electrostatic potential. Taking into account the charge conservation law and the assumptions that no cation diffusion occurs in YSZ at 1300 K and polarization effects can become negligible, we obtain [12]: with Equation (56) and the relation according to Equations (12)-(15), the following transport equations apply to the electrolyte [12,41]: with φ = φ(y) as the electric potential. By the relation from Equation (56), two of three transport equations are independent of each other, which completely describe the coupled transport of heat and oxygen ions. The transport equations for the electrolyte finally result as a function of the empirical coefficients: For the Peltier coefficient in the electrolyte, Equation (22) applies: The gradient of the heat flux density is determined from an energy balance in a time independent state:

Exergy Analysis
For a full evaluation of the efficiency of the system, the exergetic efficiency ζ is used. ExergyĖx is a thermodynamic state variable, and it describes the part of energy that can be converted into any other form of energy. An energy conversion dissipates energy; so, some exergyĖx is converted to anergyḂ. This fraction is called exergy lossĖx L . In this dissipative process, anergy cannot be converted into any other form of energy. For the exergetic efficiency of a fuel cell, the fraction of the useful outgoing exergy to the incoming exergy is: withĖx i ph as physical exergy of a material flow [42]: andĖx i ch as chemical exergy of a material flow [42]: The equations are valid for ideal gas mixtures. The chemical exergy of a substance Ex θ m,k (T 0 ) at T 0 = 298.15 K can be taken from Reference [43].

Simulation
The model equations derived in the previous section are used to calculate onedimensional profiles of partial pressures, temperature, heat flux, electric potential, and local entropy productions rate. MATLAB ® version R2020a software is used to perform these calculations. The operating conditions of the cell are fixed, which include the operating temperature T, the operating pressure p, the composition of the inlet gases, and the electric current density j. These operating conditions also determine the boundary conditions of the cell. The cell temperature agrees with the operating temperature, both when y = 0 and when y = ∆y a + ∆y e + ∆y c . At the anode side of the cell, the electric potential φ (y = 0) is set to zero. The other boundary conditions are determined by the partial pressures of the inlet gases. As start values for the iteration, the heat flux density J q (y = 0) and the partial pressures p R O 2 and p R N 2 of the cathode gases in the reaction zone are estimated first. The differential equation systems for the anode GDL, the electrolyte, and the cathode GDL, which result from the transport and balance equations, are solved step by step using the Runge-Kutta method. The profiles at the reaction layers are calculated using the Equations from Section 2.2.9. If the boundary conditions are not met after the first calculation round, the starting values are varied, and the cell is calculated again. The starting values are then adjusted using Newton's method to start the next iteration step.

General Parameters
The parameters for a specific example case are summarized in Table 1. The dimensions of the modeled SOFC single cell correspond to a type KeraCell III from the company Kerafol GmbH (Eschenbach i. d. Opf., Germany) [44]. The validation is carried out using the known U, j characteristic curve of the cell under consideration. The U, j-characteristic curve determined experimentally by the manufacturer can be taken from the data sheet. Furthermore, a U, j-characteristic curve was recorded at the Institute with an existing SOFC/SOEC test bench (Evaluator C1000-HT) from HORIBA FuelCon (Magdeburg-Barleben, Germany) under the same conditions.
For an electrolyte of YSZ08, many references exist in which the thermal conductivity or the ionic conductivity were determined experimentally. In the publications by Schlichting et al. [45], Sasaki et al. [46], and Radovic et al. [47], measured thermal conductivities are presented as a function of temperature. Here, the thermal diffusivity a is measured using laser flash analysis. The density of the test specimen is determined by Schlichting et al. using Archimedes' principle. In the work of Sasaki et al., the mass and volume of the test specimen are measured in order to subsequently calculate the density. Radovic et al. do not make any statements about the measurement of the density. The specific heat capacity C of the ceramic is measured in all three papers using dynamic differential calorimetry. Sasaki et al. determine the thermal diffusivity of a non-porous ceramic by linear extrapolation of the values of a porous ceramic. The ionic conductivity of YSZ08 was measured experimentally in the work of Gibson and Irvine [48], Antunes et al. [49], and Artemov et al. [50]. In these papers, the impedance spectrum of the ceramic is measured at different temperatures using impedance spectroscopy to subsequently determine the ionic conductivity. Similar results are obtained in all three papers. Due to the larger data range, the data of Gibson and Irvine are chosen for our modeling.

Parameters for Mass Transport
For a detailed investigation of the influence of different mass transfer equations, the mass flow of all individual components are plotted as a function of the current density in the anode and cathode GDL. The required parameters and their origin can be found in Table 2.  With increasing current density, the deviation grows until it reaches a value of 0.78% at j = 4000 A/m 2 . Reasons for this are the calculation of the losses in the reaction layers and/or in the electrolyte. The dominance of the ohmic losses is due to the limited ionic conductivity of solid oxide ceramics and the strong dependence of the ohmic losses on the electric current density. From the calculation of the overvoltages at the electrodes due to the reaction kinetics, high exchange current densities occur due to the high operating temperature. However, in contrast to the ohmic losses, there is no linear relationship to the increasing electric current density. For this reason, it is assumed that slightly increased ohmic losses are assumed in the model.  The deviation of the simulated OCV is 0.84% compared to the experimentally determined value. This can be explained by the leakages occurring in the experiments and slightly varying gas concentrations. With increasing current density, the deviations also increase. At a current density of j = 4000 A/m 2 , there is a deviation of 0.93%. The deviations can be explained by several factors. These include activation and concentration overvoltages at the electrodes, charge transport losses across the electrolyte, and ohmic resistances in the electrical contacts. The deviations of the model given are within a reasonable range, such that further investigations will be continued with the parameters present.

Partial Pressures
In Figure 3, the partial pressures of all components in the reaction layers are shown as a function of the electric current density at T = 1123.15 K. The calculated various diffusion coefficients are summarized in Table 3. The Stefan-Maxwell diffusion and Fick's law graphs are described as pure mass transport for better comparability. All further graphs include the previously implemented mechanism.  With increasing electric current density, the chemical reaction rates also increase, resulting in a decrease in the partial pressure of the hydrogen on the anode side and an increase in the partial pressure of water. Due to the same effective binary diffusion coefficient, the curves of the pure Stefan-Maxwell diffusion show the same slopes in terms of magnitude. Furthermore, the graphs for the Knudsen diffusion-Fick's law and convection-thermal diffusion overlap. Since Knudsen diffusion additionally takes into account the collisions of the molecules with the pore wall, there is a significant difference to pure Stefan-Maxwell diffusion. The overlaid Knudsen diffusion shows hardly any noticeable differences in the partial pressures compared to Fick's law, although Fick's law does not take into account the force of friction caused by the intermolecular collisions. The mass flow densities of the hydrogen and the water are the same due to the equimolar composition, which means that no differences are apparent in the superimposed Knudsen diffusion and Fick's law. However, the different components are affected differently by Knudsen diffusion. On the one hand, this is due to the different molecular sizes of the components, which results in the largest effective Knudsen diffusion coefficient for hydrogen. Thus, the effect of Knudsen diffusion on the mass transport of hydrogen is the smallest. As the partial pressure gradients of the individual components are not equal on both the anode and cathode sides due to the different effect of Knudsen diffusion, a gradient is established in the total pressure, which leads to some additional convection of the gas mixtures. However, it can be seen that the effect of convection is small in comparison to the other transport mechanisms, although all gases have low dynamic viscosities. The reason for this lies in the low gradients in the total pressure. In the anode reaction layer, there is a higher total pressure than in the gas channel, which is why the gas mixture moves in the direction of the gas channel as a result of the convection, thus impeding the mass transport of the hydrogen and improves the mass transport of the water. Thermal diffusion shows no influence on mass transfer under the operating conditions considered. The reason for this is that the possible prevailing temperature gradients are too small to give thermal diffusion a significant influence. The temperature gradients in the anode, electrolyte, and cathode are explained in more detail in the following sections.  The partial pressure of the nitrogen remains unchanged in the steady state due to J N 2 = 0. For this reason, the partial pressure curves of the Stefan-Maxwell and Knudsen diffusion also overlap. Since the total pressure in the gas channel is higher in the cathode than in the reaction layer, convection improves the oxygen transport. This leads to a larger gradient in the partial pressure of nitrogen. These partial pressure gradients of nitrogen are necessary to counteract the drag force due to the molecular collisions with the oxygen and the convection, thus ensuring J N 2 = 0. The influence of Stefan-Maxwell diffusion on the cathode gases is greater, since a clearly smaller effective binary diffusion coefficient results for the cathode-side gas mixture. As oxygen is the largest of the molecules present, it has the lowest effective Knudsen diffusion coefficient. However, the effect of the superimposed Knudsen diffusion is smaller than for water due to the lower mass flow density of oxygen. The improved mass transfer of oxygen due to convection is very small, as the gas mixture consists mostly of nitrogen on the cathode side. The mass transport via Fick's law shows the greatest changes in the gradient of the partial pressure for oxygen. By taking into account the mole fraction of oxygen in the Stefan-Maxwell diffusion, the mass transport is improved by a significant amount. On the cathode side, thermal diffusion also has no influence on mass transport under the operating conditions considered here. Figure 4 shows the partial pressure profiles of the gas components, determined via the thermal diffusion approach according to Section 2.2.7, and the TiP approach according to Section 2.2.8. In the anode, the partial pressures of hydrogen and water determined by the TiP approach behave as the superimposed Knudsen diffusion, as the thermal diffusion approach includes the extension to convection. Due to the prevailing low temperature gradients, the thermodiffusion in the TiP approach according to Equation (10) does not take influence on the partial pressure curves. Even imposing a temperature gradient of ∆T = 5 K over the electrolyte thickness ∆y e shows a change in the range of 10 − 3 kPa only. Only by imposing ∆T = 50 K over ∆y e can two essential facts be clarified. First, the mass transfer for both components is worsened by imposing of a temperature gradient, i.e., the effect of thermal diffusion. Secondly, the consideration of convection is highlighted, which worsens the mass transfer of hydrogen and improves that of water. The differences between the two approaches can be explained using Equations (18) and (40). In Equation (40), the thermal diffusion coefficient is divided by the binary diffusion coefficient, whereas, in Equation (18), the thermal diffusion coefficient is divided by the binary and Knudsen diffusion coefficients after rearranging the equation. As the Knudsen diffusion coefficient of hydrogen is higher than that of water, the term has a smaller effect on hydrogen, so that smaller differences are seen in the partial pressure of hydrogen between the two approaches. Due to the very small thermal diffusion coefficients of oxygen and nitrogen, thermal diffusion does not affect mass transfer despite the imposition of a temperature gradient (cf. Figure 4). The similar molar masses of oxygen and nitrogen result in small thermodiffusion factors, which means that mass transfer of a mixture of these two components due to a temperature gradient has no effect.   Figure 5 shows the simulated heat flux density curves and the temperature fields at T = 1123.15 K for different electric current densities as a function of the spatial coordinate y. It can be seen that the heat flow J q always has a negative sign, i.e., heat flows from the cathode toward the anode. The jumps in the course of the heat flow at the reaction layers follow, on the one hand, from the irreversible activation overvoltages. In addition, at the cathode, heat is released due to the negative reaction entropy (∆ R S c m (1123.15 K) = −80.47 J/mol K), whereas the anode consumes heat due to the positive reaction entropy (∆ R S a m (1123.15 K) = 21.65 J/mol K). This leads to a big amount of heat being released at the cathode due to these two effects. Net heat is absorbed at the anode because, here, the heat consumption due to the reaction entropy exceeds the heat released due to the irreversibilities. Thus, in the steady state, there is an increased temperature at the cathode and a decreased temperature at the anode. Within the electrolyte, the heat generated due to irreversible ion transport can be recognized by heat flux densities increasing in magnitude in the direction of the charge transport. Considering the temperature curve along the spatial coordinate, the relevance of the Peltier effects in the cell becomes clear. At the GDL of the cathode, the heat flow flows from a lower to a higher temperature accordingly. At the GDL of the anode, the Peltier effect shows a smaller impact due to the smaller coefficient (π a (1123.15 K) = −0.433 J/C, π c (1123.15 K) = −0.742 J/C). Furthermore, it can be observed that both the heat flux density and the temperature gradients become larger with increasing electric current density. This can be explained by increasing irreversibilities at higher electric current densities. The case of equal temperatures at the anode and cathode sides of the cell is motivated by the cooling gas within the bipolar plates. Here, the heat transfer on the anode side is much more intense than on the cathode side. If the same heat transfer coefficient is assumed, the temperature at the anode side will be higher. Thus, a new case is studied.
In Figure 6, the resulting temperature and heat flux density curves are shown for an imposed temperature difference of 5 K between the anode and cathode. An almost linear temperature profile is obtained, wh ich is due to the small temperature differences during normal operation. The heat flux flows continuously from the anode toward the cathode, as the impact of the Peltier effects on the total heat flux decreases, and the heat flux is primarily transported by the temperature gradient. However, the imposed temperature difference of 5 K only manifests itself in the course of the temperature and the heat flux density and has no significant influence on the overall operating behavior of the cell.  Figure 6. Simulated profile of (a) heat flux density and (b) temperature along spatial coordinate y for a imposed temperature gradient ∆T = 5 K between cathode and anode for j = 8000 A/m 2 . Figure 7 shows the simulated electric potential along the spatial coordinate as a function of the electric current density for T(y = 0) = T(y = ∆y a + ∆y e + ∆y c ). In the course of the electrical potential, jumps can be observed at the reaction layers. These result from the equilibrium electrode potentials of the electrodes and the overvoltages that reduce the potential differences at the reaction layers. The electrical potential is greatest in the electrolyte due to the oxygen ion transport that occurs there. With increasing electric current density, greater losses in the electric potential can be observed there, which follow from the irreversible ion transport. The voltage losses due to the charge transport in the GDLs can hardly be observed due to the low electrical resistance. For j = 8000 A/m 2 overvoltages at the electrodes and the ohmic losses across the electrolyte of η a = 0.063 V, |η c | = 0.050 V and ∆φ e = 0.2069 V are obtained. It can be seen that the ohmic losses in the electrolyte are dominant. The fraction of the total voltage losses accounted for by the ohmic losses increases with increasing electric current density, as the activation overvoltages do not increase linearly with increasing electric current density. This dominance of the ohmic losses is due to the limited ionic conductivity of solid oxide ceramics and the strong dependence of the ohmic losses on the electric current density. The overvoltages at the anode, as well as at the cathode, are in the same order of magnitude, but the overvoltage at the anode is slightly greater. This can be explained by the larger exchange current density of the cathode reaction (j c 0 = 6061.6 A/m 2 ) as compared to the anode reaction (j a 0 = 5711.2 A/m 2 ) at the operating temperature considered here.

Entropy Production Rate
The local entropy production rates for the anode GDL, the cathode GDL, and the electrolyte are shown in Figure 8 along the spatial coordinate at an operating temperature of 1123.15 K. In order to present the curves of the local entropy production rates for different electric current densities in one diagram, the local entropy production rates are each related to their value at the left edge of the respective layer. These absolute values are listed for the electric current densities j = 2000 A/m 2 and j = 8000 A/m 2 in Table 4.  Table 4. Reference values of the local entropy production rate at y = 0 (anode GDL), y = 0.4 · 10 −4 m (electrolyte), and y = 1.9 · 10 −4 m (cathode GDL). The absolute values of the local entropy production rates of all sections increase in amount with increasing electric current density. The individual contribution of the heat flux density to the entropy production in the cathode GDL are negative and assume larger values in magnitude as the electric current density increases. This is due to the large influence of the Peltier effect, whereby heat is transported toward higher temperatures as a result of charge transport in the GDLs, which leads to negative entropy production rates. The overall entropy production is, of course, always positive. As can already be seen in Figure 5, the Peltier effect has a smaller influence at the anode. With increasing electric current density and, thus, increasing temperature gradient, there is a stronger influence of the Fourier effect. The reduced partial pressure of the oxygen compared to the reactants in the anode GDL leads to a similar magnitude. From Figure 8, it can be seen that, in the anode GDL, the local entropy production rate decreases toward the reaction zone. This fact can be attributed to the increasing partial pressure of the of the water, which causes the local entropy production rate to decrease. Due to the higher total pressure existing in the reaction zone, thus, the local entropy production rate as a result of the water transport predominates. In the cathode GDL, the decreasing partial pressure of oxygen toward the reaction layer leads to increasing local entropy production rates. Despite a constant electrical potential in the GDLs, a change in the local entropy productions due to the electrical current density becomes apparent. The decrease of the absolute temperature leads to an increase of the corresponding force and, thus, to an increase of the local entropy productions. The local entropy productions in the GDLs due to the heat flow can be explained by the different influences of the corresponding forces. In the anode GDL, there is an increase in local entropy production at an electric current density of j = 2000 A/m 2 . Here, the increasing corresponding force prevails due to the negatively increasing temperature gradient. In contrast, for j = 8000 A/m 2 , the local entropy production decreases due to the predominant decreasing corresponding force due to the increasing absolute temperature. For the electric current densities of j = 4000 A/m 2 and j = 6000 A/m 2 , larger changes are seen. For these two operating points, very small temperature gradients and, thus, very small absolute local entropy productions are obtained. In the cathode GDL, the decreasing local entropy production is due to the decreasing magnitude of the heat flux density along the spatial coordinate.

Anode GDL
From Table 4, it can be seen that the local entropy production due to the heat flux density is significantly lower than that due to the charge transport. However, with increasing electric current density, a greater deviation from the reference value becomes apparent. At the right edge of the electrolyte, a large part of the heat is transported by the Peltier effect. In the direction of the ion flow, however, the heat to be transported increases due to the irreversibilities, causing the temperature gradient to rise. Since both the flux and its corresponding force increase in the direction of the charge transport, this leads to a strong increase in the local entropy production rate. In the course of the local entropy production rate due to ion transport, a minimum near the left edge of the electrolyte is created. As the electric current density is constant, the cause lies in a variation of the corresponding force. Along the spatial coordinate, the temperature increases monotonically. In contrast, the gradient in the electric potential in the vicinity of the anode side reaction layer indicates smaller potential differences than at the cathode side reaction layer, whereby the increase in absolute temperature initially leads to a decrease in local entropy production, until the corresponding force of the electric potential is predominated by the increase in the amount of the gradient in the electric potential. However, it can be seen from the ordinate that these deviations along the spatial coordinate are quite small. . Simulated curve of the local entropy production rate in the electrolyte due to (a) heat flux density and (b) electrical potential along the spatial coordinate y for different current densities j and temperatures T. In Figure 9, the curves of entropy production rates in the electrolyte for different operating temperatures at low and high electric current density are shown. The reference values are summarized in Table 5. Increasing operating temperature has a positive effect on the irreversibilities in the electrolyte. In particular, the dominant losses due to ion transport decrease significantly with increasing temperature, which can be explained by the strong temperature dependence of the ionic conductivity of the ceramic electrolyte. Figure 10 shows the exergetic efficiencies as a function of operating temperature for a current density of j = 2000 A/m 2 and j = 8000 A/m 2 . A low electric current density, the exergetic efficiency is initially nearly constant. It decreases with increasing temperature, whereas, at high electric current density, the exergetic efficiency increases monotonically with temperature. This can be explained by two effects. As the operating temperature increases, the ionic conductivity increases, while the dominant losses due to ion transport decrease, resulting in an improvement of the exergetic efficiency. In addition, the physical exergy of the supplied material currents increases with increasing temperature, resulting in a decrease in exergetic efficiency. At low electrical current densities, the increase in the supplied exergy predominates, since the losses due to charge transport are low. This causes the exergetic efficiency to decrease as the operating temperature increases, whereas the positive influence of ionic conductivity predominates at high electric current densities and causes the monotonically increasing curve. A direct comparison with values from the literature is difficult due to the different assumptions, parameters and operating conditions. Zouhri et al. [53] investigate the exergetic efficiency of a single SOFC cell by varying the parameters for the anode and operating conditions, and they conclude that the porosity and tortuosity of the anode have no effect. With a temperature of T = 1273 K, a tortuosity of τ = 5, and a porosity of = 0.3, the exergetic efficiency as a function of current density results in ζ(j = 2000 A/m 2 ) = 0.76 and ζ(j = 8000 A/m 2 ) = 0.61. Since the same ionic conductivity was assumed in Zouhri et al.'s study as in this study, one possible reason is the thickness of the electrolyte. At 40 µm, the thickness of the electrolyte is only one third of the case considered here, which leads to lower losses. This fact is also pointed out by Midilli et al. [54], where the exergetic efficiency for an SOFC stack is investigated under variation of the electrolyte parameter. With a temperature of T = 1273 K, an electrolyte thickness of 150 µm, and feeding 97% pure hydrogen, the exergetic efficiency is ζ(j = 2000 A/m 2 ) = 0.48, which, again, is lower than the values gained in this study.

Conclusions
The aim of this study is to theoretically determine the local entropy production rate of an SOFC single cell at an operating temperature of T = 1123.15 K and an operating pressure of p = 1 bar. Local entropy production rates can be identified as loss mechanisms, so that this knowledge is of great importance for cell development and optimization of operating strategies. For this purpose, the mass transport in the GDLs of a 1D SOFC model from Reference [12] is described by improved transport equations. The modeling of the GDLs and the electrolyte are based on the NET approach, whereas the reactivity layers are described with the Butler-Volmer approach. The mass transport starting from pure Stefan-Maxwell diffusion is gradually superimposed with further transport mechanisms (Knudsen diffusion, convection, and thermal diffusion) in order to investigate the influences of the individual transport mechanisms. For a comparison of all approaches, the mass transport is additionally modeled with Fick's law. Furthermore, the transport equations are described as a function of the simulated phenomenological coefficients from Reference [12], with the definitions from Reference [20] as a function of the empirical coefficients. The data for thermal conductivities and ionic conductivity are taken from the literature. The Peltier coefficient is calculated with the approach of Ratkje et al. [14]. For the validation, data for a U, j-characteristic curve of the single cell KeraCell II from the manufacturer Kerafol GmbH and, additionally, self-measured data for the same single cell are used. The calculated curves of temperature, heat flow, electrical potential, and local entropy production rates are discussed. An exergetic analysis is used to assess the efficiency of the cell.
Due to the equimolar composition of the hydrogen and water in the anode, the DGM shows no advantages over the simplified Fick's law in this particular case. The effect of convection is small compared to the other transport mechanisms because there are low gradients in the total pressure in the anode. There is a higher total pressure in the anode reaction layer than in the gas channel, which is why the gas mixture moves in the direction of the gas channel as a result of convection, thus hindering the mass transport of the hydrogen and improving the mass transport of the water. The influence of Stefan-Maxwell diffusion on the cathode gases is larger, as a significantly smaller effective binary diffusion coefficient results for the cathode-side gas mixture. The improved mass transfer of oxygen by convection is very small in its effect since the gas mixture on the cathode side consists predominantly of nitrogen. On the anode and cathode sides, thermal diffusion has almost no influence on mass transport under the operating conditions considered. The differences between the DGM and the NET approach are negligible in this case.
The existing model is able to represent the U, j-characteristics with a minimum error of 0.08% and a maximum error of 0.93% in a range of the electric current density of j = 0-4000 A/m 2 . Due to the negative reaction entropy of the cathode reaction and the irreversibilities, a temperature maximum occurs at the cathode reaction layer. In the cathode GDL, the heat flows toward higher temperatures, which illustrates the significance of the Peltier effects. At the anode reaction layer, net heat is absorbed by the positive reaction entropy. Under an electric current density of j = 4000 A/m 2 , Peltier heat absorption outweighs Fourier heat absorption. The electric potential field is significantly influenced by the ohmic losses in the electrolyte. These are always greater than the overvoltages at the electrodes and cause more than half of the voltage losses. This is also shown in the curves of the local entropy production rates. The entropy production rates as a result of the charge transport in the electrolyte are significantly larger in amount than all other losses. The irreversibilities cause 64.44% of the total losses in the cell in the electrolyte at an electric current density of j = 8000 A/m 2 and an operating temperature of T = 1123. 15.
The irreversibilities at the electrodes, together, cause one third of the losses. As the entropy productions in the GDLs are very low, the primary need is to investigate alternative electrolyte materials.
The exergetic efficiency is fundamentally influenced by the operating temperature and the electrical current density. An increase in exergetic efficiency can be achieved through low current densities and higher operating temperatures.