The new design of the MS-SOEC is shown in Figure 1A. A 5 μm GDC layer is used as the main electrolyte while a 1 μm 10Sc1CeSZ layer is deposited on both sides of the GDC layer. As a comparison, the cell structure of the MS-SOFC developed by Ceres Power (Leah et al., 2005) is displayed in Figure 1B. MS-SOEC has a high reliability and can endure frequent quick starts. The ferritic substrate normally has a thickness of 100–300 μm. Small holes with a diameter of 10–30 μm are drilled by lasers to form a porous region. The cathode layer is deposited on the porous substrate, which is composed of Ni/GDC with a thickness of 15 μm. The main GDC electrolyte layer is fabricated by screen printing or tape casting. For the case of Figure 1B, the thickness is set to 5.8 μm. To prevent electron leakage, a 1 μm thick YSZ layer is deposited on the surface of GDC. Finally, a 0.2 μm thin GDC layer is deposited on YSZ layer as a barrier layer to prevent the direct solid reaction of YSZ with the anode material (LSCF) and the formation of low-conductivity impure phase. The anode (LSCF/GDC) and cathode (Ni/GDC) are the same with that of the MS-SOFC of Ceres Power.

During the working process of the MS-SOEC, high-temperature steam-hydrogen mixture is supplied to the cathode and diffused to the triple phase boundary (TPB). Therein, steam is electrolyzed as Eq. 1 shows. The generated oxygen ions are transported to the anode through the dense electrolyte layer. In the anode, Reaction (2) occurs and oxygen is generated. H 2 O + 2 e − → O 2 − + H 2 (1) O 2 − − 2 e − → O 2 (2)

An electrochemical model is setup to describe the performance of the MS-SOEC. The overall input voltage is determined by E = E 0 + η a c t , a + η a c t , c + η c o n , a + η c o n , c + η o h m (3) where E is the input voltage, E₀ is the thermodynamically reversible voltage that the electrochemical reactions need, η a c t , a and η a c t , c are the activation overpotentials of the anode and the cathode, η c o n , a and η c o n , c are the corresponding concentration overpotentials, and η o h m is the ohmic loss.

The reversible voltage is predicted by the following Nernst equation according to the temperature and pressure of the species. E 0 = Δ G H 2 0 2 F + R m T 2 F ln ( P H 2 O i n P H 2 i n ( P O 2 i n ) 0.5 ) (4) where ∆ G H 2 0 ( J / m o l ) is the Gibbs free energy of steam electrolysis reaction, T (K) is the operating temperature, P H 2 i n and P O 2 i n are the initial partial pressures of the hydrogen and the oxygen, respectively.

The activation overpotentials are the voltage losses controlled by the kinetics of the electrochemical reactions. It is expressed by the Butler-Volmer equation as a function of the current density.

For the anode, J I = j a 0 { exp ( α n e F η a c t , a R m T ) − exp [ − ( 1 − α ) n e F η a c t , a R m T ] } (5)

For the cathode, J I = j c 0 { exp ( α n e F η a c t , c R m T ) − exp [ − ( 1 − α ) n e F η a c t , c R m T ] } (6) where J _I (A/m²) is the operating current density, j a 0 and j c 0 are the exchange current densities of the anode and the cathode, α is the charge transfer coefficient and normally set to 0.5, n _e is the electron number of the reaction, η a c t , a and η a c t , c are the activation overpotentials in a unit of V. The exchange current density of the anode is denoted by j a 0 = P O 2 R   (0.25) K a R m T 2 F exp ( − E a , a c t R m T ) (7) where P O 2 R is the interface partial pressure of oxygen, K _a is the pre-exponential factor of the anode, E _a,act is the activation energy. The exchange current density of the cathode is expressed as j c 0 = P H 2 O L P H 2 L K c R m T 2 F exp ( − E c , a c t R m T ) (8) where P H 2 O L and P H 2 L are the interface partial pressures of steam and hydrogen, K _c is the pre-exponential factor of the cathode, E _c,act is the activation energy.

The concentration overpotential occurs due to the mass transport inside the porous electrodes and cannot be neglected under high operating currents. The concentration losses of the anode and the cathode are denoted by Eqs 9, 10, respectively (Ni et al., 2007). η c o n c , a = R m T 2 F ln [ ( P O 2 R P O 2 o u t ) 1 / 2 ] (9) η c o n c , c = R m T 2 F ln [ P H 2 L P H 2 O i n P H 2 i n P H 2 O L ] (10) where η c o n c , a and η c o n c , c are the concentration overpotentials, P O 2 R is the partial oxygen pressure at the interface between the anode and the electrolyte layer, P O 2 o u t is the partial pressure of oxygen at the interface between the anode and the interconnecting plate, P H 2 L   a n d   P H 2 O L are the partial pressures at the interface between the cathode and the electrolyte layer, P H 2 i n   a n d   P H 2 O i n are the initial partial pressures of the hydrogen and the steam at the inlet of the cathode.

The dust gas model is used to describe the mass transport processes inside the electrodes. There are two diffusion mechanisms inside the porous electrode: molecular diffusion and Knudsen diffusion. When the hole diameter is much greater than the mean free path of the molecule, the actions among molecules dominate the diffusion process. Otherwise, the collision between the molecule and the hole wall is the main part (Knudsen mechanism). For the mass transport of the electrodes, both mechanisms are important. The effective diffusion coefficient is expressed by (Andreassi et al., 2007; Shi et al., 2022) 1 D i , j e f f = ξ ε [ 1 D i , k + 1 D i , j ] (11) where D i , j e f f is the effective diffusion coefficient of the gas, ε is the electrode porosity, ξ is the electrode tortuosity, D _i,k is the Knudsen diffusion coefficient, D _i,j is the molecule binary diffusion coefficient of the species i.

The Knudsen diffusion coefficient is determined by D i , k = r ¯ 2 3 8 R m T π M i (12) where r ¯ is the average diameter of the holes in the metal substrate, M _i (kg/kmol) is the relative molecular weight.

The binary diffusion coefficient is denoted by (Todd and Young, 2002) D i , j = 0.00143 T 1.75 P M i , j 1 / 2 ( V i 1 / 3 + V j 1 / 3 ) 2 (13) where M _i,j is the relative molecular weight of the mixture, V _i is the special diffusion volume of the species i. M i , j = 2 / ( 1 M i + 1 M j ) (14) where M _i and M _j are the relative molecular weights of the species i and j.

The mass transport of the steam and the hydrogen in the porous cathode is important for the kinetics of the electrochemical reaction. It is assumed that the mean diffusion path is estimated by the thickness of the small holes of the substrate and the hole pitch. The partial pressures of the hydrogen and oxygen at the TPB are determined by (Kim, 1991; Virkar et al., 2000; Chan et al., 2001) P H 2 L = P H 2 i n + R m T ( τ h o l e p i t c h + τ s u b s t r a t e ) 2 F D e f f , c a t h o d e J I (15) P H 2 O L = P H 2 O i n − R m T ( τ h o l e p i t c h + τ s u b s t r a t e ) 2 F D e f f , c a t h o d e J I (16) where τ _holepitch is the hole spacing of the substrate, τ _subatrate is the thickness of the metal substrate, D _{eff, cathode} is the effective diffusion coefficient of the gases in the cathode.

Similarly, the partial pressure of the oxygen in the TPB of the anode is denoted by (Laurencin et al., 2011) P O 2 R = 1 + ( P O 2 o u t − 1 ) exp ⁡ ⁡ ( − R m T J I τ a n o d e 4 F D e f f , a n o d e P ) (17) where τ_anode is the thickness of the anode, D _{eff, anode} is the thickness of the metal substrate, P is the operating pressure of the anode.

The ohmic overpotential considers the mass transfer of oxygen ions and electrons inside the electrolyte and the electrodes and is expressed by η o h m = ( τ a n o d e σ a n o d e + τ c a t h o d e σ c a t h o d e + τ e l i σ e l i ) J I + η c o n t (18) where η _ohm is ohmic overpotential, τ _anode , τ _cathode , and τ _ei are the thicknesses of the anode, the cathode and the electrolyte, respectively, σ _anode , σ _cathode , and σ _ei are the conductivities of the anode, cathode and electrolyte, respectively, R _cont is the overall contact resistance between the adjacent layers. The ionic conductivities of GDC, YSZ, and 10Sc1CeSZ are determined according to the equations in (Sahibzada et al., 1997; Leah et al., 2005).

Oxygen ion flow is the main part of the mass transport in the multilayer electrolyte. However, for the MIEC material such as GDC, the electronic transport should also be considered. Equivalent circuit model can be used to model the equilibrium potential inside the electrolyte and hence the partial pressure of the oxygen ions can be estimated (Virkar, 1991). Figure 2 shows the equivalent circuit for the designed multilayer electrolyte. The left side is the cathode and the right side is the anode. From the cathode to the anode, the layers of the multilayer electrolyte are labelled with 1 through 3. Point L is the contact site between the first layer of the electrolyte and the cathode. Point C is the contact site between the second layer and the first layer of the electrolyte. Point D is the contact site between the third layer and the second layer of the electrolyte. Point R is the contact site between the anode and the third layer of the electrolyte.

According to Nernst-Plank equation, the oxygen chemical potential can be expressed as a function of the oxygen partial pressures on the two sides of a single electrolyte layer. E i = R m T 4 F ln [ p O 2 ( L ) p O 2 ( R ) ] (19) where E _i is the electromotive force determined by Nernst-Planck equation, p O 2 ( L ) is the partial pressure of oxygen at the interface between the hydrogen electrode and the electrolyte layer, p O 2 ( R ) is the partial pressure of oxygen at the interface between the oxygen electrode and the electrolyte layer.

In fact, there are two flow paths in parallel inside the electrolyte: the ionic transport path and the electrical path. The ionic and electrical resistances of the kth layer can be determined by R i ( k ) = τ ( k ) σ i ( k ) (20) R e ( k ) = τ ( k ) σ e ( k ) (21) where R _i (k) and R _e (k) are the ionic and electronic resistances of the electrolyte at the kth layer in the equivalent circuit, respectively, τ(k) is the thickness of the kth layer, σ _i (k) and σ _e (k) are the ionic and electronic conductivities of the kth layer, respectively. Herein, k equals to 1, 2, or 3 counted from the left to right according to the sequence of the multilayer electrolyte shown in Figure 2.

According to the Kirchhoff law, the following equations are obtained. J i ( 1 ) = J i ( 2 ) = J i ( 3 ) = J I   (22) J e ( 1 ) = J e ( 2 ) = J e ( 3 ) = J e   (23) E = ϕ R − ϕ L   (24) ϕ C − ϕ L R e ( 1 ) = ϕ D − ϕ C R e ( 2 ) = ϕ R − ϕ D R e ( 3 )   (25) J = J I + J e   (26) ϕ C − ϕ L − E i ( 1 ) R i ( 1 ) = ϕ D − ϕ C − E i ( 2 ) R i ( 2 ) = ϕ R − ϕ D − E i ( 3 ) R i ( 3 )   (27) where J _i is the current density of the oxygen ion, J _e is the current density of the electron, Φ _R is the potential at the position between the anode and the electrolyte layer, Φ _L is the potential at the position between the cathode and the electrolyte layer, E is the voltage between the positions of L and R, Φ _C and Φ _D are the potentials at the positions of C and D, just shown as Figure 2, J is the overall current density, E _i (k) is the oxygen chemical potential at the kth electrolyte layer.

When the circuit is in equilibrium, the electrical current is denoted by J e = E ∑ k = 1 3 R e ( k ) (28)

If the potential at the left side is assumed zero, the potentials at points C and D are denoted by ∅ C = ∅ L + R e ( 1 ) ∑ k = 1 3 R e ( k ) E (29) ∅ D = ∅ R − R e ( 3 ) ∑ k = 1 3 R e ( k ) E (30)

The oxygen chemical potentials of the three layers of the electrolyte are expressed by E i ( 1 ) = E ( R e ( 1 ) + R i ( 1 ) ) ∑ k = 1 3 R e ( k ) − J R i ( 1 )   (31) E i ( 2 ) = E ( R e ( 2 ) + R i ( 2 ) ) ∑ k = 1 3 R e ( k ) − J R i ( 2 )   (32) E i ( 3 ) = E ( R e ( 3 ) + R i ( 3 ) ) ∑ k = 1 3 R e ( k ) − J R i ( 3 )   (33)

Finally, the oxygen partial pressures at the interfaces between the adjacent layers inside the electrolyte are determined according to Eq. 19.

In this study, the mathematical model of the MS-SOEC was implemented in the MATLAB R2018a software. The corresponding key parameters and boundary conditions are listed in Tables 1 and 2, respectively. The data for the oxygen ionic conductivity are given in Table 1. The data for the electrical conductivity of GDC are the same as Ref. (Leah et al., 2005). The electrical conductivity of YSZ is set to one thousandth of the corresponding oxygen ionic conductivity of this material. Similar assumption is used for 10Sc1CeSZ. The established mathematical model is validated based on the MS-SOFC shown in Figure 1B. The results of the polarization curve are compared with that of Ceres Power at 873 K (Leah et al., 2005), and displayed in Figure 3A. The porosity is set to 40% and the tortuosity is 4.5. The tendency of the simulated profile is consistent with that in (Leah et al., 2005). The relative deviations are shown in Figure 3B. When the operating current density is less than 0.65 A/cm², the relative deviation is less than 2%. The maximum relative deviation is 3.05% in the region with a large current density. The precision of the built model is acceptable and can be used for the performance estimation of the designed MS-SOEC. Since the data of SOEC mode similar to such a design cannot be found, the results of SOFC mode are used to validate the established mathematical model. For the SOEC mode, the polarization losses are assumed as the same with that of the SOFC mode. It is assumed that the default porosity is 40% and the tortuosity is 4.5 for the following analysis.

The effects of porosity and tortuosity on the performance of MS-SOEC are evaluated at first. These parameters will affect the mass transports in the cathode and anode. An inappropriate configuration will increase the concentration losses. The operating temperature is set to 600°C and a mixture of steam and hydrogen is supplied to the cathode with a steam mole concentration of 0.8. Air is input to the anode. First, Eq. 3 through (18) are used to determine the performance of the SOEC. The effects of porosity and tortuosity are evaluated based on Eqs. (11) through (14). The results are displayed in Figures 4A,B, respectively, when the porosities are set to 20 and 60%. The concentration and ohmic losses increase almost linearly with the increase of the current density. The activation loss ascends rapidly in the small current density region and increase slowly when the current density is large. The variation of the porosity almost has no effect on the ohmic losses. However, as the porosity rises, the concentration loss declines evidently. When the current density is 0.7 A/cm², the concentration loss is 43.72 mV with a porosity of 20% while it is 16.85 mV when the porosity is 60%, which is decreased by 61.46%.

The concentration losses of the electrodes are compared as the porosity increases from 0.2 to 0.6. The results are shown in Figure 5. Both the concentration losses decrease as the porosity increases, especially in the regions with a large current density. When the current density is 0.7A/cm², the concentration losses of the cathode and anode are decreased by 61.41 and 66.36%, respectively, as the porosity increases from 0.2 to 0.6. Meanwhile, the concentration overpotential of the anode is obviously less than that of the cathode. The overpotential of the cathode is 81–107 times greater than the anode and this proportion increases with the decrement of the current density or the increment of the porosity. This is because of the existence of the metal substrate on the cathode. As the current density increases, the mass flow rates of the species need to be increased as well. This requires a large mass transport capacity and a large porosity is beneficial for the improvement of the diffusion process. This is because the probability of the collision of the molecules against the hole walls is decreased for a larger porosity. The effective diffusion coefficients of the electrodes are shown in Figure 5C. The effective diffusion coefficients increase linearly as the current density rises. When the porosity is increased from 0.2 to 0.6, the effective diffusion coefficient of the cathode is increased from 0.712 × 10⁻⁶ m²/s to 2.137 × 10⁻⁶ m²/s while it is increased from 2.460 × 10⁻⁶ m²/s to 7.381 × 10⁻⁶ m²/s in the anode. The effective diffusion coefficient of the anode is about 3.45 times of the cathode, indicating that the metal substrate in the cathode has an obvious influence on the mass transport.

The effects of the tortuosity are analyzed on the same working conditions with a porosity of 40%. Figure 6 shows the polarization curves and the different overpotentials when the tortuosity is 2.5 and 6.5, respectively. For a fixed tortuosity, the variations of the overpotentials versus the current density are similar with that in Figure 4. However, the proportion of the concentration loss increases apparently as the tortuosity rises. When the current density is 0.7A/cm², a concentration overpotential of 14.29 mV is obtained for a tortuosity of 2.5 while it increases to 33.14 mV when the tortuosity is 6.5, which is increased by 131.9%.

A comparison of the concentration losses of the cathode and the anode is shown in Figure 7 when the tortuosity increases from 2.5 to 6.5. The concentration overpotentials of the anode and cathode are increased gradually with the increment of the tortuosity. As the current density rises, this augmentation tendency is more significant. The reason is attributed to a poorer mass transport due to a larger tortuosity. As the current density rises, the amounts of the reactants need to be increased, the mass flow rates must be increased accordingly, leading to a larger flow resistance. When the current density is 0.7A/cm², the concentration overpotentials of the anode and the cathode are increased by 158.4 and 131.7%, respectively, as the tortuosity rises from 2.5 to 6.5. Figure 7C shows the effective diffusion coefficients of the electrodes as a function of the tortuosity. As the tortuosity increases, these coefficients decrease. The effects of the tortuosity on the cathode are more significant than that on the anode, especially when the tortuosity is small.

During the working process of an SOEC, the operating temperature, pressure, and the steam concentration in the cathode are three important parameters. The influences of these parameters on the steam electrolysis performance are analyzed in this section.

The electrolyte GDC is suitable for intermediate-temperature conditions (550–650°C). It has been successfully used in MS-SOFC systems. However, the electronic conductivity is serious under a reducing atmosphere due to the reduction of Ce⁴⁺ to Ce³⁺. A thin and dense barrier layer is employed to prevent the internal electronic short circuit. For such a design with multilayer electrolyte, a very high oxygen partial pressure may occur at the interfaces and inside the electrolyte (Zhang et al., 2020), which exceeds the mechanical limit of the material and causes delamination and degradation under the operation conditions with a high current density. The degradation of the electrolyte in an electrochemical device is important and must be considered seriously (Liao, 2020). Current investigations only considered the cases with two layers of electrolyte, for example, using YSZ as the main electrolyte and GDC as the barrier layer. It is necessary to give an evaluation of the oxygen partial pressure for the two designs in Figure 1 with three layers of electrolyte. In this study, the oxygen partial pressure issue is explored based on the equivalent circuit model.

When the operating temperature is 600°C, the oxygen partial pressures at the interfaces between two adjacent layers inside the designed multilayer electrolyte are determined along the thickness direction. Herein, Case 1 represents the new design as 10Sc1CeSZ|GDC|10Sc1CeSZ displayed in Figure 1A and Case 2 is the conventional design of Ceres Power as GDC|YSZ|GDC displayed in Figure 1B. The multilayer electrolytes for the Cases 1 and 2 are composed of three layers shown as Figures 1, 2. With regard to the oxygen ion transport path, each layer has an ohmic potential decrement and an oxygen chemical potential increment E_i(k). When the current density is 0.17 and 0.7 A/cm², the results of these two different potentials are compared in Figure 11. The number of x-axes labels the layer of the electrolyte. Layer one denotes the left layer in the proximity of the cathode and Layer 3 is the one close to the anode. When the current density is 0.17A/cm², the results of the ionic ohmic potential are compared between the Case 1 and Case 2 in Figure 11A. Figure 11B gives the results of the oxygen chemical potential. Meanwhile, Figures 11C,D are the results with a current of 0.7A/cm².

When 10Sc1CeSZ is used as Case 1 shows, the ionic conductivity of 10Sc1CeSZ approximates to the value of GDC according to the data in Table 1. Meanwhile, the oxygen ohmic potential is also affected by the thickness of each layer. As a result, the ohmic potential of the GDC layer is evidently greater than that of 10Sc1CeSZ. When YSZ is used as the barrier layer in Case 2, because the GDC Layer one is much thicker than the other two layers, the ohmic potential of GDC is close to the YSZ Layer 2. On the other hand, the ohmic potential of each layer is much smaller than the associate chemical potential by almost 2 orders of magnitude. The oxygen chemical potentials of the 10Sc1CeSZ and YSZ are significantly greater than that of GDC layer. This is mainly because of the constraints of the electronic path. The difference of the electronic conductivity of the electrolytes has a great influence on the oxygen chemical potential distribution. With regard to the transport path of the electrons in Figure 2, the potential is mainly imposed on the YSZ layer and the 10Sc1CeSZ layers because GDC can be treated as a conductor compared with these two materials according to the data of electrical conductivity. Then, considering the oxygen ion transport path in Figure 2, this potential is the sum of the oxygen ohmic potential and the value of E_i(k). Because the E_i(k) is evidently greater than the ohmic potential shown as Figures 11A,B, the distribution of E_i(k) is almost the same as the potential in the electron transport path. Therefore, the oxygen chemical potential increment of GDC in these two cases are almost zero. Generally, as the current density increases, the ohmic potential of each layer increases apparently. The oxygen chemical potential is basically imposed on the YSZ layer for the Case 2. This is not helpful for the suppression of delamination. As a contrast, the overall chemical potential is evenly distributed between the two 10Sc1CeSZ layers for Case 1. If the oxygen chemical potential is too large, oxygen molecules may precipitate and enter into the nano voids that exist in the electrolyte. The gas pressure in the nano voids will continue to increase and exceed the mechanical limit of the electrolyte, causing the delamination. The results indicate that the new design using 10Sc1CeSZ as the barrier layer can alleviate the oxygen partial pressure issue of GDC significantly and the durability can be improved accordingly.

Figure 12 shows the results of the partial oxygen pressure and the corresponding potential at the interfaces of the multilayer electrolyte. The x-axis is the thickness direction of the electrolyte. The interface between the cathode and the multilayer electrolyte is denoted by 0 (corresponding to the point L in Figure 2) and the interface between the anode and the electrolyte is at 7 μm (corresponding to the point R in Figure 2). Points C and D are the two interfaces inside the multilayer electrolyte of Case 1 labelled as Figure 2. For Case 2, these two points are replaced with Points C′ and D’. The results of the partial pressure are shown in Figures 12A,B when the current densities are 0.17 and 0.7A/cm². For Case 1 where 10Sc1CeSZ is deposited on both sides of the GDC layer, the values of the oxygen partial pressure on the two interfaces (corresponding to the points C and D in Figure 2) are 9.43 × 10⁻¹¹ and 1.11 × 10⁻¹⁰ bar, respectively. In contrast, when the barrier layer YSZ is used as Case 2 shown, the oxygen partial pressures at the two interfaces of GDC layers (corresponding to the points L and D′ in Figure 2) are 9.15 × 10⁻²¹ and 1.004 bar. If the oxygen partial pressure of GDC is lower than the order of 10⁻¹⁶ bar, decomposition of GDC may occur according to the thermodynamic equilibrium. It is evident that the new design can improve the oxygen partial pressure of the GDC layer close to the cathode and a too small partial pressure is avoided, just as Case 2 shows. Meanwhile, for Case 1, the oxygen partial pressure of the GDC layer adjacent to the 10Sc1CeSZ layer close to the anode (Point D) is also evidently lower than 1 bar, which is beneficial for the relief of delamination. However, if YSZ is used as the barrier layer, the oxygen partial pressure at the interface between YSZ and the thinner GDC layer (Point D’) is slightly greater than 1bar. During the transient working process, this oxygen partial pressure may exceed the mechanical limit and causes delamination. If the potential at the point L is set to 0, the potentials at the interfaces can be determined according to the mathematical model. The results are shown in Figures 12C,D. The profiles for both cases are very similar with that of the oxygen partial pressures. The potential variations on the GDC layers are very small because the electronic conductivity of GDC is significantly greater than that of 10Sc1CeSZ and YSZ.