Next Article in Journal
Online Adaptive Parameter Estimation of a Finite Control Set Model Predictive Controlled Hybrid Active Power Filter
Next Article in Special Issue
A System-Level Modeling of PEMFC Considering Degradation Aspect towards a Diagnosis Process
Previous Article in Journal
Implementation and Evaluation of a Complex Pumped-Storage Hydropower Plant with Four Units, Common Penstock, and Surge Tank in a Real-Time Digital Simulator
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Numerical Modeling and Simulation of the Solid Oxide Cell Stacks and Metal Interconnect Oxidation with OpenFOAM

1
Institute of Energy and Climate Research, Fundamental Electrochemistry (IEK-9), Forschungszentrum Jülich GmbH, D-52425 Jülich, Germany
2
Institute of Physical Chemistry, RWTH Aachen University, D-52074 Aachen, Germany
*
Author to whom correspondence should be addressed.
Energies 2023, 16(9), 3827; https://doi.org/10.3390/en16093827
Submission received: 3 April 2023 / Revised: 24 April 2023 / Accepted: 27 April 2023 / Published: 29 April 2023
(This article belongs to the Special Issue Hydrogen and Fuel Cell Technology, Modelling and Simulation II)

Abstract

:
Solid oxide cells are capable of efficiently converting various chemical energy carriers to electricity and vice versa. The urgent challenge nowadays is the faster degradation rate compared with other fuel cell/electrolyzer technologies. To understand the degradation mechanisms, simulation of a solid oxide cell is helpful. Since most previous research developed models using commercial software, such as COMSOL and ANSYS Fluent, a gap for knowledge transfer is being gradually formed between academia and industry due to licensing issues. This paper introduces a multiphysics model, developed by a computational code, openFuelCell2. The code is implemented with an open-source library, OpenFOAM. It accounts for momentum transfer, mass transfer, electrochemical reactions and metal interconnect oxidation. The model can precisely predict I–V curves under different temperatures, fuel humidity and operation modes. Comparison between OpenFOAM and COMSOL simulations shows good agreement. The metal interconnect oxidation is modeled, which can predict the thickness of the oxide scale under different protective coatings. Simulations are conducted by assuming an ultra-thin film resistance on the rib surface. It is revealed that coatings fabricated by atmospheric plasma spraying can efficiently prevent metal interconnect oxidation, with a contribution of only 0.53 % to the total degradation rate.

1. Introduction

The transition to a renewable energy system where solar and wind energy prevail leads to a challenge for the grid balance. Applying solid oxide cells (SOCs) as energy storage and conversion devices may be part of the solution to these problems. Their operation at high temperatures provides several advantages, such as high efficiencies, the possibility of sector coupling on the heat side, the capability of being reversibly operated, and the tolerance to various fuels. In order to increase the lifetime of SOCs, their degradation behavior should be studied, as it obstructs the long-term operation. Significant degradation mechanisms, such as Ni agglomeration and Cr poisoning, are related to the local distribution of the overpotential and chemical species. However, measuring local values for these during the operation of SOCs is not only a technical challenge but also a potentially dangerous task, as gas leakages, burns and electric shocks could happen. In contrast, simulation can be an efficient and safe method to visualize the local physical values. Therefore, a three-dimension model that can capture the physical processes in an SOC is necessary.
Although there have been a number of scientific findings due to simulation in the past decade [1,2,3,4,5,6,7,8], simulations are usually steady-state simulations, making it impossible to trace the degradation over time. It can be seen that many numerical investigations of SOCs were carried out with commercial software, for example, DETCHEM, ANSYS Fluent, COMSOL, etc. Mathematical models used to describe the multiphysical and multiscale transport processes in the SOCs have been incorporated into a ‘black-box’. It is usually impossible to visualize the implementation of these models. In addition, the knowledge transfer may be limited within the research community due to the license issues. In this work, the authors employed a computational code, openFuelCell2 [9,10], which was implemented with an open-source library, OpenFOAM [11]. The code witnessed reliable performance in low-temperature electrochemical applications, for example, the proton exchange membrane fuel cells [12,13,14,15]. This work aims to extends the capability of openFuelCell2 to high-temperature applications, i.e., SOFC and SOEC. The comparison between the simulation results with openFuelCell2 and COMSOL is conducted for the basic steady-state model. A degradation model to describe the performance degradation with the evolution of time due to metal interconnect (MIC) oxide scale growth is additionally proposed and implemented into openFuelCell2. The effects of different protective coatings on MIC oxidation are discussed. Numerical simulations are performed for the F10 SOC stack design [16] in Forschungszentrum Jülich GmbH (FZJ), the detailed description of which is given in the next section.

2. Geometry of the Model

It is helpful to introduce the general information of the F10 SOC stack design in FZJ first. The simplified geometry of the repeating unit of the SOC stack is shown at the left-hand side of Figure 1. The standard design is an fuel-electrode supported cell (FESC), where a thick support layer at the fuel side, made of Ni/8YSZ (8 mol% yttria-stabilized zirconia), provides necessary mechanical strength. The repeating unit consists of roughly 1000 µm thick Ni-mesh as the distributor and contacting element at the fuel side, a 500 µm thick support layer Ni/8YSZ [17], a 7 µm thick fuel electrode Ni/8YSZ [17], a 10 µm thick 8YSZ electrolyte [17], a 45 µm thick air electrode La0.58Sr0.4Co0.2Fe0.8O3-δ (LSCF) [17], and a 100 µm thick contact layer LSCF (or LCC12, La0.97Mn0.4Co0.3Cu0.3O3-δ). A barrier layer (GDC, gadolinia-doped ceria) is added between the electrolyte and the air electrode to prevent interfacial reaction. Furthermore, there is a protective coating on the MIC to alleviate issues of metal oxidation and Cr poisoning. In FZJ, the protective coating on the MIC surface is usually made of M n C o 1.9 F e 0.1 O 4 ( M C F ) applied by atmospheric plasma spraying (APS). Another widely used coating is made of MnOx and is added onto MIC via wet powder spraying (WPS). The drawback of WPS, in contrast to APS, is that it can only generate a porous layer that cannot efficiently prevent metal oxidation and chromium poisoning. Consequently, stacks with APS coatings have better performance than stacks with WPS coatings [18]. The MIC is made by either Crofer 22 APU or intermediate temperature metal (ITM). For sealing the stack glass, ceramic sealants are used on all bonding surfaces.
The numerical simulation of the whole geometry usually means a rather large calculation effort. With the purpose of accelerating the computation while capturing the main features of the repeating unit of a real SOC stack, the following assumptions are made:
  • The barrier layer and protective coating are ignored in the geometry of the model. Their effects on the resistance are considered in the electronic conductivity of the air electrode and air-side contact layer.
  • The metal frame is ignored in the geometry of the model.
Therefore, the computational domain in this study refers to a representative three-dimensional slice of one channel in a real SOC stack as shown on the right-hand side of Figure 1.

3. Governing Equations of the Basic Model

The SOC stack is a complex electrochemical device that contains different physics processes, which are coupled with each other. This interplay strongly affects the performance of the SOC stacks. Therefore, suitable mathematical models should be implemented; otherwise, the results are questionable. Although it is recommended to involve as many details as possible, several plausible simplifications can be made to alleviate the difficulty in multiphysics modeling:
  • Properties of materials, such as permeability, porosity and electrical conductivity, are homogeneous and isotropic.
  • All fluid is in the laminar flow regime.
  • Heat transfer is ignored. In other words, isothermal assumption is applied in the model. The assumption holds well for stacks listed in this work, as the temperature difference is 10∼15 C according to thermocouple measurements.
  • The quality of the inlet gas is perfect, meaning the air only has oxygen and nitrogen, and the fuel only has hydrogen and steam.
In this section, models of major physical processes, including electrochemical reaction, momentum transfer and mass transport, are present. A summary of the models used in this work and corresponding variables to be solved is shown in Table 1.

3.1. Electrochemical Reaction

Electrochemical reaction is a key part in modeling SOCs due to its direct relation to I–V curves and strong coupling with other physics. Commonly, it is simulated by the Butler–Volmer equation in Equation (1) [19]:
i = i 0 exp β a n e F R T η exp β c n e F R T η
where i is the current density, i 0 is the exchange current density, β a ( β c ) is the anodic (cathodic) charge transfer coefficient, n e is the number of transferred electrons, F is the Faraday constant, R is the gas constant, T is the temperature, and η is the overpotential. In spite of the popularity of Equation (1), Mogensen et al. [20] raised several arguments against using Equation (1) for SOCs, such as the reaction rate probably not being limited by the charge transfer. Similarly, Adler et al. [21] found that, except for charge transfer, the surface diffusion of chemical species and dissociative chemisorption were also important for reaction rates. As illustrated by Maria et al. [22], it could be that the high operation temperature of SOC makes adsorption, diffusion and other processes more rate-limiting than the charge transfer. Additionally, in most cases, linear I–V curves instead of exponentially changed I–V curves are observed from the measurement in FZJ.
In light of the above reasons, in this work, a linear expression of electrochemical kinetics is used, shown in Equations (2) and (3) [22,23]:
i v , fuel = η fuel R fuel exp ( E a , fuel R T ) a H 2 0.1 a H 2 O 0.33
i v , air = η air R air exp ( E a , air R T ) a O 2 0.25
where i v is the volumetric reaction current density, η is the overpotential, R is a pre-factor, E a is the energy barrier, and a is the activity of the reactant. The overpotential can be calculated by [24]
η air = ϕ ele , air ϕ io , air ϕ eq , air
η fuel = ϕ ele , fuel ϕ io , fuel ϕ eq , fuel
where ϕ io is the ionic potential, ϕ ele is the electromotive potential and ϕ eq is the equilibrium potential defined as [22]
ϕ eq , air = R T 4 F ln ( Y O 2 )
ϕ eq , fuel = Δ G water splitting 2 F R T 2 F ln Y H 2 Y H 2 O
where Y O 2 , Y H 2 and Y H 2 O are molar fractions of oxygen, hydrogen and water, respectively. Δ G water splitting is the Gibbs free energy change of water splitting.
Two types of current exist in SOC: ionic current i io and electronic current i ele . The current and potential obey Ohm’s law:
i io = σ io ϕ io
i ele = σ ele ϕ ele
where σ io and σ ele are the ionic conductivity and electronic conductivity that can be found in Table 2. In a porous and composite medium, the conductivity needs to be corrected by considering the effect of the solid volume fraction χ of and tortuosity τ of different materials. If the number of different material properties in the medium is n, the effective conductivity is [25]
σ ele,eff = i = 1 n = σ i ele χ i τ i
σ io,eff = i = 1 n = σ i io χ i τ i
The conservation of current is guaranteed by
i io = S current
i ele = S current
where S current is a current source defined as the following:
S current = i v in electrochemical active regions 0 in electrochemical inactive regions
and i v is given in Equations (2) and (3).
The electronic and ionic conductivity is present in Table 2. Other parameters, such as the porosity and volume factions of the porous medium, are given in Table A1.

3.2. Momentum Transfer

Momentum transfer of fluid happens in gas channels and porous medium. In gas channels, it is a common practice to apply Navier–Stokes (NS) equations to describe the momentum transfer as shown in Equation (15) [28]:
( ρ f U ) t + ρ f ( U · ) U = ( P I ) + ρ f g + μ U + ( U ) T 2 μ 3 ( · U ) I
where ρ f is the density of the fluid, P is the pressure, U is the velocity, μ is the viscosity, and I is the identity matrix.
In a porous medium, momentum transfer can be simulated by several models. The easiest one is Darcy’s law, considering that the velocity of fluid should linearly respond to the pressure gradient:
U = κ μ P
where κ is the permeability of the porous medium. In spite of its conciseness, it should be kept in mind that it is an empirical relation that originated from the experimental data of water passing through sands at low velocity [29]. Therefore, there are some limitations to Equation (16). For example, Darcy’s law neglects viscous forces. Assuming that Darcy’s law is implemented in a porous medium and the NS equation is used in gas channels, two calculation domains need to be defined because Equation (15) has second-order differential terms, while Equation (16) misses it. This inherent conflict leads to trouble in setting boundary conditions at the interface between channels and the porous medium [30].
Provided by the drawbacks in Darcy’s law, it is recommended to use the Darcy’s modified NS equation, which is able to be implemented in both gas channels and porous media, as given in Equation (17) [31,32]:
( ρ f U / ε ) t + ρ f ( U ε · ) U ε = ( P I ) + ρ f g + μ ε U + ( U ) T 2 μ 3 ε ( · U ) I μ κ U
where ε is the porosity of the porous medium. In gas channels where permeability approaches infinity and porosity becomes one, Equation (17) can automatically become Equation (15), which avoids setting another boundary condition at the interface between the gas channels and porous medium. The continuity equation is
( ε ρ f ) t + ( ρ f U ) = S mass
where S mass is the mass source due to reactions.
The data of porosity and permeability are given in Table A1.

3.3. Mass Transport

Mass transport is governed by Equation (19):
( c Y i ) t + ( c U Y i ) + N i = S mole , i
where c is the molar density, Y i is the molar fraction, N i is the molar diffusion flux of species i and S mole , i is the source term of moles of chemical specie i due to electrochemcial reactions, and it is described by the following Equation (20):
S mole , i = i v n e F in electrochemical active regions 0 in electrochemical inactive regions
where i v is the reaction current density, and n e is the number of electrons transferred in electrochemical reactions. To solve Equation (19), N i needs to be calculated. Several diffusion models can be chosen. For instance, Fick’s law can be used, which is the most popular model due to its simplicity, and it has been successfully implemented in SOC numerical modeling [33,34]. The mathematical description is shown in Equation (21)
N i = D i j B ( c Y i )
where D i j B is the binary diffusion coefficient for species i and j. For a multi-component system, it is more plausible to apply the Stefan–Maxwell equation to describe the transport phenomena, which is the extension of Fick’s law. The Stefan–Maxwell equation is given in Equation (22)
j = 1 , j i n Y j N i Y i N j D i j B , eff = P R T Y i Y i R T P
where D i j B , eff is the effective binary diffusion coefficient. It should be noted that both Fick’s law and the Stefan–Maxwell equation are factually applied for the open space instead of the porous medium [35]. To simulate the behavior in the porous medium, effective diffusion coefficients are introduced, which are related to the porosity and tortuosity of the porous medium. In the model, Equation (23) is used to define the effective diffusion coefficient [25]:
D i j B , eff = ε k τ k D i j B
where ε is the porosity and τ is the tortuosity. k is the component, such as the fuel electrode and the air electrode. The tortuosity is calculated through following relationship:
τ k = ( 1 ε k ) 2 2 D s 2 D s
where D s is usually chosen as 1.2 for SOC systems [36]. The binary diffusion coefficient D i j B is evaluated by the Chapman–Enskog theory [37]:
D i j B = D j i B = 3.16 × 10 8 T 1.75 P ( υ i 1 / 3 + υ j 1 / 3 ) 2 ( 1 M i + 1 M j ) 1 / 2
where υ is the diffusion volume and M is the molar mass.
Fick’s law and the Stefan–Maxwell model function well regarding mass transfer in a porous medium, but they will fail if the interaction between the species and solids cannot be ignored. In this case, Knudsen diffusion should be taken into consideration in addition to molecular diffusion. One model that includes Knudsen diffusion is the dusty gas model as given in Equation (26):
j = 1 , j i n Y j N i Y i N j D i j B , eff + N i D i K , eff = P R T Y i Y i R T 1 + κ P μ D i K , eff P
where κ is the permeability, D i K , eff is the effective Knudsen diffusion coefficient for the species i, and it can be calculated through Equation (27) [25]:
D i K , eff = ε k τ k D i K = d p ε k 3 τ k 8 R T π M i
where d p is the diameter of the pore size.
In brief, three models, Fick’s law, the Stefan–Maxwell model and the dusty gas model, can be chosen to simulate the mass transport. In this paper, Fick’s law is chosen after several considerations. First, it is easiest to achieve numerically, as it has a simpler form than the others. Secondly, Fick’s law actually can approximate the other two models well, particularly when the number of components in the system is less than four [38]. Thirdly, to improve the accuracy, Fick’s law can be modified to include the Knudsen diffusion [22,39,40] by employing the Bosanquet formula as shown in Equation (28):
1 D i j eff = 1 D i j B , eff + 1 D i j K , eff
where D i j K , eff is a mean Knudsen diffusion coefficient that is obtained by replacing M i with M i j in Equation (27), and M i j is defined as
1 M i j = 1 M i + 1 M j / 2
Given the above, Fick’s law (Equations (21), (27)–(29)) is applied in the model.

4. Governing Equations of the Degradation Model of MIC Oxidation

The typical material used for interconnects of SOC stacks are metals (e.g., Crofer 22 APU). It provides strong mechanical support and high electronic conductivity. However, it always suffers from oxidation [41], leading to the formation of the oxide scale. During the growth of the oxide scale, the resistance at the interface increases, which serves as one part of the total degradation. The model in this work only considers the MIC oxidation at the air side. On the fuel side, MIC oxidation is believed to have little influence on the performance of SOCs. The Ni mesh results in good contact with the fuel support layer, providing spot welding is used to connect MIC and Ni mesh, which provides a strong and reliable contact.

4.1. Thickness of the Oxide Scale on a Clean Surface

To simulate the growth of the oxide scale, it is a common practice to employ the parabolic law [42,43] given in Equation (30)
d 2 = Ξ · t
where d is the thickness of the oxide scale, t is time, and Ξ is the parabolic constant. The “constant” means that its value does not depend on the thickness of the oxide scale. Before applying Equation (30), it is necessary to prove that, for SOCs, the growth of the C r 2 O 3 fulfills the parabolic format. It was pointed out that the growth of C r 2 O 3 obeys the parabolic law when the temperature is above 623 K [44], which is lower than the operating temperature of SOCs. There are works in the literature [42,45], where the parabolic law was used to ensure good curve fitting. In brief, it is reasonable to assume that Equation (30) is satisfied during the whole operation of SOCs.
The next step is to obtain the parabolic constant Ξ . The typical method is Wagner’s theory, which relates ions’ tracer diffusion coefficients to the parabolic constant. According to Wagner’s theory, the parabolic constant can be written as follows:
Ξ = Y O 2 metal oxide Y O 2 oxide gas Z C r 3 + / Z O 2 D C r 3 + eff + D O 2 eff d ln Y O 2
where Z is the absolute value of the valency, D eff is the effective tracer diffusion coefficients, which should include the bulk diffusion and grain-boundary diffusion [46,47], the superscript oxide–gas and metal–oxide refer to the interface of oxide–gas and the interface of metal–oxide, respectively. Although Equation (31) is widely used for the simulation of alloy oxidation, it is not suitable for the model of SOCs due to the following reasons:
  • In order to obtain tracer diffusion coefficients, the concentration of defects (e.g., V C r that serve as the path for the ions needs to be known. However, there remain debates [47, 48, 49] as to which defects dominate inside Cr2O3.
  • No literature about tracer diffusion coefficients inside the oxide scale on Crofer 22 APU can be found. If the data of a different alloy are used (e.g., Ni-Cr alloy [49]),the parabolic constant obtained from Equation (31) will possibly deviate from the experimental findings.
In light of the difficulties in applying Equation (31), an empirical relation is used and introduced in the following section, which was proposed by FZJ [42].

4.2. Thickness of the Oxide Scale on the Surface Covered by a Porous Coating

It should be noted that Equation (30) only works for a clean surface without coatings. However, the MIC is usually protected by WPS M n O x or APS MCF. Considering the thickness of oxide scale is proportional to the mole of oxygen atoms adsorbed on the MIC, the following equation is obtained:
d poro n O , poro = n ˙ O , poro d A poro d t
= n ˙ O , poro d ε coating A d t
= ε coating n ˙ O , poro d A d t
= ε coating n ˙ O d A d t
= ε coating n O
which means
d poro = ε coating d
where d poro is the thickness of the oxide scale under porous coatings (e.g., the protective coating), n O , poro is the mole of oxygen atoms adsorbed on the MIC surface covered by coatings, n ˙ O , poro is the adsorption rate of oxygen atoms per area over the MIC surface covered by coatings, A poro is the free surface area under a porous coating that can be reached by oxygen gas, A is the total surface area of MIC, ε coating is the porosity of the porous coatings, n ˙ O is the adsorption rate of oxygen atoms over the MIC surface without coatings, and n O is the mole of oxygen atoms adsorbed on the MIC surface without coatings. Two underlying assumptions are used for derivation of Equation (37). First, the homogenization assumption is applied in Equation (). It is usually hard for CFD models to capture micro-structure values. Secondly, it is assumed in Equation () that the adsorption rate over the MIC surface is not affected by the coating. The most evident difference of physical phenomena between the cases with and without coating is the gas transport of oxygen, it means the oxygen partial pressure affects insignificantly on the adsorption rate. If the porosity of the coating is higher than 0.3, this assumption is reasonable, as it was found that oxygen could freely reach the surface [45]. At present, there is no direct evidence to support the second assumption regarding the application of a dense coating. However, it is thought that the rate-limiting step during metal oxidation is usually the diffusion of ions instead of gas transport over the metal surface [43]. Therefore, the second assumption is generally valid for the growth of the oxide scale under a dense coating.
In brief, Equations (30) and (37) are used to simulate the growth of the oxide scale on MIC in SOCs.

4.3. Voltage Drop Due to MIC Oxidation

To simulate the voltage drop due to MIC oxidation, an internal boundary condition is needed. In the model, this boundary is called the ultra-thin film resistance. The effect of the oxide scale is represented by an internal surface boundary condition at the rib, as shown in Figure 2. The mesh of the rib surface is split into two sub-surfaces, and Ohm’s law is applied to each of them as follows:
n up i ele = σ oxide scale d poro ϕ ele , up ϕ ele , low
n low i ele = σ oxide scale d poro ϕ ele , low ϕ ele , up
where n is the unit surface normal vector, i ele is the electronic current density, σ oxide scale is the conductivity of the oxide scale, and ϕ ele is the electronic potential. The subscripts “up” and “low” refer to the upper surface (with higher electronic potential) and lower surface (with lower electronic potential), respectively. The directions of the above vectors can be found in Figure 2. The difference between ϕ ele , up and ϕ ele , low arises from the resistance of the growing oxide scale, which can be obtained through
Δ ϕ ele = ϕ ele , up ϕ ele , low = | i ele d poro σ oxide scale |
where Δ ϕ ele is the voltage degradation due to the MIC oxidation.

4.4. Parameters

According to the experimental data [42], the parabolic constant of Crofer 22 APU (700∼900 C) is summarized in Equation (41):
Ξ Crofer = 5.5075 × 10 9 exp 1.9992 × 10 5 R T
It was found that under 800 C, the growth rate of the oxide scale on ITM was half of that on Crofer 22 APU [50]. Assuming this is valid for the temperature range of 700∼900 C, the parabolic constant for ITM is,
Ξ ITM = 0.25 × Ξ Crofer = 1.3769 × 10 9 exp 1.9992 × 10 5 R T
The conductivity of the oxide scale on Crofer 22 APU, σ oxide scale , Crofer , is summarized as follows with the help of experimental data [43]:
σ oxide scale , Crofer = 1 T 10 3228 / T + 6.7661
The available conductivity data of the oxide scale on ITM are very limited. It was found that the increasing rate of ASR on ITM during oxidation was about 3.5 times as much as that on Crofer 22 APU [51], even if the thickness was half. The explanation was that if the coating is not a spinel phase, the oxide scale generated on ITM is porous, while that on Crofer 22 APU is quite dense [51]. ITM was usually combined with M n O x coating (not a spinel phase) in FZJ SOC stacks design, it is assumed that the conductivity of the oxide scale on ITM, σ oxide scale , ITM , is
σ oxide scale , ITM = 1 7 σ oxide scale , Crofer = 1 7 T 10 3228 / T + 6.7661
The porosity of WPS coating is 0.45 [52], while the porosity of APS coating is 0.03 [53].

5. Numerical Procedure

The mesh in the model is presented in Figure 3, where, in total, 74,800 elements are made. Mesh refinement is performed at the interface to capture electrochemical reactions. In other areas, mesh is coarser.
The basic model is steady state, while MIC oxidation is a time-dependent model. Discretization of the time derivative is achieved by the backward Euler method. To increase the convergence of the time-dependent simulation, the basic model is solved first. The time-dependent simulation of MIC oxidation is carried out with the initial values being results of the basic model. The time step for the time-dependent simulation is 10 h. The discretization of the gradient terms and divergence terms is performed with the Gauss linear method; the discretization of all Laplacian terms is Gauss linear, except the conservation of the current, which is discretized by the Gauss harmonic. In addition, a SIMPLE solution algorithm [54] is employed to solve the Navier–Stokes equation. The coupling between each physics is summarized in Figure 4, where the arrows mean the passed variables to the other physics. For example, to solve the governing equation of the mass transport, pressure and velocity fields are needed, which are taken from the solution of the momentum transfer.
The solver and discretization used in COMSOL and OpenFOAM are different. The former employs a coupled solution manner, which solves all variables simultaneously. During time intervals, a direct solver, which solves the inverse of the matrix, is used. In contrast, the segregated solver is applied in OpenFOAM that subdivides the numerical calculation into several steps. For each iteration, the iterative solver is applied. Moreover, COMSOL implements the finite element method, while OpenFOAM utilizes the finite volume method. Despite the above differences in the numerical calculation, it is expected that both software will give similar results when it comes to the same model with the same boundary conditions and parameters. All the simulations are conducted on a single core of a computer with CPU i7-9700K and 32 GB RAM.

6. Experiment

SOC stacks were operated in a constant current mode with the gases supplied in a counter-flow regime. The characterization and operation of SOC stacks were carried out in an oven test rig. Inside the oven, there were electric heating elements to maintain the temperature of stacks through heat radiation. Deionized water with a resistivity of 10 MΩ cm [55] was supplied to the electrically heated steam generator. Hydrogen with a minimal purity of 99.9% was used for the inlet fuel. Dehumidified air with a dew point below −40 C was sent as the inlet air. The mass flow of inlet air and inlet fuel was controlled by the mass flow controllers (MFCs). Mica sealants between the test bench interface and the stack in combination with a compressive force of 1 kN were applied to prevent gas leakage. Table 3 shows details of the MFCs, operation conditions and stack designs of several stacks used for the validation in this work. The temperature values in Table 3 were measured by thermocouples that were placed roughly 10 mm deep into the intermediate interconnector (details of the locations can be found in reference [56]).
F1004-115 consists of a 400 µm thick support layer, 10 µm thick fuel electrode and a 6 µm thick electrolyte. The geometry of the rest of the parts of F1004-115 and other stacks can be found in Figure 1.

7. Validation

7.1. Validation of I–V Curves of the Basic Model

The validation of the OpenFOAM simulation has two parts: (1) The first shows a comparison between the OpenFOAM simulation and the COMSOL simulation. The same mesh, geometry, model, boundary conditions and parameters are used in both simulations. Although OpenFOAM and COMSOL use different discretization methods and solvers, their results are expected to be close. (2) The second shows a comparison between the OpenFOAM simulation and the experimental results. The comparison is carried out under different operation modes (e.g., fuel cell and electrolysis modes) and various operation conditions (e.g., different operation temperature and humidity of inlet fuel).
The validation by I–V curves is shown in Figure 5. Figure 5a,b present the I–V curves only in the fuel cell mode. With varied humidity and temperature, very good matches between the OpenFOAM simulation and the experimental data can be found, where the maximum voltage deviation is only 4.5%. Figure 5c presents the I–V curves for reversible SOCs. The negative current density means the electrolysis operation, while the positive current density indicates the fuel cell operation. There is good agreements between OpenFOAM and the experimental measurement providing the maximum voltage deviation is just 4.1%. The OpenFOAM simulation closely approximates the COMSOL simulation in all cases shown in Figure 5, which proves an additional validity of the numerical model in OpenFOAM. The deviations between the simulated and experimentally measured I–V curves could be due to the fluctuating temperature (5∼10 C) during the I–V curve characterizations.

7.2. Validation of MIC Degradation Model

The core equation of the MIC oxidation model is Equation (37), which describes the thickness growth of the oxide scale under a porous coating. In the following, examples are given to demonstrate the validity of Equation (37).

7.2.1. Result from Persson: Crofer 22 APU with and without ( L a , S r ) M n O 3 (LSM) Coating

Persson [45] conducted a series of oxidation tests on Crofer 22 APU. The experimental data were taken from thermogravimetric analysis (TGA), which made it possible to track the growth of the oxide scale during operation. Assuming that Equation (37) provides a good approximation of the experimental results, we would expect the ratio of the thickness of the oxide scale with coating to the thickness of the oxide scale without coating to be close to the porosity of the coating, which is 0.3. According to the results in Figure 6, after 1000 hours of operation, the ratio is close to 0.3, which confirms the validity of the model. However, at the beginning of the operation, the ratio is higher than 0.3. This could result from the sintering process of the coating, as Persson used spray guns to add coatings on the MIC at room temperature and then started the oxidation tests directly [45]. However, in general, it can be claimed that the model is capable of predicting an accurate thickness of the oxide scale.

7.2.2. Result from FZJ: F1002-95 and F1002-97

Two stacks designed by FZJ are used for validation, F1002-97 and F1002-95. The former used ITM as the MIC and WPS MnOx as the protective coating, while the latter employed Crofer 22 APU as the MIC and WPS MnOx as the protective coating. Details of the stacks can be found in Table 3.
Unlike the TGA data in Figure 6, the validation data in this section are manually obtained from SEM figures. It should be kept in mind that the model predicted an average thickness. As shown in Figure 7 where average thickness is compared, good agreement is observed, although there exist slight deviations between the simulation and the experimentally measured data. The deviations can be due to the limited choices of parameters (e.g., the porosity of the protective coating and properties of ITM) and the errors when processing data from SEM images.

8. Results and Discussions

8.1. Spatial Distribution of Mole Fraction and Overpotential

Distribution of mole fraction and overpotential can affect the local degradation. Here, only analysis of the SOEC is given.
The spatial distribution of the mole fraction and overpotential are shown in Figure 8 and Figure 9, respectively. The simulation works for the electrolysis operation of F1004-115 under −0.5 A/cm 2 , fuel inlet composition H2/H2O = 50:50 and air inlet composition O2/N2 = 21:79. In the simulation, the flow field is the counter flow, where fuel flows toward the x-coordinate, while air flows against the x-coordinate.
Figure 8 shows the mole fraction distribution in the electrolysis mode of F1004-115, at the interface between the fuel electrode and the electrolyte (FEL/ELEC) and the interface between the air electrode and the electrolyte (AEL/ELEC). From Figure 8a, it is found that along the x-coordinate, which is the fuel flow direction, Y H 2 gradually increases due to the electrochemical reaction. Similarly, as indicated by Figure 8c, Y O 2 increases along the flow direction (against the x-coordinate). In addition, it should be noted that O 2 tends to segregate at the rib zones. As a comparison, on the fuel side, H 2 is more uniformly distributed. The main reason behind this phenomenon is that a gas channel is used for air distribution, while the Ni mesh is used for sending the fuel. Unlike the Ni mesh, the air channel does not directly contact the rib zones. Hence, O 2 transportation at the rib zones mainly relies on the diffusion process. The diffusion coefficient of O 2 is not high, as it is one order of magnitude lower than that of hydrogen and water. Consequently, the lack of enough of a driving force of transportation leads to O 2 segregation at the rib zones. When it comes to the mole fraction, the OpenFOAM simulation is close to the COMSOL simulation. There is a very slight difference between Figure 8a,b, which can be caused by different solvers applied in OpenFOAM and COMSOL.
Figure 9 shows the overpotential distribution in the electrolysis mode of F1004-115, where η fuel is negative and η air is positive. According to Figure 9a,b, both overpotentials η fuel and η air reach the maximum absolute values close to the electrolyte as expected because the overpotential is the driving force for electrochemical reactions, which majorly happen near the electrolyte interface. Compared to the COMSOL simulation (right column), the evolution tendency is the same, while the values are slightly different.

8.2. Voltage Degradation Due to MIC Oxidation in F1002-97

The details of the stack F1002-97 can be found in Table 3. Figure 10 shows stack voltage evolution predicted by simulation at four time points. The inserted figures show the distributions of the electronic potential in the contact layer and the MIC. At time = 0 h, when no oxide scale exists, the electronic potential is uniformly distributed because the electronic conductivity of materials is high. With time going on, the oxide scale gradually grows, leading to degradation. It can be seen that the ultra-thin film resistance leads to a sudden potential drop across the rib in the inserted figures. Figure 11 gives the evolution of the voltage of the whole lifetime and the share of degradation rate only due to MIC oxidation. Because of the unstable water supply, the voltage measured in experiments fluctuate all the time. With the help of linear fitting, the degradation rate is determined to be 7.3 mV/kh before 40 kh and 1.6 mV/kh after 40 kh. Figure 11a indicates that, overall, the MIC oxidation contributed ∼11% of the total degradation, which is not negligible. In addition, as shown in Figure 11b, its influence could be high in the first 5 kh, contributing 20–31% to the total degradation.

8.3. Voltage Degradation Due to MIC Oxidation in F1004-67

Crofer 22 APU was used as the MIC, and APS MCF was applied as the protective coating for the stack F1004-67. It was operated at 730 C and 0.5 A/cm 2 with a degradation rate of 3 mV/kh over 25 kh operation. With a dense coating, the simulation predicts that the voltage decrease at a rate of 0.016 mv/kh due to MIC oxidation, which only contributes ∼0.53% to the total degradation. Therefore, degradation due to the MIC oxidation can be safely neglected. The degradation mainly arise from other processes, such as Ni agglomeration and Cr poisoning.

9. Conclusions and Outlook

A three-dimensional multiphysics SOC model is built by using a computational code, openFuelCell2. The code is constructed inside an open-source library, OpenFOAM. The model is proven to be able to carry out not only steady-state but also time-dependent simulations. The model is intensively validated for the steady-state simulation to precisely predict I–V curves under different operation temperatures, fuel humidities and operation modes with a maximum error below 4.5%. It is indicated that the model is able to predict the SOC performance, as the spatial distribution of overpotential and mole fractions simulated by OpenFOAM and COMSOL are very close. A time-dependent model of MIC oxidation is proposed. Different experiments validated that the model can predict accurate thicknesses of the oxide scale. According to the simulation, it is found that if MCF protective coating fabricated by APS is applied, the degradation due to the MIC oxidation can be generally neglected since it only contributes ∼0.53% to the total degradation rate. In contrast, if WPS MnOx protective coating is used, the MIC oxidation cannot be ignored, especially at the beginning of the operation.
This work establishes a framework for further sophisticated modeling studies. For example, although an empirical formulation is currently used to calculate the parabolic constant, models accounting for more detailed physics can be proposed and implemented. More degradation mechanisms, such as Ni agglomeration and Cr poisoning, will be integrated in the model in the near future to reveal how the performance of the SOC stack evolves with time and how to optimize the stack design.

Author Contributions

Conceptualization, S.Y.; data curation, S.Y.; investigation, S.Y.; methodology, S.Y. and S.Z.; project administration, D.S. and F.K.; supervision, D.S. and R.-A.E.; visualization, S.Y.; writing the original draft, S.Y.; writing—review and editing, S.Z., D.S., R.P. and S.Y. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by the German Federal Ministry of Education and Research (BMBF) grant numbers FKZ 03SF0621A.

Data Availability Statement

he data presented in this study are available on reasonable request to the corresponding author. Some data are not publicly available due to confidentiality.

Acknowledgments

The authors would like to thank all their colleagues engaged in SOC development at Jülich, including the institutes: IEK-1, IEK-2, IEK-9, and ZEA-3.

Conflicts of Interest

The authors have no conflict of interest to declare. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

Abbreviations

8YSZ8 mol% yttria-stabilized zirconia
AEL/ELECThe interface between the air electrode and the electrolyte
APSAtmospheric plasma spraying
FEL/ELECThe interface between the fuel electrode and the electrolyte
FESCFuel electrode-supported cell (FESC)
FZJForschungszentrum Jülich GmbH
GDCGadolinia-doped ceria
ITMIntermediate temperature metal
LCC12La0.97Mn0.4Co0.3Cu0.3O3-δ
LSCFLa0.58Sr0.4Co0.2Fe0.8O3-δ
MCFMnCo1.9Fe0.1O3-δ
MFCMass flow controller
MICMetal interconnect
NSNavier–Stokes
PVDPhysical vapor deposition
SOCSolid oxide cell
SPScreen printing
WPSWet powder spraying
Greek symbols
χ Volume fraction[1]
κ Permeability[m2]
ϕ ele Electronic potential[V]
ϕ io Ionic Potential[V]
σ i ele Electronic conductivity of phase i[S/m]
σ i io Ionic conductivity of phase i[S/m]
σ oxide scale Electronic conductivity of the oxide scale[S/m]
τ Tortuosity[1]
υ Diffusion volume[m 3 /mol]
ε Porosity[1]
Ξ Parabolic constant[m 2 /s]
Roman symbols
U Velocity[m/s]
cMolar density[mol/m 3 ]
dThickness of the oxide scale[m]
D i j B Binary diffusion coefficient[m 2 /s]
D i j K Knudsen diffusion coefficient[m 2 /s]
i ele Electronic current density[A/m 2 ]
i io Ionic current density[A/m 2 ]
i v Reaction current density[A/m 3 ]
MMolar mass[g/mol]
NMolar diffusion flux[mol/(m 2 s)]
PPressure[Pa]
S mass Source term of mass[g/(m 3 s)]
S mole Source term of mole[mol/(m 3 s)]
TTemperature[K]
YMolar fraction[1]
RGas constant[(J mol)/K]

Appendix A

Table A1. Other physical parameters used in the model.
Table A1. Other physical parameters used in the model.
ParameterValueUnitReference
Mean pore diameter of the fuel electrode ( d p , fel )641 nm  [59]
Mean pore diameter of the air electrode ( d p , ael )306 nm  [59]
Permeability in the fuel electrode ( κ fuel )1 × 10−13 m 2
Permeability in the air electrode ( κ air )1 × 10−13 m 2
Hydrogen diffusion volume ( υ H 2 )6.12 × 10−6 m 3 / mol  [60]
Water diffusion volume ( υ H 2 O )1.31 × 10−5 m 3 / mol  [60]
Nitrogen diffusion volume ( υ N 2 )1.79 × 10−5 m 3 / mol  [22]
Oxygen diffusion volume ( υ O 2 )1.66 × 10−5 m 3 / mol  [22]
Volume fraction of Ni in solid part of fuel electrode ( χ Ni )0.41 [25]
Volume fraction of YSZ in solid part of fuel electrode ( χ YSZ )1- χ Ni 1 [25]
Volume fraction of LSCF in solid part of air electrode ( χ LSCF )11 [57]
Porosity of the Ni mesh ( ε Ni mesh )0.81
Porosity of the fuel electrode ( ε fel )0.21 [57]
Porosity of the fuel electrode support layer ( ε felsup )0.4151 [25]
Porosity of the air electrode ( ε ael )0.451 [25]
Porosity of the contact layer ( ε ocl )0.451
Porosity of the WPS coating ( ε coating , WPS )0.451 [52]
Porosity of the APS coating ( ε coating , APS )0.031 [53]
Tortuosity of LSCF ( τ LSCF )2.21 [25]
Tortuosity of Ni ( τ Ni )4.641 [25]
Tortuosity of YSZ ( τ YSZ )2.181 [25]

References

  1. Hawkes, G.; O’Brien, J.; Stoots, C.; Hawkes, B. 3D CFD model of a multi-cell high-temperature electrolysis stack. Int. J. Hydrog. Energy 2009, 34, 4189–4197. [Google Scholar] [CrossRef] [Green Version]
  2. Nohl, M.; Mazumder, J.; Vinke, I.C.; Eichel, R.A.; de Haart, L.G.J. Modeling of Multi-Physics Phenomena for High-Temperature Co-Electrolysis. ECS Trans. 2021, 103, 797. [Google Scholar] [CrossRef]
  3. Zheng, J.; Xiao, L.; Wu, M.; Lang, S.; Zhang, Z.; Chen, M.; Yuan, J. Numerical Analysis of Thermal Stress for a Stack of Planar Solid Oxide Fuel Cells. Energies 2022, 15, 343. [Google Scholar] [CrossRef]
  4. Navasa, M.; Miao, X.Y.; Frandsen, H.L. A fully-homogenized multiphysics model for a reversible solid oxide cell stack. Int. J. Hydrog. Energy 2019, 44, 23330–23347. [Google Scholar] [CrossRef]
  5. Miao, X.Y.; Rizvandi, O.B.; Navasa, M.; Frandsen, H.L. Modelling of local mechanical failures in solid oxide cell stacks. Appl. Energy 2021, 293, 116901. [Google Scholar] [CrossRef]
  6. Mozdzierz, M.; Berent, K.; Kimijima, S.; Szmyd, J.S.; Brus, G. A Multiscale Approach to the Numerical Simulation of the Solid Oxide Fuel Cell. Catalysts 2019, 9, 253. [Google Scholar] [CrossRef] [Green Version]
  7. Wehrle, L.; Schmider, D.; Dailly, J.; Banerjee, A.; Deutschmann, O. Benchmarking solid oxide electrolysis cell-stacks for industrial Power-to-Methane systems via hierarchical multi-scale modelling. Appl. Energy 2022, 317, 119143. [Google Scholar] [CrossRef]
  8. Wang, Y.; Wehrle, L.; Banerjee, A.; Shi, Y.; Deutschmann, O. Analysis of a biogas-fed SOFC CHP system based on multi-scale hierarchical modeling. Renew. Energy 2021, 163, 78–87. [Google Scholar] [CrossRef]
  9. Zhang, S. Modeling and Simulation of Polymer Electrolyte Fuel Cells. Ph.D. Thesis, Forschungszentrum Jülich GmbH Zentralbibliothek, Verlag, Jülich, Germany, 2020. [Google Scholar]
  10. Zhang, S.; Hess, S.; Marschall, H.; Reimer, U.; Beale, S.B.; Lehnert, W. openFuelCell2: A New Computational Tool for Fuel Cells, Electrolyzers, and other Electrochemical Devices and Processes. Comput. Phys. Commun. 2023. In preparation. [Google Scholar]
  11. Ltd, O. OpenFOAM. Available online: https://www.openfoam.com (accessed on 18 February 2023).
  12. Zhang, S.; Reimer, U.; Beale, S.B.; Lehnert, W.; Stolten, D. Modeling Polymer Electrolyte Fuel Cells: A High Precision Analysis. Appl. Energy 2019, 233–234, 1094–1103. [Google Scholar] [CrossRef]
  13. Zhang, S.; Beale, S.B.; Reimer, U.; Andersson, M.; Lehnert, W. Polymer Electrolyte Fuel Cell Modeling—A Comparison of Two Models with Different Levels of Complexity. Int. J. Hydrogen Energy 2020, 45, 19761–19777. [Google Scholar] [CrossRef]
  14. Zhang, S.; Beale, S.B.; Shi, Y.; Janßen, H.; Reimer, U.; Lehnert, W. Development of an Open-Source Solver for Polymer Electrolyte Fuel Cells. Ecs Trans. 2020, 98, 317. [Google Scholar] [CrossRef]
  15. Zhang, S. Low-Temperature Polymer Electrolyte Fuel Cells. In Electrochemical Cell Calculations with OpenFOAM; Beale, S., Lehnert, W., Eds.; Lecture Notes in Energy; Springer International Publishing: Cham, Switzerland, 2022; pp. 59–85. [Google Scholar] [CrossRef]
  16. Steinberger-Wilckens, R.; De Haart, L.; Vinke, I.; Blum, L.; Cramer, A.; Remmel, J.; Blass, G.; Tietz, F.; Quadakkers, W. Recent Results of Stack Development at Forschungszentrum Jülich. In Proceedings of the Fuel Cell Technologies: State and Perspectives; Sammes, N., Smirnova, A., Vasylyev, O., Eds.; Springer: Dordrecht, The Netherlands, 2005; pp. 123–134. [Google Scholar]
  17. Geisler, H.I. Finite Element Method (FEM) Model and Performance Analysis of Solid Oxide Fuel Cells. Ph.D. Thesis, Karlsruher Institut für Technologie (KIT), Karlsruhe, Germany, 2019. [Google Scholar] [CrossRef]
  18. Blum, L.; Fang, Q.; Groß-Barsnick, S.M.; de Haart, L.B.; Malzbender, J.; Menzler, N.H.; Quadakkers, W.J. Long-term operation of solid oxide fuel cells and preliminary findings on accelerated testing. Int. J. Hydrog. Energy 2020, 45, 8955–8964. [Google Scholar] [CrossRef]
  19. Chan, S.H.; Khor, K.A.; Xia, Z.T. A complete polarization model of a solid oxide fuel cell and its sensitivity to the change of cell component thickness. J. Power Sources 2001, 93, 130–140. [Google Scholar] [CrossRef]
  20. Mogensen, M.; Frandsen, H.; Nielsen, J.; Li, W.; Jacobsen, T.; Graves, C. Current density—Overvoltage relations for solid oxide electrodes. In Proceedings of the Smart Energy Conversion and Storage Conference: IV Polish Forum, Krynica, Poland, 1–4 October 2013. [Google Scholar]
  21. Adler, S.B.; Bessler, W.G. Chapter 29—Elementary Kinetic Modeling of Solid Oxide Fuel Cell Electrode Reactions. In Handbook of Fuel Cells; John Wiley & Sons, Ltd.: Hoboken, NJ, USA, 2010; Volume 5, p. p. 441. [Google Scholar] [CrossRef]
  22. Navasa, M.; Graves, C.; Chatzichristodoulou, C.; Løye Skafte, T.; Sundén, B.; Lund Frandsen, H. A three dimensional multiphysics model of a solid oxide electrochemical cell: A tool for understanding degradation. Int. J. Hydrog. Energy 2018, 43, 11913–11931. [Google Scholar] [CrossRef]
  23. Chang, H.; Jaffé, G. Polarization in Electrolytic Solutions. Part I. Theory. J. Chem. Phys. 1952, 20, 1071–1077. [Google Scholar] [CrossRef]
  24. Mogensen, M.; Jacobsen, T. The Course of Oxygen Partial Pressure and Electric Potentials across an Oxide Electrolyte Cell, the Electrochemical Society. In Proceedings of the 213th ECS Meeting, Phoenix, AZ, USA, 18–23 May 2008; Volume 13, pp. 259–273. [Google Scholar]
  25. Joos, J. Microstructural Characterisation, Modelling and Simulation of Solid Oxide Fuel Cell Cathodes. Ph.D. Thesis, Karlsruher Institut für Technologie (KIT), Karlsruhe, Germany, 2017. [Google Scholar] [CrossRef]
  26. VDM. Crofer 22 APU. Available online: https://www.vdm-metals.com/fileadmin/user_upload/Downloads/Data_Sheets/Data_Sheet_VDM_Crofer_22_APU.pdf (accessed on 25 November 2022).
  27. Guo, M.; Ru, X.; Lin, Z.; Xiao, G.; Wang, J. Optimization Design of Rib Width and Performance Analysis of Solid Oxide Electrolysis Cell. Energies 2020, 13, 5468. [Google Scholar] [CrossRef]
  28. Kundu, P.K.; Cohen, I.M.; Dowling, D.R. (Eds.) Chapter 4—Conservation Laws; Academic Press: Boston, MA, USA, 2012; pp. 95–169. [Google Scholar] [CrossRef]
  29. Zeng, Z.; Grigg, R. A Criterion for Non-Darcy Flow in Porous Media. Transp. Porous Media 2006, 63, 57–69. [Google Scholar] [CrossRef]
  30. Brinkman, H.C. A calculation of the viscous force exerted by a flowing fluid on a dense swarm of particles. Flow Turbul. Combust. 1949, 1, 27. [Google Scholar] [CrossRef]
  31. Le Bars, M.; Worster, M.G. Interfacial conditions between a pure fluid and a porous medium: Implications for binary alloy solidification. J. Fluid Mech. 2006, 550, 149–173. [Google Scholar] [CrossRef] [Green Version]
  32. COMSOL Multiphysics v. 5.6 CFD User Guide, Chapter 7—Porous Media and Subsurface Flow. Available online: https://doc.comsol.com/5.6/doc/com.comsol.help.cfd/CFDModuleUsersGuide.pdf (accessed on 25 November 2022).
  33. Ni, M.; Leung, M.K.H.; Leung, D.Y.C. A modeling study on concentration overpotentials of a reversible solid oxide fuel cell. J. Power Sources 2006, 163, 460–466. [Google Scholar] [CrossRef]
  34. Ni, M. Computational fluid dynamics modeling of a solid oxide electrolyzer cell for hydrogen production. Int. J. Hydrog. Energy 2009, 34, 7795–7806. [Google Scholar] [CrossRef]
  35. Webb, S.W.; Pruess, K. The Use of Fick’s Law for Modeling Trace Gas Diffusion in Porous Media. Transp. Porous Media 2003, 51, 327–341. [Google Scholar] [CrossRef]
  36. Reiss, G.; Frandsen, H.L.; Persson, H.; Weiß, C.; Brandstätter, W. Numerical evaluation of oxide growth in metallic support microstructures of Solid Oxide Fuel Cells and its influence on mass transport. J. Power Sources 2015, 297, 388–399. [Google Scholar] [CrossRef]
  37. Fuller, E.N.; Schettler, P.D.; Giddings, J.C. New method for prediction of binary gas-phase diffusion coefficients. Ind. Eng. Chem. 1966, 58, 18–27. [Google Scholar] [CrossRef]
  38. Beale, S.B.; Andersson, M.; Boigues-Muñoz, C.; Frandsen, H.L.; Lin, Z.; McPhail, S.J.; Ni, M.; Sundén, B.; Weber, A.; Weber, A.Z. Continuum scale modelling and complementary experimentation of solid oxide cells. Prog. Energy Combust. Sci. 2021, 85, 100902. [Google Scholar] [CrossRef]
  39. Grondin, D.; Deseure, J.; Brisse, A.; Zahid, M.; Ozil, P. Simulation of a high temperature electrolyzer. J. Appl. Electrochem. 2010, 40, 933–941. [Google Scholar] [CrossRef]
  40. Serincan, M.F.; Pasaogullari, U.; Singh, P. Controlling reformation rate for a more uniform temperature distribution in an internal methane steam reforming solid oxide fuel cell. J. Power Sources 2020, 468, 228310. [Google Scholar] [CrossRef]
  41. Jian, P.; Jian, L.; Bing, H.; Xie, G. Oxidation kinetics and phase evolution of a Fe–16Cr alloy in simulated SOFC cathode atmosphere. J. Power Sources 2006, 158, 354–360. [Google Scholar] [CrossRef]
  42. Froitzheim, J.H. Ferritic Steel Interconnectors and Their Interactions with Ni Base Anodes in Solid Oxide Fuel Cells (SOFC). Ph.D. Thesis, Lehrstuhl für Werkstoffe der Energietechnik (FZ Jülich), Jülich, Germany, 2008. [Google Scholar]
  43. Asensio Jimenez, C. Effect of Composition, Microstructure and Component Thickness on the Oxidation Behaviour of Laves Phase Strengthened Interconnect Steel for Solid Oxide Fuel Cells (SOFC). Ph.D. Thesis, Lehrstuhl für Werkstoffe der Energietechnik (FZ Jülich), Jülich, Germany, 2013. [Google Scholar]
  44. Hope, G.A.; Ritchie, I.M. Mechanism of chromium oxidation. J. Chem. Soc. Faraday Trans. 1 1981, 77, 2621–2631. [Google Scholar] [CrossRef]
  45. Persson, Å.H. High Temperature Oxidation of Slurry Coated Interconnect Alloys. Ph.D. Thesis, Technical University of Denmark, Lyngby, Denmark, 2012. [Google Scholar]
  46. Unutulmazsoy, Y. Oxidation Kinetics of Metal Films and Diffusion in NiO for Data Storage. Ph.D. Thesis, University of Stuttgart, Stuttgart, Germany, 2016. [Google Scholar]
  47. Tsai, S.C.; Huntz, A.M.; Dolin, C. Diffusion of 18O in massive Cr2O3 and in Cr2O3 scales at 900 °C and its relation to the oxidation kinetics of chromia forming alloys. Oxid. Met. 1995, 43, 581–596. [Google Scholar] [CrossRef]
  48. Holt, A.; Kofstad, P. Electrical conductivity and defect structure of Cr2O3. I. High temperatures (>∼1000 C). Solid State Ionics 1994, 69, 127–136. [Google Scholar] [CrossRef]
  49. Tsai, S.C.; Huntz, A.M.; Dolin, C. Growth mechanism of Cr2O3 scales: Oxygen and chromium diffusion, oxidation kinetics and effect of yttrium. Mater. Sci. Eng. A 1996, 212, 6–13. [Google Scholar] [CrossRef]
  50. Glatz, W. P/M Processing and Properties of High Performance Interconnect Materials and Components for SOFC Applications. ECS Proc. Vol. 2005, 2005, 1773–1780. [Google Scholar] [CrossRef]
  51. Megel, S.; Girdauskaite, E.; Sauchuk, V.; Kusnezoff, M.; Michaelis, A. Area specific resistance of oxide scales grown on ferritic alloys for solid oxide fuel cell interconnects. J. Power Sources 2011, 196, 7136–7143. [Google Scholar] [CrossRef]
  52. Ruder, A.; Buchkremer, H.P.; Jansen, H.; Malléner, W.; Stöver, D. Wet powder spraying—a process for the production of coatings. Surf. Coatings Technol. 1992, 53, 71–74. [Google Scholar] [CrossRef]
  53. Grünwald, N.; Lhuissier, P.; Salvo, L.; Villanova, J.; Menzler, N.H.; Guillon, O.; Martin, C.L.; Vaßen, R. In situ investigation of atmospheric plasma-sprayed Mn-Co-Fe-O by synchrotron X-ray nano-tomography. J. Mater. Sci. 2020, 55, 12725–12736. [Google Scholar] [CrossRef]
  54. Patankar, S.V.; Spalding, D.B. A calculation procedure for heat, mass and momentum transfer in three-dimensional parabolic flows. In Numerical Prediction of Flow, Heat Transfer, Turbulence and Combustion; Elsevier: Amsterdam, The Netherlands, 1983; pp. 54–73. [Google Scholar]
  55. Schäfer, D.; Queda, L.; Nischwitz, V.; Fang, Q.; Blum, L. Origin of Steam Contaminants and Degradation of Solid-Oxide Electrolysis Stacks. Processes 2022, 10, 598. [Google Scholar] [CrossRef]
  56. Fang, Q.; Blum, L.; Stolten, D. Electrochemical Performance and Degradation Analysis of an SOFC Short Stack Following Operation of More than 100,000 Hours. J. Electrochem. Soc. 2019, 166, F1320–F1325. [Google Scholar] [CrossRef]
  57. Menzler, N.H.; Sebold, D.; Sohn, Y.J.; Zischke, S. Post-test characterization of a solid oxide fuel cell after more than 10 years of stack testing. J. Power Sources 2020, 478, 228770. [Google Scholar] [CrossRef]
  58. Menzler, N.H.; Batfalsky, P.; Groß, S.; Shemet, V.; Tietz, F. Post-Test Characterization of an SOFC Short-Stack after 17,000 Hours of Steady Operation. ECS Trans. 2011, 35, 195–206. [Google Scholar] [CrossRef]
  59. Geisler, H.; Kromp, A.; Weber, A.; Ivers-Tiffée, E. Stationary FEM Model for Performance Evaluation of Planar Solid Oxide Fuel Cells Connected by Metal Interconnectors. J. Electrochem. Soc. 2014, 161, F778–F788. [Google Scholar] [CrossRef]
  60. Kong, W.; Zhu, H.; Fei, Z.; Lin, Z. A modified dusty gas model in the form of a Fick’s model for the prediction of multicomponent mass transport in a solid oxide fuel cell anode. J. Power Sources 2012, 206, 171–178. [Google Scholar] [CrossRef]
Figure 1. Simplification of the geometry of a standard F10 stack [16]. The left figure is a representive geometry of the active volume. The right figure is the cross-section geometry used in the simulations. The barrier layer and protective coating are not considered in the model. Not to scale.
Figure 1. Simplification of the geometry of a standard F10 stack [16]. The left figure is a representive geometry of the active volume. The right figure is the cross-section geometry used in the simulations. The barrier layer and protective coating are not considered in the model. Not to scale.
Energies 16 03827 g001
Figure 2. Schematic diagram of the ultra-thin film resistance model. To more clearly illustrate the model, only the contact layer and MIC are shown in the figure.
Figure 2. Schematic diagram of the ultra-thin film resistance model. To more clearly illustrate the model, only the contact layer and MIC are shown in the figure.
Energies 16 03827 g002
Figure 3. The mesh of the model.
Figure 3. The mesh of the model.
Energies 16 03827 g003
Figure 4. A schematic diagram showing how each of the physics is coupled with each other.
Figure 4. A schematic diagram showing how each of the physics is coupled with each other.
Energies 16 03827 g004
Figure 5. Validations of I–V curves for (a) F1002-97, 12% humidified fuel, SOFC, (b) F1004-67, 20% humidified fuel, SOFC and (c) F1004-115, 50% humidified fuel, SOFC and SOEC.
Figure 5. Validations of I–V curves for (a) F1002-97, 12% humidified fuel, SOFC, (b) F1004-67, 20% humidified fuel, SOFC and (c) F1004-115, 50% humidified fuel, SOFC and SOEC.
Energies 16 03827 g005
Figure 6. Validation of MIC oxidation model. The experimental data are taken from TGA results of Persson [45]. Crofer 22 APU was used as the MIC. One stack was operated without coatings, while another stack was coated with 15 µm LSM, whose porosity was measured as 0.3 after 4000 h of operation [45].
Figure 6. Validation of MIC oxidation model. The experimental data are taken from TGA results of Persson [45]. Crofer 22 APU was used as the MIC. One stack was operated without coatings, while another stack was coated with 15 µm LSM, whose porosity was measured as 0.3 after 4000 h of operation [45].
Energies 16 03827 g006
Figure 7. Validations of MIC oxidation model. The experimental data (21 data points for F1002-97 and 29 data points for F1002-95) are obtained from a scanning electron microscope (SEM) image of FZJ SOC stacks after they were shut down [57,58].
Figure 7. Validations of MIC oxidation model. The experimental data (21 data points for F1002-97 and 29 data points for F1002-95) are obtained from a scanning electron microscope (SEM) image of FZJ SOC stacks after they were shut down [57,58].
Energies 16 03827 g007
Figure 8. Spatial distribution of mole fractions in the x-y plane of the fuel electrode and the air electrode. The simulation is the electrolysis operation of F1004-115, where fuel flows toward the x-coordinate direction and air flows against the x-coordinate. The geometry is scaled by the transformation vector (x, y, z) = (1, 20, 1).
Figure 8. Spatial distribution of mole fractions in the x-y plane of the fuel electrode and the air electrode. The simulation is the electrolysis operation of F1004-115, where fuel flows toward the x-coordinate direction and air flows against the x-coordinate. The geometry is scaled by the transformation vector (x, y, z) = (1, 20, 1).
Energies 16 03827 g008
Figure 9. Spatial distribution of overpotential on the outer surface of the fuel electrode and the air electrode in the x-z plane (see geometry in Figure 1). The simulation is the electrolysis operation of F1004-115, where fuel flows toward the x-coordinate direction and air flows against the x-coordinate. The geometry is scaled by the transformation vector (x, y, z) = (0.05, 1, 60).
Figure 9. Spatial distribution of overpotential on the outer surface of the fuel electrode and the air electrode in the x-z plane (see geometry in Figure 1). The simulation is the electrolysis operation of F1004-115, where fuel flows toward the x-coordinate direction and air flows against the x-coordinate. The geometry is scaled by the transformation vector (x, y, z) = (0.05, 1, 60).
Energies 16 03827 g009
Figure 10. Simulation of the stack voltage evolution of F1002-97 with time. Only MIC oxidation is considered to be the degradation. The inserted figures are the distribution of the electronic potential across the contact layer and the MIC.
Figure 10. Simulation of the stack voltage evolution of F1002-97 with time. Only MIC oxidation is considered to be the degradation. The inserted figures are the distribution of the electronic potential across the contact layer and the MIC.
Energies 16 03827 g010
Figure 11. (a) Voltage evolution from experiments and simulation. The red texts refer to the degradation rates: (b) Ratio of the degradation rate due to MIC oxidation to the total degradation rate. Comparison of the experimental data and simulation results of F1002-97: (a) shows the voltage evolution only due to MIC oxidation and (b) exhibits the share of voltage degradation rate caused by MIC oxidation. A sudden increase is observed in (b) because the total degradation rate decreases from 7.3 mV/kh to 1.6 mV/kh at 40 kh.
Figure 11. (a) Voltage evolution from experiments and simulation. The red texts refer to the degradation rates: (b) Ratio of the degradation rate due to MIC oxidation to the total degradation rate. Comparison of the experimental data and simulation results of F1002-97: (a) shows the voltage evolution only due to MIC oxidation and (b) exhibits the share of voltage degradation rate caused by MIC oxidation. A sudden increase is observed in (b) because the total degradation rate decreases from 7.3 mV/kh to 1.6 mV/kh at 40 kh.
Energies 16 03827 g011
Table 1. Summary of main models used in this work.
Table 1. Summary of main models used in this work.
Physical PhenomenonDomain or SurfaceMain EquationVariables to Be Solved
Electrochemical reactionElectrodesChang–Jaffe kineticElectronic and ionic potential
Momentum transferChannelsNavier–Stokes equationsPressure and velocity
Porous mediumDarcy–Forchheimer equationPressure and velocity
Mass transportChannels and porous mediumFick’s lawMolar fractions of species
MIC oxidationRib surfaceParabolic lawElectronic potential difference due to the oxide scale
Table 2. Conductivity of material.
Table 2. Conductivity of material.
Component σ ele (S/m) [22,26,27] σ io (S/m) [22]
Ni3.27 × 106 −1065.3T
8YSZ 6.25 × 104 exp(−10,300/T)
LSCF22,591–1.6 × 106 exp(−6024/T)5.5 × 109 exp(−9050/R/T)
LCC12 122,591–1.6 × 106 exp(−6024/T)5.5 × 109 exp(−9050/R/T)
Crofer 22 APU9.09 × 105
1Assume LCC12 has the same electronic and ionic conductivity as LSCF.
Table 3. Experimental details of SOC stacks.
Table 3. Experimental details of SOC stacks.
F1002-97F1002-95F1004-67F1004-115
InstrumentMFCBronkhorst, F-201CVBronkhorst, F-201CVBronkhorst, F-201CVBrooks, 0254
Stack designContact layerLCC12LCC12LSCFLSCF
MICITMCrofer 22 APUCrofer 22 APUCrofer 22 APU
Protective coatingWPS MnOxWPS MnOxAPS MCFAPS MCF
Operation conditionsTemperature720 C710 C730 C
Fuel mass flow1.18 × 10−7 kg/s1.15 × 10−7 kg/s1.15 × 10−7 kg/s
Molar ratio in fuelH2/H2O = 79/21H2/H2O = 80/20H2/H2O = 80/20
Air mass flow1.9 × 10−6 kg/s1.18 × 10−6 kg/s1.18 × 10−6 kg/s
Molar ratio in airO2/N2 = 21/79O2/N2 = 21/79O2/N2 = 21/79
Current density0.5 A/cm 2 0.5 A/cm 2 0.5 A/cm 2
Time100 kh17 kh25 kh
Conditions of I-V
curves characterization
Temperature700∼800 C 700∼800 C700∼800 C
Fuel mass flow1.52 × 10−7 kg/s 1.85 × 10−7 kg/s5.53 × 10−7 kg/s
Molar ratio in fuelH2/H2O = 88/12 H2/H2O = 80/20H2/H2O = 50/50
Air mass flow2.38 × 10−6 kg/s 1.79 × 10−6 kg/s2.83 × 10−6 kg/s
Molar ratio in airO2/N2 = 21/79 O2/N2 = 21/79O2/N2 = 21/79
Current density0∼0.8 A/cm 2 0∼0.8 A/cm 2 −0.5∼0.5 A/cm 2
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Share and Cite

MDPI and ACS Style

Yu, S.; Zhang, S.; Schäfer, D.; Peters, R.; Kunz, F.; Eichel, R.-A. Numerical Modeling and Simulation of the Solid Oxide Cell Stacks and Metal Interconnect Oxidation with OpenFOAM. Energies 2023, 16, 3827. https://doi.org/10.3390/en16093827

AMA Style

Yu S, Zhang S, Schäfer D, Peters R, Kunz F, Eichel R-A. Numerical Modeling and Simulation of the Solid Oxide Cell Stacks and Metal Interconnect Oxidation with OpenFOAM. Energies. 2023; 16(9):3827. https://doi.org/10.3390/en16093827

Chicago/Turabian Style

Yu, Shangzhe, Shidong Zhang, Dominik Schäfer, Roland Peters, Felix Kunz, and Rüdiger-A. Eichel. 2023. "Numerical Modeling and Simulation of the Solid Oxide Cell Stacks and Metal Interconnect Oxidation with OpenFOAM" Energies 16, no. 9: 3827. https://doi.org/10.3390/en16093827

APA Style

Yu, S., Zhang, S., Schäfer, D., Peters, R., Kunz, F., & Eichel, R. -A. (2023). Numerical Modeling and Simulation of the Solid Oxide Cell Stacks and Metal Interconnect Oxidation with OpenFOAM. Energies, 16(9), 3827. https://doi.org/10.3390/en16093827

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop