Next Article in Journal
Governing China’s Clean Energy Transition: Policy Reforms, Flexible Implementation and the Need for Empirical Investigation
Next Article in Special Issue
Electrochemical Mechanism for FeS2/C Composite in Lithium Ion Batteries with Enhanced Reversible Capacity
Previous Article in Journal
Genetic Algorithm-Based Design Optimization of Electromagnetic Valve Actuators in Combustion Engines
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Integrated Planar Solid Oxide Fuel Cell: Steady-State Model of a Bundle and Validation through Single Tube Experimental Data

1
Thermochemical Power Group (TPG)-Department of Civil, Chemical and Environmental Engineering (DICCA), Polytechnic School, University of Genoa, Via Opera Pia 15, Genoa 16145, Italy
2
Rolls-Royce Fuel Cell Systems Limited, SinA-7, PO Box 31, Derby DE24 8BJ, UK
3
Thermochemical Power Group (TPG)–Department of Mechanics, Energetics, Management and Transportation (DIME), Polytechnic School, University of Genoa, Via Montallegro 1, Genoa 16145, Italy
*
Author to whom correspondence should be addressed.
Present address: Carestream Health Italia S.r.l., Palazzina S. Lorenzo, Via al Porto Antico, Genova 16128, Italy
Energies 2015, 8(11), 13231-13254; https://doi.org/10.3390/en81112364
Submission received: 3 July 2015 / Revised: 2 November 2015 / Accepted: 6 November 2015 / Published: 20 November 2015
(This article belongs to the Special Issue Reacting Transport Phenomena in Electrochemical Cells)

Abstract

:
This work focuses on a steady-state model developed for an integrated planar solid oxide fuel cell (IP-SOFC) bundle. In this geometry, several single IP-SOFCs are deposited on a tube and electrically connected in series through interconnections. Then, several tubes are coupled to one another to form a full-sized bundle. A previously-developed and validated electrochemical model is the basis for the development of the tube model, taking into account in detail the presence of active cells, interconnections and dead areas. Mass and energy balance equations are written for the IP-SOFC tube, in the classical form adopted for chemical reactors. Based on the single tube model, a bundle model is developed. Model validation is presented based on single tube current-voltage (I-V) experimental data obtained in a wide range of experimental conditions, i.e., at different temperatures and for different H2/CO/CO2/CH4/H2O/N2 mixtures as the fuel feedstock. The error of the simulation results versus I-V experimental data is less than 1% in most cases, and it grows to a value of 8% only in one case, which is discussed in detail. Finally, we report model predictions of the current density distribution and temperature distribution in a bundle, the latter being a key aspect in view of the mechanical integrity of the IP-SOFC structure.

1. Introduction

Solid oxide fuel cells (SOFCs) are electrochemical reactors that directly convert the energy of the reactants into electrical energy, avoiding the intermediate thermodynamic cycles that characterise the methods traditionally employed for energy conversion. The high efficiency resulting from the direct conversion process, together with interesting features from the point of view of the environmental compatibility (for example, zero NOx emissions), are the main qualities that presently make SOFCs the subject of intensive studies.
SOFCs are now divided into two families, depending on the operating temperature. High temperature SOFCs (HT-SOFCs), operating at 800–1000 °C, are considered for systems in the 100 s of kW to MW power range [1,2,3]. The high operating temperature provides high quality waste heat to recover in cogeneration and bottoming cycles, but places restrictions on the materials that can be used in both the fuel cell and the balance of the plant. To overcome this problem, intermediate temperature SOFCs (IT-SOFCs) have been proposed, operating at 550–800 °C. These are considered for smaller scale applications, where integration with a heat engine is not appropriate. This lower operating temperature is expected to solve various material problems, to extend the range of the acceptable material selection and to improve cell durability [3]. However, lowering the operating temperature also lowers the fuel cell performance, since the kinetics of the electrochemical reaction decreases rapidly as the temperature decreases. Furthermore, the electrode and electrolyte materials become less conductive. A complete review of both HT-SOFC and IT-SOFC features is reported in [1,2,3,4,5,6,7].
The present work focuses on HT-SOFCs. A typical HT-SOFC is formed by a solid electrolyte, which is usually a fluorite-structured oxide ion conductor, i.e., yttria-stabilised zirconia (ZrO2)1-x(Y2O3)x (often referred to as YSZ), coupled to two porous electrodes usually formed by a mixture of YSZ and electronic conducting electrocatalyst particles (usually Ni for the anode and a perovskite-structured oxide lanthanum strontium manganite La1-xSrxMnO3 for the cathode). In order to achieve high voltages, cells are usually connected in series through an electron conducting interconnection, often strontium-doped lanthanum chromite (La0.84Sr0.16CrO3), or chromium-based alloys or ferritic steels. Typically, H2 is the anodic and O2 the cathodic reactant. In each electrode, the electrochemical reaction occurs at the three-phase boundary (TPB) between the electron conductor, ionic conductor and gas:
½ O2 + 2 e → O − −     cathode
H2 + O − − → H2O + 2 e     anode
H2 + ½ O2 → H2O     overall reaction
In the case that the fuel feeding is a partially-reformate mixture of natural gas, then CH4, CO, CO2 and H2O are also present at the anode side, and thus, the water-gas shift (WGS) and methane steam reforming (MSR) reactions can occur:
H2O + CO ↔ CO2 + H2     WGS
CH4 + H2O ↔ CO + 3H2     MSR
Generally speaking, the commonly-adopted shape for fuel cells of any type is a planar geometry. However, since at the high operating temperature of HT-SOFCs the mismatch between the thermal expansion coefficients of the materials forming the layered structure of a planar stack is causing problems in terms of sealing and of mechanical stability (stresses and breakages), other geometries have been developed for HT-SOFCs, in particular the tubular and the integrated planar SOFC (IP-SOFC). The tubular geometry [8], where only one end of the tube is constrained and the other one is free to expand, has the advantage of achieving higher mechanical stability than the planar geometry, but there are two main drawbacks, i.e., that the current paths are rather long, causing significant internal losses, and that the manufacturing process is expensive. The IP-SOFC geometry [9] relies on to the presence of a porous ceramic support (MgAl2O4) onto which the single cells are deposited, as displayed in Figure 1 and Figure 2. As well as in the tubular case, also in this case only one end of the support is constrained, and the other one is free to expand, which ensures good mechanical stability. On the other hand, in IP-SOFCs, current paths are short due to the sequence of striped cells deposited on the support, and fabrication costs are low due to the screen printing deposition technique. As displayed in Figure 2a, several single IP-SOFCs deposited on the support form the so-called multi-cell tube, and several tubes coupled together form the so-called bundle. Further details about the IP-SOFC geometry are given in Section 2 “Experimental Setup”.
Figure 1. Scheme of one single IP-SOFC cell. Charge paths are displayed.
Figure 1. Scheme of one single IP-SOFC cell. Charge paths are displayed.
Energies 08 12364 g001
Figure 2. (a) IP-SOFC scale-up: single cell, planar tube and bundle; (b) IP-SOFC bundle as the basic repeat unit of the experimental strip and block; (c) fuel flow scheme in the IP-SOFC bundle.
Figure 2. (a) IP-SOFC scale-up: single cell, planar tube and bundle; (b) IP-SOFC bundle as the basic repeat unit of the experimental strip and block; (c) fuel flow scheme in the IP-SOFC bundle.
Energies 08 12364 g002
In this paper, we present a mathematical model of the IP-SOFC. After the first seminal SOFC models [10,11], several different approaches have been presented in the literature, which have been reviewed by [12,13,14,15,16,17]. Concerning IP-SOFCs, electrochemical models for the evaluation of voltage and current density distributions have been developed [18,19]. Furthermore, three-dimensional models based on mass, energy and momentum balances within the IP-SOFC structure have been presented [20,21]. However, only in rare cases has an attempt at comparing simulation results and experimental data been made. In this paper, a steady-state model is presented, and its results are discussed in connection with the performance reported experimentally. In particular, this model embeds a previous model developed for the evaluation of the local electrochemical kinetics of the IP-SOFC [18] into a macroscopic model of a full-sized IP-SOFC bundle. In turn, this model is currently included in a plant model [22,23], which is being used to simulate the IP-SOFC hybrid system [24]. The present work focuses on IP-SOFC bundle simulation, and hybrid system simulation results are reported elsewhere [22,23,24].

2. Experimental Setup

The experimental setup has been designed, built and operated by Rolls-Royce Fuel Cell Systems Ltd. The basic element of the overall device is the single cell; a scheme of the geometry of the IP-SOFC concept, along with a description of the current paths through the cell, is reported in Figure 1 (more details can be found in a previous paper [18]). The main geometrical data of the specific experimental IP-SOFC under consideration in this work are reported in Table 1. In particular, each cell measures 9 mm × 60 mm; fifteen stripe-shaped single cells are then printed on each face of a rectangular porous supporting structure (thickness of about 3.5 mm) to give a multi-cell planar tube (Figure 2a). The porous material supporting the active cells allows the fuel to diffuse from the channels up to the anodes, where the electro-oxidation reaction takes place (Figure 1). All of the cells printed over each side of a planar tube are electrically connected in series, by means of conductive interconnections. As a result, on each side, all of the cells are crossed by the same current. The distance between one cell and the subsequent one is equal to the width of the interconnection, which is 6.25 mm. The active cells and the interconnections printed on the tubes do not cover the totality of the available surface, so that inactive edges are left along the contour of each tube. The areas of the tubes where the electrochemical reaction does not occur, i.e., the interconnections and the inactive edges, are covered with a glaze layer in order to prevent a direct chemical reaction between fuel and oxidant.
Table 1. Geometrical details of the IP-SOFC planar tube and bundle.
Table 1. Geometrical details of the IP-SOFC planar tube and bundle.
GeometryDatum
Number of tubes in the bundle (-)10
Number of cells on each face of the tube (-)15
Tube size (×10−3 m)280 × 68
Tube thickness (×10−3 m)3.5
Tube porosity (-)0.3
Active cell size (×10−3 m)9 × 60
Interconnection size (×10−3 m)6.25 × 60
Number of anodic channels15
Width of one anodic channel (×10−3 m)3.47
Thickness of one anodic channel (×10−3 m)1.8
Width of one cathodic channel (×10−3 m)280
Thickness of one cathodic channel (×10−3 m)2
Ten standing tubes are finally placed one in front of the other and connected in series from both an electrical and fuel-feeding point of view: the resulting structure, which is called a bundle (Figure 2a), represents the complete experimental set taken into account in this work. The bundle is then the basic repeat unit of the strip and of the block (Figure 2b), which is, in turn, the main component of the plant. In this work, attention is focused on the single tube and bundle, and Figure 2c displays the fuel flow scheme: fuel is fed into the first tube of the bundle through a number of parallel channels (anodic channels) bored inside the porous supporting tube, and then, it flows in series through all ten tubes, passing from one to the subsequent one through a ceramic connection (visible in Figure 2b). The fuel path along the bundle is up-and-down (serpentine) through the ten tubes that form it. Going more into detail, the fuel is fed into the first tube from the bottom to the top, and then, thanks to the ceramic connection, it flows into the second tube in the opposite direction. This scheme is repeated also for the other tubes, so that all of the odd tubes are crossed by the fuel upwards, while the even ones are crossed downwards. The electrical current has a similar path through the tubes that form the bundle, due to the presence of additional electrical connections positioned at the extremities of the tubes. These electrical connections are manufactured through highly conductive stripes (visible in Figure 2b), adhering on the first and last interconnection on both faces of each single tube. Then, for each tube, the two stripes at each end are welded to the corresponding stripes of the subsequent tube, creating an electrical link. This allows a serpentine path for the electrical current along the sequence of tubes forming the bundle, and thus, all of the tubes of the bundle are in series from the electrical point of view.
As far as the oxidant is concerned, air is swept in the spaces between the tubes (which are referred to as cathodic channels), in the direction perpendicular to the direction of the fuel flow. In this way, the width of the cathodic channels is equal to the length of the tube, i.e., 280 mm, and the thickness is equal to the distance between two adjacent tubes, i.e., 2 mm. To sum up, electricity and fuel cross the ten tubes in series, whereas the air flows in parallel among them.
The subsequent Section 4.1 “Model Validation” reports experimental data collected from single tubes housed in a muffle furnace kept at fixed temperature and adapted to enable fuel and air to be fed to the test tube, as well as to allow electrical measurements to be made. Concerning the operating conditions, here, operation at atmospheric pressure is considered. The fuel is humidified until saturation at ambient conditions (water molar fraction ξH2O = 0.03) before being preheated and fed into the experimental setup. The air is also preheated before being fed into the cathode compartment.

3. Mathematical Model

3.1. Model Development

The model is based on the following hypotheses:
  • For each tube, a steady-state 2D model has been developed along the x and y coordinates (Figure 3).
  • Heat exchange by radiation (i) between gases and solid in each single tube and (ii) between facing tubes inside the bundle has been neglected. Heat exchange of any type between tubes or bundles and the surrounding environment is neglected.
  • Gravitational effects are not taken into account.
  • A uniform distribution of the anodic flow rate among the inner channels of the tubes is assumed.
  • The air flows only in the y direction with no transversal flows; moreover, the total inlet flow rate is equally subdivided among the bundle cathodic channels.
  • All of the single cells supply the same electrical current. Both voltage and current are zero over the inactive edges of the tubes.
  • Voltage can be considered uniform over each interconnection and active cell, due to the presence of elements with high electrical conductivity.
  • Mass and energy balances of the gas streams are written in plug-flow form along their flow directions, because, even if the flow regime is laminar, the Péclet number is at least 20.
  • MSR and WGS reactions are assumed at equilibrium in the anode channels.
  • No carbon deposition is considered in the anode.
  • The electrochemical oxidation of H2 is the only electrochemical reaction occurring at the anode.
  • The bundle is simulated as a sequence of a number of tubes, placed one after the other.
Figure 3. Discretization mesh for model integration.
Figure 3. Discretization mesh for model integration.
Energies 08 12364 g003
Concerning Hypothesis 1, Figure 3 shows that each tube has been divided into mesh elements along the fuel flow direction (x direction) and the air flow direction (y direction), while all of the quantities are assumed to be constant along its thickness. The assumption is justified by the fact that the tubes are extremely thin (3.5 mm). As a consequence of this hypothesis, the two faces of each tube are considered to have identical temperature. Concerning Hypothesis 2, heat exchange of any type between tubes or bundles and the surrounding environment is neglected, including any possible heat conduction toward surrounding structures, due, for example, to the presence of pipes or manifolds for the fuel feeding, devices to support or fix the tubes and wires for the electrical connections. This ideal adiabatic operating condition is approached when the bundle is operated in the plant, surrounded by a thick insulation layer. According to Hypothesis 9, MSR and WGS reactions are considered at equilibrium within the fuel flowing into the anode channels [25]. A possible undesired chemical reaction, which is neglected in the present work as stated in Hypothesis 10, is carbon deposition, which can take place when carbonaceous components are present in the anode gas. A complete treatment of carbon deposition in SOFCs is given by [26]. Concerning Hypothesis 11, also CO and CH4, in addition to H2, can be electrochemically oxidised in SOFCs; however, their electrochemical reaction rates are considerably slower than that of H2 [1,2], and therefore, these reactions are neglected in the present work. Finally, according to Hypotheses 1 and 12, the model equations presented below form the mathematical model of each single tube, and then, the bundle is simulated as a number of tubes, placed one after the other.

3.1.1. Local Electrochemical Kinetics

As far as the simulation of the local electrochemical kinetics is concerned, its scope is the evaluation of the voltage distribution all over the bundle when a fixed electrical current is flowing through it. To this end, an electrochemical model has been developed, which is based on the evaluation of the local thermodynamic reversible potential (Nernst potential) [12]:
V N e r n s t = Δ G n e F = Δ G 0 n e F R g T n e F ln [ i ( p i 0 ) ν i ]
where pi° are the partial pressures of the reactants and products in the bulk anodic and cathodic flows; Rg is the universal gas constant; T the operating temperature and F Faraday’s constant. For SOFCs, the open circuit voltage (OCV, i.e., the voltage difference between the electrodes when no current is flowing through the fuel cell) is equal to the Nernst potential VNernst [12], indicating that reversibility is obtained at the open circuit (which is not the case for other types of fuel cells, namely PEMFCs, i.e., proton exchange membrane fuel cells). Then, when current flows, irreversibilities arise and the operating voltage decreases [12]:
V = V N e r n s t η
where the overall overpotential η includes ohmic, activation and concentration losses [12]:
η = η o h m + η a c t + η c o n c
A previous paper [18] reports a detailed mathematical treatment of the different kinds of losses in IP-SOFCs: (i) ohmic loss in the anode, electrolyte and cathode; (ii) activation loss in the anode; (iii) activation loss in the cathode; (iv) concentration loss in the anode; and (v) concentration loss in the cathode. Concerning the ohmic term, the model includes an evaluation of the ohmic resistance of the interconnections that link one cell to another in a single tube [18]. A contact resistance term, accounting for defects in the adhesion among different layers of the fuel cells [27], which is often treated as an adjustable parameter, is set at zero in this model. Further, the resistances of the highly-conductive stripes forming the electrical connections between the subsequent tubes of the bundle are neglected. Concerning activation losses, they are expressed through the Butler–Volmer equation, which embeds, as an electrochemical parameter, the exchange current density io. The complete mathematical formulation can be found in [18], and here, we wish to recall only the equation adopted for the anodic exchange current density, which is expressed as a function of the Arrhenius law and of the composition of the reacting gases [18,28]:
i 0 , a = γ a ( p H 2 p r e f ) ( p H 2 O p r e f ) exp ( E a c t , a R g T )
where γa is the anodic pre-exponential factor; pH2, pH2O, pref are hydrogen, water and reference partial pressures and Eact,a is the anodic activation energy. Several different expressions have been proposed in the literature for io,a [28]; in particular, with some correlations, io,a decreases when decreasing pH2O, while others display the opposite behaviour. Here, an expression with io,a directly proportional to pH2O is chosen (Equation (4)), where the anodic exchange current density is small when either the H2 partial pressure or the H2O partial pressure is low. When the anodic exchange current density is low, the anodic activation loss is high, and also, it shows a deviation from linearity. This deviation from linearity of activation losses is analogous to that frequently observed at operating temperatures below 900 °C and has been highlighted and discussed on theoretical grounds [29]. Marked bending of I-V curves caused by non-linear activation losses has often been observed experimentally [30].
Equations (1)–(3), together with the equations developed in [18], which express ohmic, activation (including Equation (4)), and concentration losses as a function of the physico-chemical parameters, form the local electrochemical kinetics of the IP-SOFC electrochemical reactor and are all embedded into the present model.

3.1.2. Mass Balances

Since the local rate of the electrochemical reaction is proportional to the electrical current, the model calculates the mass balances as a function of the current density, according to Faraday’s law. At the cathodic side, oxygen is the only reacting species, and the mass balance reads as follows [12]:
d W O 2 d y = J 2 n e F b c
where WO2 denotes the molar flow of oxygen; J the local value of current density (being zero in those areas of the tube where no current is supplied) and bc the thickness of the cathodic channel.
As for the anodic side, since the fuel flows inside a number of channels (nch), it is possible to refer the mass balances to one single channel, assuming that the consumption of fuel in that channel is related to the electrochemical reaction occurring on a vertical stripe of width L, equal to 1/nch of the overall width of the active cells. With this assumption, the anodic mass balances read as follows [12]:
d W i d x = r i
where the subscript i denotes the reacting species; while r represents the reaction rate. The expressions of the reacting rates for the anodic species are reported below:
r H 2 O = J n e F b a L a a r M S R r W G S
r H 2 = J n e F b a L a a + 3 r M S R + r W G S
r C O = r M S R r W G S
r C O 2 = r W G S
r C H 4 = r M S R
where aa and ba are the anode channel width and thickness. Equations (7)–(11) relate the reaction rates of H2O, H2, CO, CO2 and CH4 to the MSR and WGS reaction rates rMSR and rWGS. The MSR and WGS reactions are considered at equilibrium (Hypothesis 9), so that rMSR and rWGS are calculated as a function of the MSR and WGS equilibrium constants and as a function of temperature and composition, which changes continuously due to the subtraction of hydrogen and the production of water related to the electrochemical reaction.

3.1.3. Energy Balances of the Gases

The temperatures of the air and of the fuel change continuously during their transit inside their respective channels, mainly due to the convective heat exchanges with the solid structure. Since all of the chemical and electrochemical reactions are catalysed by the electrode materials, the related heat dissipation term is included in the energy balance of the solid. At the cathodic side, the energy balance for the air is written as follows [12]:
k W k c p , k T c y + c p , O 2 W O 2 y ( T c T s ) + h c b c ( T c T s ) = 0
In Equation (12), the first term gives the variation of temperature of the cathodic gas (Tc) along the y direction; the second term takes into consideration the mass transfer of oxygen from the cathodic gas to the solid due to the electrochemical reaction (calculated through the related mass balance, Equation (5)). The third term accounts for convective heat exchange between gas and solid, hc being the convective heat transfer coefficient. hc is calculated as a function of the Nusselt number (Nu), the thermal conductivity of the gas (λ) and a characteristic length (dh):
h = N u λ d h
In general, the Nusselt number is a function of the channel geometry, the Reynolds number (Re) and the Prandtl number (Pr). In this case, as established laminar flow holds for both the anodic and cathodic gases, the Nusselt number is assumed to be a constant; dh is computed as the hydraulic diameter of the channel.
As far as the fuel is concerned, it is possible to refer its energy balance to one single channel, as previously done for the anodic mass balances (Equation (6)). The anodic energy balance reads [12]:
i W i c p , i T a x + i c p , i W i x ( T a T s ) + h a B b a ( T a T s ) = 0
where Ta is the temperature of the anodic gas. The second term, related to mass transfer between solid and gas, here includes the contributions of electrochemical, MSR and WGS reactions. Again, this term is calculated through the related mass balance (Equation (6)). The coefficient B appearing in the third (convective) term is the ratio between the gas-solid heat exchange area and the cell area, and it is zero in those regions of the tubes where there are no channels (edges). ha is calculated through the same procedure applied for the calculation of hc.

3.1.4. Energy Balance of the Solid

The energy balance of the solid structure reads as follows [12]:
1 s h c ( T c T s ) + 1 s h a B ( T a T s ) 1 s ( j Δ H j n e F + V ) J + K x 2 T s x 2 + K y 2 T s y 2 = 0
This balance takes into account the convective heat exchanges with the cathodic and anodic gases (first and second terms, respectively), the thermal contribution due to electrochemical, MSR and WGS reactions (third term) and the heat conduction within the solid structure (fourth and fifth terms). Further details about the calculation of the thermal conductivities Kx and Ky for the layered IP-SOFC structure are reported in Section 3.3 “Modelling Heat Conduction through the Solid Structure”.

3.2. Geometry Discretization

A 2D discretisation of the tubes has been operated (Figure 3): each tube of the bundle has been subdivided into a number of elements both in the x and y direction. As far as the fuel flow direction (x direction) is concerned, the integration step has been chosen in such a way to consider, at each step, an active cell or an interconnection. As for the edges at the top and at the bottom of the tubes, they have been meshed, as well, with an integration step equal to the average between the previous two steps. Thus, the integration step along the fuel flow direction is not fixed, but in general, three different steps (edges, active cells and interconnections) alternate, depending on the position on the tube. As regards the air flow direction (y direction), the active area of each tube, where cells and interconnections are present, has been meshed with a number of steps equal to the number of anodic channels. For the side borders, an integer number of steps, with a size as close as possible to that used in the active region, has been used.
Following the discretization of the geometry, the equations are discretised, as well, through a finite difference approach, coupled to a relaxation method for the energy balance of the solid structure.

3.3. Modelling Heat Conduction through the Solid Structure

According to Hypothesis 1, the energy balance of the solid (Equation (15)) includes only the in-plane heat conduction along the x and y directions. When discretizing Equation (15) through a finite difference approach, for each integration element, the conductive thermal fluxes must be calculated toward (or from) the two adjacent elements along the x direction and the two ones along the y direction. Since the elements on the perimeter of the tubes are supposed not to exchange heat by conduction toward the surroundings, in our simulation, these elements exchange heat only with three or two adjacent elements, depending on their position. In the model, it has been assumed that heat flows from the centre of one element to the centre of the adjacent element, through a structure composed of a superimposition of layers characterised by different materials and different thickness (supporting tubes, electrodes, electrolyte and current collectors). Thus, the thermal flow has to cross different resistances depending on the direction of propagation; for this reason, it is necessary to evaluate, in general, four different effective thermal resistivities for each integration element (one for each pair of adjacent elements, except for the elements on the perimeter of a tube). Each integration element, except for those placed on the inactive edges, is characterised by a different arrangement of the layers depending on whether it is a part of an active cell or of an interconnect. The elements over the inactive edges are ceramic porous supporting tubes (covered by a thin layer of glaze). As a consequence, there are three different kinds of elements that can come into contact with each other in a number of different ways. For each possible arrangement, the in-plane thermal resistivity Rth,12 from the centre of the first element to the centre of the second one is calculated. If the effective in-plane thermal conductivities of the two elements of a pair are supposed to be known (K1,eff and K2,eff), the thermal resistivity Rth,12 is calculated as follows:
R t h , 12 = d 1 k 1 , e f f 1 + d 2 k 2 , e f f 1 d 1 + d 2
where d1 and d2 are the distances covered by the thermal flux within Elements 1 and 2, respectively. In turn, the effective in-plane thermal conductivity Keff of an element formed by N layers, each one with thermal conductivity Ki and thickness ti, is evaluated through a weighted mean:
k e f f = i = 1 N t i k i i = 1 N t i
For the elements over the inactive edges, where only the material of the supporting tube is supposed to be present, the effective thermal conductivity is calculated as:
k s t , e f f = k s t ( 1 ε s t )
where εst is the porosity of the supporting tube.
The procedure described above allows one to evaluate the resistivity Rth,12 against conductive heat transfer through two adjacent integration elements. If the two elements are located along the x direction, then:
K x , 12 = R t h , 12 1
Instead, if the two elements are located along the y direction, then:
K y , 12 = R t h , 12 1
and the resulting thermal conductivities Kx,12 and Ky,12 are then embedded into the discretised form of Equation (15).

3.4. Model Integration

The method used for the discretization of the equations is based on finite differences for the mass and energy balances (Equations (5)–(20)), coupled to a relaxation method additionally applied to the energy balance of the solid structure (Equation (15)). The relaxation method consists of simulating the evolution of the solid temperature over time, starting from a uniform distribution given as input data, up to a steady-state condition. The temperature variation between two iterations is calculated, over each element, assuming a time step dt equal to 0.2 s and a specific heat for the solid structure equal to 700 J/kg·°C. The relaxation method is the outer iterative loop implemented in the tube simulation code, with an inner iterative loop nested inside. The inner loop is developed based on the hypotheses that both the voltage and current are zero over the inactive edges of the tubes (Hypothesis 6), and that voltage can be considered uniform over each interconnection and active cell (Hypothesis 7). Since the electrical conductivity of the interconnection is practically independent of temperature, also the current distribution is uniform over the interconnections. On the contrary, the electrical current generated by the electrochemical reaction on each cell is not uniform since: (i) the Nernst voltage is non-uniform, due to the non-uniformity of the temperature and reactant concentration; and (ii) voltage losses (ohmic, activation and concentration) are all dependent on temperature; activation and concentration losses also depend on the reactant composition. All simulations are run by imposing an overall operating current, which is the same for all of the cells, while the voltage is a result of the calculations and varies from one cell to another. While the overall current supplied by each cell is fixed, the current distribution within each single cell is evaluated within the inner loop. To this end, the inner loop implements an iterative procedure that varies the operating voltage of each cell (uniform over the cell itself), until its value is consistent with the overall current generated by the cell itself (fixed).
Finally, according to Hypothesis 12, the bundle is simulated as a sequence of a number of planar tubes, placed one after another, with the fuel composition at the outlet of the previous tube being the inlet composition of the subsequent one. Clearly, the numerical tool can simulate either a single tube or a full-sized bundle.
Both isothermal operation and adiabatic thermal boundary conditions can be accounted for. In the previous case, the temperatures of the solid and gases are assumed to be uniform and constant over time, and the energy balances (Equations (12)–(20)) are not solved through the outer iteration loop. From the experimental point of view, this corresponds to the situation where the single tube or full bundle is operated in a furnace kept at a fixed temperature. If adiabatic simulation is selected, the energy balances (Equations (12)–(20)) are actually solved by the simulation code. The boundary conditions for the energy balances are coherent with Hypothesis 2, i.e., heat exchange of any type between tubes or bundles and the surrounding environment is neglected.
The model equations are integrated through an in-house developed code, written in Fortran 90, which features high execution speed (tens of run-time seconds for a model containing about 7000 grid elements running on an Intel CORE i7 CPU running at 2.8 GHz with 12 GB of RAM). High execution speed is a key feature of this model in view of the integration into the whole plant model [24].

3.5. Model Inputs and Outputs

The inputs of the model are physical parameters, geometrical data and operating conditions. The main geometrical data of the simulated bundle are reported in Table 1 and include the dimensions of single cells and tubes of the simulated bundle. As far as the physical parameters of the model are concerned, they have been chosen on the basis of previously-published work [18,28] or literature values [31]. In particular, the data regarding the electrochemical kinetics and the electrical conductivities of the material are not reported here, as they are the same as used in a previous work on the electrochemical kinetics of IP-SOFCs [18], where they were also validated against experimental data. Concerning the parameters regarding the thermal part of the model, the Nusselt number is assumed constant and equal to 4, and the thermal conductivities of the solid components of the IP-SOFC are taken from the literature [31]: being all ceramic materials, they range between 2 and 3 W/m·°C, and they do not vary significantly in the range of temperatures under consideration. Finally, the thermal conductivities of the gases are modelled through literature temperature-dependent equations [32]. In this way, the model does not contain any adjustable parameter. Table 2 reports the base case operating conditions used in the simulations and includes inlet temperature and flow rate for the fuel (prior to humidification). Various fuels with different compositions are employed, which are duly specified case by case in the following Section 4 “Results and Discussion”. Table 2 reports also the oxidant (pure air) overall flow rate, as well as the operating pressure.
Table 2. Base case operating conditions.
Table 2. Base case operating conditions.
Operating parameterDatum
Fuel flowrate (NL/min)6
Oxidant (air) flowrate (NL/min)500
Fuel inlet temperature (°C)930
Air inlet temperature (°C)910
Pressure (atm)1
The results of the model are distributions of current density and gas and solid temperatures, actual voltage and Nernst voltage over the tubes, as well as all of the voltage losses (ohmic, activation and concentration).

4. Results and Discussion

Two types of results are reported: (i) isothermal single tube simulations, compared to experimental data (Section 4.1 “Model Validation”); (ii) adiabatic bundle simulations (Section 4.2 “Model Predictions”). In any case, when not explicitly stated, base case operating conditions apply (Table 2).

4.1. Model Validation

An extensive validation of the model has been carried out comparing the simulation results of current-voltage (I-V) curves to the experimental data provided by Rolls-Royce Fuel Cell System Ltd. In order to carry out this validation, several different anodic mixtures have been used; in addition, also the effects of a variation in both the anodic flow rate and operating temperature have been investigated. Since the experimental data have been collected from single tubes operated in an oven kept at a fixed temperature, the simulations have been performed in the isothermal modality.

4.1.1. Validation with Binary Mixtures at the Anodic Side

The first validation of the model has been carried out using binary mixtures at the anodic side. In the following, for simplicity, all of the compositions reported are the compositions of the dry anodic mixtures. However, since in experimental practice, all of the anodic mixtures are humidified, humidification is considered in the simulations (ξH2O = 0.03 in the feeding fuel). In Figure 4, a comparison is reported between experimental and simulated results in the case of a fuel composed of H2 and N2, varying the composition, the temperature and the flow rate. Figure 4a shows the I-V curves obtained with a fuel flow rate equal to 4 NL/min and a tube temperature maintained ideally constant at 950 °C in an oven. The flowrate of air fed at the cathodic side is considered to be 200 NL/min in the simulations; in the experimental tests, it is not measured with accuracy, because, due to the large excess compared to the stoichiometric values, the oxygen molar fraction is practically uniform and equal to 0.21 all over the tube. The agreement obtained between simulation and experimental results is very satisfactory (overall average error <0.8%), especially considering that the model does not contain any adjustable parameter. As expected, both experimental and simulated I-V curves improve with increasing H2 content in the fuel. An increase of H2 concentration in the fuel is expected to increase the OCV and to decrease both anodic activation and anode concentration losses, with all of the other loss terms (ohmic, cathodic activation and concentration) remaining unchanged. The fact that Figure 4a shows a good agreement between experimental and simulation results for all of the compositions under consideration confirms in particular the reliability of the simulation of the previous three terms (i.e., OCV, anodic activation and concentration). Starting from the general consideration that Figure 4a demonstrates a good agreement, some comments can be made regarding some specific aspects. In particular, the OCVs are a little higher in the simulations compared to the experimental values; we believe that this discrepancy is due to the fact that the tubes may have some defects or some porosities through which a leakage of gas between the cathodic and anodic side may occur, which is not taken into account in the model. In particular, due to the properties of hydrogen and, in particular, its high diffusivity, we believe that the main effect of leakage is a loss of hydrogen from the anodic side, with a decrease of the hydrogen molar fraction and the consequent increase of the water molar fraction at the anodic side, which tends to lower the experimental value of the OCV (Equation (1)). Variations in the fuel composition, even if small, have a significant effect on the OCV, since the argument of the logarithmic term in the Nernst law (Equation (1)) is close to zero, because of the low partial pressure of water. The hypothesis that hydrogen leakage is responsible for the discrepancy between simulated and experimental OCVs reported in Figure 4a is supported by the fact that the discrepancy is higher at high H2 contents of the fuel (average error about 2% for ξH2 = 0.8 in the fuel), while it almost disappears when the molar fraction of H2 in the fuel is low (average error about 0.1% for ξH2 = 0.2 in the fuel), since leakage is enhanced by high H2 molar fractions in the fuel. Another remark regards the curvature of the I-V curves. Figure 4a shows that the curvature becomes more visible when decreasing ξH2 in the fuel, and this can be explained in terms of increased anodic activation loss. Indeed, an additional analysis (not reported here) of the relative importance of all of the loss terms in the operating conditions of Figure 4a shows that the anodic activation loss exhibits a marked increase when decreasing the hydrogen content, with concentration losses remaining lower. Indeed, in the operating region of Figure 4a, the voltage is above 0.4 V, which means that we are far away from the limiting current, and thus, it is not surprising that concentration losses are lower than the other sources of loss. The increase of anodic activation mentioned above, occurring when decreasing the H2 molar fraction in the fuel, is also associated with a more marked non-linear behaviour as a function of current density, as already observed in Section 3.1.1 “Local Electrochemical Kinetics”. Thus, in the case displayed in Figure 4a, we explain the increased deviation from the linearity of I-V curves when decreasing the hydrogen content, in terms of increased anodic activation. For IP-SOFCs, we have already demonstrated that our electrochemical kinetics subroutine captures this effect often displayed by the experimental data [18], and this is further demonstrated by Figure 4a. A further remark is that in the high current density range (but still far away from the regime dominated by concentration limitation), Figure 4a displays some deviations between experimental data and simulations, and we believe this is due to the Joule effect, which, at high current densities, brings the temperature of the tubes a little bit over the temperature of the furnace, with a consequent decrease in the internal resistances and an improvement of the performance. This thermal effect is not simulated by the model, since the simulation is assumed to be perfectly isothermal.
Figure 4b shows a comparison between the simulation and experiments at 890 °C, again using a fuel flow rate equal to 4 NL/min. As is reasonable to expect, the fuel cell performance decreases at a lower operating temperature, due to globally higher internal losses [18]. In this case, the agreement between calculated and measured OCVs is better, since, being that the temperature is lower, the gaseous diffusion between the anodic and cathodic side has a less significant effect. The curvature of the I-V curves, both experimental and simulated, is more pronounced in this case compared to 950 °C, because of the anodic activation losses, which increase further and further deviate from linearity as a function of the current density [18,29] because of the lower temperature. Furthermore, the upward deviation of the experimental data at high current densities is more evident here compared to 950 °C, since the Joule effect, responsible for such a deviation, is higher because of the higher internal resistances due to the lower temperature.
Figure 4. Validation of I-V curves obtained with H2/N2 fuel mixtures. Symbols: experimental. Lines: simulations. Captions indicate molar fractions of H2 and N2 in the dry fuel for cases (a) operating temperature 950 °C, fuel flow rate 4 NL/min (b) operating temperature 890°C, fuel flow rate 4 NL/min and (c) operating temperature 890°C, fuel flow rate 2 NL/min.
Figure 4. Validation of I-V curves obtained with H2/N2 fuel mixtures. Symbols: experimental. Lines: simulations. Captions indicate molar fractions of H2 and N2 in the dry fuel for cases (a) operating temperature 950 °C, fuel flow rate 4 NL/min (b) operating temperature 890°C, fuel flow rate 4 NL/min and (c) operating temperature 890°C, fuel flow rate 2 NL/min.
Energies 08 12364 g004
In Figure 4c, a validation of the model is shown in the case in which the temperature is 890 °C and the fuel flow rate is 2 NL/min. Both the simulation and the experimental data show very small differences compared to the case in Figure 4b, and this is due to the superabundant fuel flow rates.

4.1.2. Further Validations

Further validations of the model have been carried out using humidified H2/CO/CO2 fuel mixtures, pre-heated at 950 °C, and keeping the fuel flow rate at 10 NL/min. The fuel compositions before and after humidification and WGS reaction are reported in Table 3, and the corresponding fuel cell performance is reported in Figure 5.
Table 3. Dry H2/CO/CO2 feeding fuel compositions (experimental), fuel compositions after ambient temperature humidification and after WGS reaction at 950 °C (calculated). Related to the results in Figure 5.
Table 3. Dry H2/CO/CO2 feeding fuel compositions (experimental), fuel compositions after ambient temperature humidification and after WGS reaction at 950 °C (calculated). Related to the results in Figure 5.
CaseFeeding fuel composition (dry)Feeding fuel composition (after humidification)Fuel composition (after humidification and WGS reaction)
ξH2ξH2OξCOξCO2ξH2ξH2OξCOξCO2ξH2ξH2OξCOξCO2
110000.970.03000.970.0300
20.66700.2530.080.6470.030.2450.0780.5990.0780.2920.031
30.66700.1330.20.6470.030.1290.1940.5170.160.2580.065
Figure 5. Further model validation with three different H2/CO/CO2 fuel compositions (Table 3). Operating temperature: 950 °C. Fuel flow rate 10 Nl/min. Symbols: experimental. Lines: simulations.
Figure 5. Further model validation with three different H2/CO/CO2 fuel compositions (Table 3). Operating temperature: 950 °C. Fuel flow rate 10 Nl/min. Symbols: experimental. Lines: simulations.
Energies 08 12364 g005
From a qualitative point of view, both the experimental and simulated I-V curves show an increase of performance by increasing the hydrogen content of the fuel. The simulated I-V curves agree quite well with the experimental ones, with an overall average error <0.7%. Going more into detail, we observe that in Case 1 of Figure 5, humidified hydrogen is used as the inlet fuel, and a behaviour similar to that reported in Figure 4a is reported, with a disagreement between experimental and simulated performance close to the OCV. The experimental OCV is significantly lower than the OCV calculated by the model, and the discrepancy is about 2%, similar to the values previously reported for the OCVs in Figure 4a. As well as for the cases reported in Figure 4a, also in this case we believe that the cause of this discrepancy is leakage, which causes a decrease of ξH2 and an increase of ξH2O at the anodic side. For an interpretation of the other results reported in Figure 5, we must take into account the WGS reaction, occurring within the anodic compartment when CO and CO2 are present in the fuel feeding mixture (Cases 2 and 3). Table 3 shows the fuel composition, calculated from the humidified feeding fuel composition considering the WGS reaction at equilibrium. Since we assume that the WGS reaction occurs very quickly at the tube entrance, this is in practice the fuel composition at the tube inlet. The data reported in Table 3 show that in Cases 2 and 3, also due to the WGS reaction, ξH2 is lower and ξH2O higher than in Case 1. In other words, the WGS reaction proceeds backwards (referring to the stoichiometry in Reaction (IV)), and consumes CO2 and H2 to produce H2O and CO. As a consequence, leakage is expected to impact the experimental results less for two reasons: (i) the hydrogen content at the anodic side is lower, and thus, leakage is expected to be lower, as well; and (ii) since the WGS reaction makes the H2O partial pressure increase at the anodic side, this results in an increase of the argument of the logarithmic term in the Nernst law (Equation (1)), so that the OCV becomes less sensitive to the small variations in the anodic gas composition due to the leakage. Indeed, Figure 5 shows that when CO and CO2 are present in the feeding fuel mixture (Cases 2 and 3), the discrepancy at OCV is smaller (less than 0.4% in both cases), and the agreement between simulated and experimental data is good, both at OCV and also under load.
Since the simulation results reported in Figure 5 for Cases 2 and 3 have been obtained considering the WGS reaction at equilibrium, then the good agreement demonstrated in Figure 5 also proves that the hypothesis of equilibrium for the WGS reaction is correct.
A further comment regards the deviation from linearity, which is visible in Case 1 in both the experimental and simulated I-V curves. Concerning the simulation, this is explained considering that the product ξH2·ξH2O in the feeding fuel after humidification and WGS reaction has values of 2.91 × 10−2, 4.67 × 10−2 and 8.27 × 10−2 in Cases 1, 2 and 3, respectively (calculated from the data reported in Table 3). Since the product ξH2·ξH2O is proportional to the product pH2·pH2O appearing in Equation (4), this means that the low water partial pressure in Case 1 causes a lower anodic exchange current density than in Cases 2 and 3 and, thus, a higher anodic activation loss, which is associated with more marked bending. Agreement between experimental data and simulation results justifies the choice of Equation (4) for the expression of activation losses for the IP-SOFC cell and, in particular, the choice of a correlation for the anodic exchange current density io,a, which decreases when decreasing pH2O. In other types of fuel cells taken into consideration in previous works [28], experimental data displayed better agreement with simulation results when the anodic exchange current density io,a was expressed through a correlation, which increased with decreasing pH2O. This different behaviour depends on the different choice of materials and manufacturing procedures in the different types of experimental SOFCs.
Simulated and experimental results obtained with fuel mixtures containing CH4 are reported in Figure 6. The first remark is that in all cases where CH4 is present in the feeding fuel, carbon deposition is expected to take place at the anodic side over time [26]; the experimental and simulation results reported here display only the initial performance. Table 4 reports the fuel compositions (i) prior to humidification and (ii) after humidification, considering both MSR and WGS reactions at thermodynamic equilibrium at 950 °C, which is, in practice, the fuel composition evaluated by the model immediately after the tube entrance. Under OCV, this composition remains unchanged all over the tube, while, under operation, water is generated at the anodic side by the electrochemical reaction, which leads to further progress of the MSR reaction and further production of H2 along the fuel direction.
Figure 6. Further model validation using three different fuel compositions, containing methane (Table 4). Operating temperature: 950 °C. Fuel flow rate 10 Nl/min. Symbols: experimental. Lines: simulations.
Figure 6. Further model validation using three different fuel compositions, containing methane (Table 4). Operating temperature: 950 °C. Fuel flow rate 10 Nl/min. Symbols: experimental. Lines: simulations.
Energies 08 12364 g006
Table 4. Dry CH4/H2/CO/CO2 feeding fuel compositions (experimental) and fuel composition after ambient temperature humidification and methane steam reforming (MSR) and WGS reactions at 950 °C (calculated). Related to the results in Figure 6.
Table 4. Dry CH4/H2/CO/CO2 feeding fuel compositions (experimental) and fuel composition after ambient temperature humidification and methane steam reforming (MSR) and WGS reactions at 950 °C (calculated). Related to the results in Figure 6.
CaseFeeding fuel composition (dry)Fuel composition (after humidification, WGS and MSR reactions)
ξCH4ξH2ξH2OξCOξCO2ξCH4ξH2ξH2OξCOξCO2
10.190.5400.2700.1470.5781.5 × 10−60.2755.2 × 10−7
20.0140.65600.3300.0030.660.0130.3190.005
300.66700.1870.1461.5 × 10−60.5520.1250.2770.046
The first remark is about Case 1 of Figure 6, where, due to the very efficient MSR reaction, the calculated data reported in Table 4 show that practically all of the humidification water is consumed at the tube entrance. From the experimental side, in Case 1 of Figure 6, the data are missing at open circuit and at low current densities, due to the OCV being extremely unstable and difficult to measure. The experimental observations are in line with the theoretical findings, since, in the absence of water, the argument of the logarithmic term in the Nernst law (Equation (1)) is close to zero, and small stochastic fluctuations of the water content in the fuel lead to large variations of the Nernst potential. An anomaly regards Case 1 in Figure 6, where the simulated I-V is lower than the experimental one, so that, beyond a certain value of current density, the I-V curve crosses the simulated I-V obtained in Case 2. This intersection between the two curves does not occur in the experimental data, since the experimental data show that the I-V curve regarding Case 1 is above the I-V curve regarding Case 2 at any value of current density. Since this intersection occurs only in the simulations and only when methane is present in the feeding fuel feedstock, we propose the explanation that the hypothesis of the equilibrium of the MSR reaction is not fully accurate. Going more into detail, marked bending is visible only in the simulated I-V curve obtained for Case 1, which is the case where the water content of the fuel is the smallest, and as a consequence, in the simulations, the anodic exchange current density io,a is the smallest and the anodic activation loss the highest. In this case, the hypothesis of equilibrium of the MSR reaction may lead to an overestimation of the degree of water consumption in the fuel feeding mixture, and this can explain the excessively high anodic activation losses calculated by the simulation model and, finally, the underestimation of the electrochemical performance, with the simulated I-V being lower than the experimental one. The hypothesis of equilibrium in SOFCs has received some criticisms in some studies on this subject [33,34], even if the deviation from thermodynamic equilibrium is expected to be greater with IT-SOFCs rather than with HT-SOFCs.
As a final remark, from a quantitative point of view the differences between experimental and simulations results for Case 1 are quite small (overall average error <0.6%). Concerning the intersection between the I-V curves obtained in Cases 1 and 2, from a qualitative point of view both the simulation and the experimental data show the same behaviour, i.e., the tendency of the two I-V curves to overlap at high current densities.

4.2. Model Predictions

The model is applied to a whole bundle composed of ten tubes. The bundle is simulated under adiabatic conditions, which is a good approximation of the real case of a bundle embedded into the whole plant. Table 1 and Table 2 show respectively the geometrical data of the simulated device and the operating conditions used for the simulation; in addition, the composition of the fuel is ξH2 = 0.72, ξN2 = 0.22, ξCH4 = 0.03 and ξH2O = 0.03 (after humidification).
Figure 7 reports the 2D temperature distributions predicted by the model for the solid structure of the tubes and for the air flowing inside the cathodic ducts. In the visualisation of the results, no interpolation tools have been used, so that the integration elements are visible. Only the temperature distributions of the first, the second and the last tube of the bundle are reported. In the first tube, the endothermic effect of the MSR reaction is evident in the first part of the tube and is the cause of a localised temperature drop occurring immediately after the inlet of the feeding fuel into the tube (with a minimum close to the right-bottom corner). Similar effects are reported also by other studies [35,36,37,38]. It is well known [39,40] that localised temperature drops can cause mechanical stresses in the HT-SOFC structure, which can lead to breakages; for HT-SOFCs, the “rule of thumb” is that a CH4 content around 5% should cause temperature gradients of about 10 °C/cm, which is considered to be at the borderline of the safety window. This is in line with our simulation results, which show that with a CH4 content of 3%, the maximum temperature gradient is about 6 °C/cm. Considering that in our simulations, we assume the MSR reaction to be at equilibrium, which might be excessive (as discussed in Section 4.1.2 “Further Validations”), then our temperature gradients might be over-estimated, and thus, we can certainly conclude that the operating condition under analysis can be considered safe. In Figure 7, in the first tube, when the methane present in the fuel is completely consumed, the cooling effect ceases, and the temperature begins to increase due to the exothermic effect of the electrochemical reaction. It can be noticed that in the first tube, the solid temperature distribution is practically one-directional in the fuel flow direction. Concerning the subsequent tubes of the bundle, the temperature maps regarding the second and the tenth tubes are quite similar to each other, and also, the temperature maps of all of the tubes between the second and the tenth (not reported here) do not manifest significant differences. All of these tubes have solid temperature maps, which are practically one-directional in the air flow direction, which demonstrates that the convective air-solid heat exchange is a major term in the heat balance of the solid. Indeed, in these tubes, fresh air enters from the left-hand side at 910 °C, and then, its temperature increases along the air flow direction due to the heat exchange with the solid. In this way, the heat released by the electrochemical reaction is removed continuously from the solid by the air flow. Figure 7 shows that the temperature of the air, after being fed at 910 °C, increases, but is always about 10–15 degrees below the solid temperature. The fuel temperature (not reported here) is practically identical to that of the solid, which demonstrates that the fuel-solid convective heat transfer is very efficient.
Figure 7. 2D temperature distribution of the solid structure and of the air (for the first, second and 10th tube of the bundle).
Figure 7. 2D temperature distribution of the solid structure and of the air (for the first, second and 10th tube of the bundle).
Energies 08 12364 g007
Figure 8, finally, shows the relationship between solid temperature and current density distributions over the single cells. According to Hypothesis 6, the current density supplied by the electrochemical reaction is zero on the inactive borders of the tubes and also in the interconnections. From Figure 8, it is possible to observe that the regions of the tubes characterised by high temperatures are also characterised by high values of current density. The maximum of current density displayed by the first tube at the left-bottom corner is due to two reasons, i.e., (i) the simultaneous presence of both freshly-supplied fuel and air in this area and (ii) since all of the cells supply the same electrical current (Hypothesis 6) and also as a consequence of the hypothesis of uniform voltage on the cell plane (Hypothesis 7), the low value of current density caused by the localised decreased temperature due to the internal MSR reaction occurring close to the right-bottom corner forces the opposite side, i.e., the left bottom-corner, to supply a high electrical current density. In any case, the maximum departure from the average value of electrical current density (280 mA/cm2) is less than 23%, which appears to be an acceptable value. As a consequence of these unbalanced currents, in spite of the high Nernst potential due to the freshly-supplied fuel and air in this area, potential losses are high, and the first tube performs an overall voltage of 11.9 V (per side), which is lower than the second tube (12.11 V per side). As already discussed, the second tube and the subsequent ones display a much smoother solid temperature distribution, practically one-directional along the air flow direction, and the same distribution is visible also in the maps of electrical current density reported in Figure 8, which, again, are clearly one-dimensional in the air-flow direction. In spite of the more uniform current distribution, the voltage of the subsequent tubes decreases due to the hydrogen consumption in the fuel, and thus, the second tube supplies 11.9 V, the fifth tube 11.5 V and the tenth tube 10.5 V (per side).
Figure 8. 2D distribution of solid temperature and current density (for the first, second and 10th tube of the bundle).
Figure 8. 2D distribution of solid temperature and current density (for the first, second and 10th tube of the bundle).
Energies 08 12364 g008

5. Conclusions

A steady-state model is presented for the simulation of the IP-SOFC tube and bundle currently developed at Rolls-Royce Fuel Cell Systems Ltd. A comparison between simulated and experimental I-V curves obtained from single tubes operated at isothermal conditions shows very good agreement for a number of operating conditions, with various N2/H2, H2/CO/CO2 and also CH4/H2/CO/CO2 humidified mixtures at the anode side and also at different operating temperatures. The comparison between simulation results and experimental data validates the choice of expressing the anodic exchange current density io,a through an equation where io,a is directly proportional to both pH and pH2O.
The simulation tool evaluates I-V curves and also the distribution of the physico-chemical parameters of the full-sized IP-SOFC bundle formed by a number of single tubes and operated at adiabatic conditions, where, overall, the bundle does not exchange heat with the surrounding environment. This condition is approached when the bundle is operated in the plant, surrounded by a thick insulation layer. The simulation results confirm that also for IP-SOFCs, as well as for other HT-SOFC concepts, attention must be given to keeping the CH4 content of the feeding fuel at low levels. When CH4 is present in the fuel feeding mixture, a steep temperature drop occurs in the area close to the edge where the fuel is fed into the tube, due to the internal endothermic MSR reaction. This might be dangerous for the mechanical integrity of the solid structure. Another result is that, when CH4 is present in the fuel feeding mixture, the first tube of the bundle has a temperature distribution that is practically one-directional along the fuel flow direction, while all of the subsequent tubes have temperature distributions very close to each other and practically one-directional along the air-flow direction. The distribution of electrical current density supplied by the single cells deposited on each tube of the bundle closely follows the distribution of the solid temperature.

Notation

Symbols
bgas channel thickness (m)
Bratio gas-solid heat exchange area/cell area (-)
cpmolar specific heat (J·mol−1·K−1)
d1, d2distances between centres of adjacent mesh elements (m)
dhhydraulic diameter (m)
dttime step (s)
Eactactivation energy (J·mol−1)
FFaraday constant (A·s·mol−1)
ΔGGibbs free energy change of the overall reaction (J·mol−1)
ΔG°standard Gibbs free energy change of the overall reaction (J·mol−1)
hheat transfer coefficient (W·m−2·K−1)
ΔHenthalpy change of the overall reaction (J·mol−1)
i0exchange current density (A·m−2)
kelectrical current density (A·m−2)
kthermal conductivity (W·m−1·K−1)
Kthermal conductivity of layered solid (W·m−1·K−1)
Lwidth of stripe of tube fed by one anodic gas channel (m)
nchnumber of channels (-)
nenumber of electrons transferred in the electrochemical reaction
NuNusselt number (-)
ppartial pressure (Pa)
Rggas constant (J·mol−1·K)
Rththermal resistivity of layered solid (W−1·m·K)
rreaction rate (mol·m−3·s−1)
soverall solid thickness (m)
tlayer thickness (m)
Ttemperature (K)
Velectrical potential (V)
Wmolar flow rate (mol·m−2·s−1)
x, ycoordinates (m)
Greek letters
εporosity (-)
γpre-exponential factor (A·m−2)
ηpotential loss (V)
λgas heat conductivity (W·m−1·K−1)
νstoichiometric coefficient (-)
ξmole fraction (-)
Subscript
aanode
actactivation
ccathode
concconcentration
effeffective
ii-th component
jj-th reaction
NernstNernst
ohmohmic
ssolid
stsupporting tube
Superscript
°bulk flow of gases
MSRmethane steam reforming reaction
WGSwater gas shifting reaction

Acknowledgments

Many thanks to Gerry Agnew and Robert Collins of Rolls-Royce Fuel Cell Systems for fruitful discussions during the work and for supplying the experimental data. The authors Paola Costamagna, Simone Grosso and Loredana Magistri gratefully acknowledge financial support from Rolls-Royce Fuel Cell Systems.

Author Contributions

All authors contributed equally to this work.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Singhal, S.C.; Kendall, K. High Temperature Solid Oxide Fuel Cells: Fundamentals, Design and Applications; Elsevier: Philadelphia, PA, USA, 2003. [Google Scholar]
  2. Srinivasan, S. Fuel Cells: From Fundamentals to Applications; Springer: Berlin, Germany, 2006. [Google Scholar]
  3. Brett, D.J.L.; Atkinson, A.; Brandon, N.P.; Skinner, S.J. Intermediate temperature solid oxide fuel cells. Chem. Soc. Rev. 2008, 37, 1568–1578. [Google Scholar] [CrossRef] [PubMed]
  4. Tietz, F. Solid Oxide Fuel Cells, Encyclopedia of Materials: Science and Technology; Elsevier: Philadelphia, PA, USA, 2008; pp. 1–8. [Google Scholar]
  5. Appleby, A.J. Fuel Cells—Overview, in Encyclopedia of Electrochemical Power Sources; Elsevier: Philadelphia, PA, USA, 2009; pp. 277–296. [Google Scholar]
  6. Wu, J.; Liu, X. Recent Development of SOFC metallic interconnect. J. Mater. Sci. Technol. 2010, 26, 293–305. [Google Scholar] [CrossRef]
  7. Behling, N. Fuel Cells: Current Technology Challenges and Future Research Need; Elsevier: Philadelphia, PA, USA, 2012. [Google Scholar]
  8. Huang, K.; Singhal, S.C. Cathode-supported tubular solid oxide fuel cell technology: A critical review. J. Power Sources 2013, 237, 84–97. [Google Scholar] [CrossRef]
  9. Gardner, F.J.; Day, M.J.; Brandon, N.P.; Pashley, M.N.; Cassidy, M. SOFC technology development at Rolls-Royce. J. Power Sources 2000, 86, 122–129. [Google Scholar] [CrossRef]
  10. Debenedetti, P.G.; Vayenas, C.G. Steady state analysis of high temperature fuel cells. Chem. Eng. Sci. 1983, 38, 1817–1829. [Google Scholar] [CrossRef]
  11. Vayenas, C.G.; Debenedetti, P.G.; Yentekakis, Y.; Hegedus, L.L. Cross-flow solid-state electrochemical reactors: A steady-state analysis. Ind. Eng. Chem. Fundam. 1985, 24, 316–324. [Google Scholar] [CrossRef]
  12. Bove, R.; Ubertini, S. (Eds.) Modeling Solid Oxide Fuel Cells: Methods, Procedures and Techniques; Springer: Berlin, Germany, 2008.
  13. Yuan, M.A.J.; Sudén, B. Review on modeling development for multiscale chemical reactions coupled transport phenomena in solid oxide fuel cells. Appl. Energy 2010, 87, 1461–1476. [Google Scholar] [CrossRef]
  14. Hajimolana, S.A.; Hussain, M.A.; Daud, W.M.A.W.; Soroush, M.; Shamiri, A. Mathematical modeling of solid oxide fuel cells: A review. Renew. Sustain Energy Rev. 2011, 15, 1893–1917. [Google Scholar] [CrossRef]
  15. Wang, K.; Hissel, D.; Péra, M.C.; Steiner, N.; Marra, D.; Sorrentino, M.; Pianese, C.; Monteverde, M.; Cardone, P.; Saarinen, J. A Review on solid oxide fuel cell models. Int. J. Hydrogen Energy 2011, 36, 7212–7228. [Google Scholar] [CrossRef]
  16. Grew, K.N.; Chiu, W.K.S. A review on modeling and simulation techniques across the length scales for the solid oxide fuel cell. J. Power Sources 2012, 199, 1–13. [Google Scholar] [CrossRef]
  17. Bhattacharyya, D.; Rengaswamy, R. A Review of Solid Oxide Fuel Cell (SOFC) Dynamic Models. Ind. Eng. Chem. Res. 2009, 48, 6068–6086. [Google Scholar] [CrossRef]
  18. Costamagna, P.; Selimovic, A.; del Borghi, M.; Agnew, G. Electrochemical Model of the Integrated Planar Solid Oxide Fuel Cell (IP-SOFC). Chem. Eng. J. 2004, 102, 61–69. [Google Scholar] [CrossRef]
  19. Cui, D.; Cheng, M. Design for segmented-in-series solid oxide fuel cell through mathematical modelling. J. Power Sources 2010, 195, 1435–1440. [Google Scholar] [CrossRef]
  20. Haberman, B.A.; Young, J.B. A detailed three-dimensional simulation of an IP-SOFC stack. J. Fuel Cell Sci. Tech. 2008, 5. [Google Scholar] [CrossRef]
  21. Mounir, H.; El Gharad, A.; Belaiche, M.; Boukalouch, M. Thermo-fluid and electrochemical modeling of a multi-bundle IP-SOFC—Technology for second generation hybrid application. Energy Convers. Manag. 2009, 50, 2685–2692. [Google Scholar] [CrossRef]
  22. Magistri, L.; Costamagna, P.; Massardo, A.F. A hybrid system based on a personal turbine (5 kW) and a solid oxide fuel cell stack: A flexible and high efficiency energy concept for the distributed power market. J. Eng. Gas Turbine Power 2002, 124, 850–857. [Google Scholar] [CrossRef]
  23. Magistri, L.; Traverso, A.; Cerutti, F.; Bozzolo, M.; Costamagna, P.; Massardo, A.F. Modelling of pressurised hybrid systems based on integrated planar solid oxide fuel cell (IP-SOFC) technology. Fuel Cells 2005, 5, 80–96. [Google Scholar] [CrossRef]
  24. Trasino, F.; Bozzolo, M.; Magistri, L.; Massardo, A.F. Modeling and performance analysis of the Rolls-Royce fuel cell systems limited: 1 MW plant. J. Eng. Gas Turbine Power 2011, 133. [Google Scholar] [CrossRef]
  25. Costamagna, P.; Arato, E.; Antonucci, P.L.; Antonucci, V. Partial oxidation of CH4 in solid oxide fuel cells: Simulation model of the electrochemical reactor and experimental validation. Chem. Eng. Sci. 1996, 51, 3013–3018. [Google Scholar] [CrossRef]
  26. Yurkiv, V. Reformate-operated SOFC anode performance and degradation considering solid carbon formation: A modeling and simulation study. Electrochim. Acta 2014, 143, 114–128. [Google Scholar] [CrossRef]
  27. Day, T.; Singdeo, D.; Bose, M.; Basu, R.N.; Ghosh, P.C. Study of contact resistance at the electrode-interconnect interfaces in planar type solid oxide fuel cells. J. Power Sources 2013, 233, 290–298. [Google Scholar] [CrossRef]
  28. Costamagna, P.; Honegger, K. Modeling of solid oxide heat exchanger integrated stacks and simulation at high fuel utilization. J. Electrochem. Soc. 1998, 145, 3995–4007. [Google Scholar] [CrossRef]
  29. Noren, D.A.; Hoffman, M.A. Clarifying the butler-volmer equation and related approximations for calculating activation losses in solid oxide fuel cell models. J. Power Sources 2005, 152, 175–181. [Google Scholar] [CrossRef]
  30. Wang, Z.; Wang, Z.; Yang, W.; Peng, R.; Lu, Y. Carbon-tolerant solid oxide fuel cells using NiTiO3 as an anode internal reformer. J. Power Sources 2014, 255, 404–409. [Google Scholar] [CrossRef]
  31. Bossel, U.G. Facts and Figures: IEA Task Report; Swiss Federal Office of Energy: Berne, Switzerland, 1992.
  32. Perry, R.H.; Green, D.W. Perry’s Chemical Engineers’ Handbook, 8th ed.; McGraw-Hill: New York, NY, USA, 2007. [Google Scholar]
  33. Mogensen, D.; Grunwaldt, J.-D.; Hendriksen, P.V.; Dam-Johansen, K.; Nielsen, J.U. Internal steam reforming in solid oxide fuel cells: Status and opportunities of kinetic studies and their impact on modelling. J. Power Sources 2011, 196, 25–38. [Google Scholar] [CrossRef]
  34. Tseronis, K.; Kookos, I.K.; Theodoropoulos, C. Modelling mass transport in solid oxide fuel cell anodes: A case for a multidimensional dusty gas-based model. Chem. Eng. Sci. 2008, 63, 5626–5638. [Google Scholar] [CrossRef]
  35. Paradis, H.; Andersson, M.; Yuan, J.; Sundén, B. CFD modeling: Different kinetic approaches for internal reforming reactions in an anode-sopported SOFC. J. Fuel Cell Sci. Technol. 2011, 8. [Google Scholar] [CrossRef]
  36. Achenbach, E.; Riensche, E. Methane/steam redorming kinetics for solid oxide fuel cells. J. Power Sources 1994, 52, 283–288. [Google Scholar] [CrossRef]
  37. Janardhanan, V.M.; Heuveline, V.; Deutschmann, O. Performance analysis of a SOFC under direct internal reforming conditions. J. Power Sources 2007, 172, 296–307. [Google Scholar] [CrossRef]
  38. Aguiar, P.; Adjiman, C.S.; Brandon, N.P. Anode-supported intermediate temperature direct internal reforming solid oxide fuel cell. I: Model-based steady-state performance. J. Power Sources 2004, 138, 120–136. [Google Scholar] [CrossRef]
  39. Clague, R.; Marquis, A.J.; Brandon, N.P. Finite element and analytical stress analysis of a solid oxide fuel cell. J. Power Sources 2012, 210, 224–232. [Google Scholar] [CrossRef]
  40. Clague, R.; Marquis, A.J.; Brandon, N.P. Time independent and time dependent probability of failure of solid oxide fuel cells by stress analysis and the Weibull method. J. Power Sources 2013, 221, 290–299. [Google Scholar] [CrossRef]

Share and Cite

MDPI and ACS Style

Costamagna, P.; Grosso, S.; Travis, R.; Magistri, L. Integrated Planar Solid Oxide Fuel Cell: Steady-State Model of a Bundle and Validation through Single Tube Experimental Data. Energies 2015, 8, 13231-13254. https://doi.org/10.3390/en81112364

AMA Style

Costamagna P, Grosso S, Travis R, Magistri L. Integrated Planar Solid Oxide Fuel Cell: Steady-State Model of a Bundle and Validation through Single Tube Experimental Data. Energies. 2015; 8(11):13231-13254. https://doi.org/10.3390/en81112364

Chicago/Turabian Style

Costamagna, Paola, Simone Grosso, Rowland Travis, and Loredana Magistri. 2015. "Integrated Planar Solid Oxide Fuel Cell: Steady-State Model of a Bundle and Validation through Single Tube Experimental Data" Energies 8, no. 11: 13231-13254. https://doi.org/10.3390/en81112364

APA Style

Costamagna, P., Grosso, S., Travis, R., & Magistri, L. (2015). Integrated Planar Solid Oxide Fuel Cell: Steady-State Model of a Bundle and Validation through Single Tube Experimental Data. Energies, 8(11), 13231-13254. https://doi.org/10.3390/en81112364

Article Metrics

Back to TopTop