Conductivities and molar composition of the SOEC.
The hydrogen production by SOECs coupled with renewable energy sources is a promising route for the sustainability hydrogen economy. Multiphysics computing simulations appear to be the most efficient approaches to analyze the coupled mechanisms of SOEC operation. Using a relevant model, it is possible to predict the electrical behavior of solid oxide electrodes considering the current collector design. The influences of diffusion and grain diameter on cell performances can be investigated through 2D simulations, current–voltage characteristics, and current source distribution through electrodes. The simulation results emphasize that diffusion is linked to a relocation of the reaction away from the interface electrolyte/electrode, in the volume of the cathode. Furthermore, the current collector proves itself to be a great obstacle to gas access, inducing underneath it a shortage of steam. Inducing gradients of grain diameters in both anode and cathode drives the current sources to occur close to the electrode/electrolyte interface, thus decreasing ohmic losses and facilitating gas access. This approach shows the crucial importance of cathode microstructure as this electrode controls the cell response.
- hydrogen production
- steam electrolysis (SOEC)
- electrochemical modeling
Usual electrolyzers employ aqueous electrolytes (alkaline water electrolyzer), and the major drawback of the electrolytic hydrogen production is its high cost in comparison with the steam-methane reforming process . In addition to alkaline-based electrolysis, there are mainly two technologies of electrolyzer, which are currently considered, one based on proton-exchange membrane (PEM) and the other based on solid oxide (SOECs). The electrochemical decomposition of steam into hydrogen and oxygen offers two advantages over the low-temperature process currently in use. High temperature electrolysis (HTE) is more efficient than traditional room-temperature electrolysis for many reasons. First, some of the energy is supplied as heat, which is cheaper than electricity. Secondly, at high temperature the electrolysis is less energy consuming due to the low theoretical decomposition voltage. In the 1980s, solid oxide fuel cell (SOFC) developments allowed steam electrolysis investigation . Thus, several experimental studies on high-temperature electrolysis (HTEs) using practical SOFC cells  have been carried out. These experiments showed encouraging results. The solid oxide electrolysis cell (SOEC) is the reverse operation of the same SOFC cell structure. Both are composed of two porous composite ceramic electrodes surrounding a gas-tight electrolyte. SOECs can rely on the interest that SOFCs have received for the past decades and thus utilize similar technology and materials. Water is reduced at the cathode Eq. (1), releasing hydrogen and oxygen ions. After crossing the electrolyte, the ions are then being oxidized to form oxygen within the anode Eq. (2). Both half reactions are balanced by general Eq. (3):
Due to the high operating temperature, SOECs do not need expensive catalysts, but these must meet strict thermal, chemical, and structural requirements imposed by temperature, and hydrogen and oxygen partial pressures. However, performance still remains limited in electrolysis mode compared to fuel cell mode .
Comparing to PEM electrolyzers for which carbon monoxide is a poison, SOECs offer the advantage to allow the co-electrolysis of water and carbon dioxide to produce syngas [5, 6]. According to AlZahrani and Dincer  a SOEC system can achieve energy and exergy efficiencies of 53 and 60%, respectively. However, the high operating temperature (>1000 K) is still considered the major limiting factor of these device. Jiang  has shown that the hydrogen production by SOECs coupled with renewable energy sources is a promising route for the sustainability of energy in the future. The solution of a real commercially competitive SOC technologies is the materials: reliability and stability of the electrode and electrolyte materials under the reversible electrolysis and fuel cell operation modes. The reversible fuel cells have the ability to switch between electrolysis cell and fuel cell modes, and it is one of the foremost features that facilitate storing/generating energy in a cost-effective manner. The optimized parameters on the designs of fell cells or of electrolyzer, solely, would not necessarily result in high performance of regenerative devices because these devices differ in electrode kinetics, gas environment, heat generation, and chemical stability. It is well known that high-temperature operation of SOECs offers inherent advantages, in terms of thermodynamics and kinetics compared with low-temperature electrolysis. In this context reversible solid oxide cells (RSOCs) are still at an early stage of development [9, 10]. Unfortunately, there is a general consensus that the performance and stability of SOECs are inferior to those of SOFCs , which is mainly due to the high-temperature operation.
Computing simulation appears to be one of the most efficient approaches to analyze the coupled mechanisms of SOEC operation. It can predict the SOEC behavior under various operating conditions. Mathematical modeling is an essential tool in the design of SOEC cells, as it is important to understand the limiting process of steam electrolysis. Recent literature shows a significant research and development effort focusing on the modeling of SOEC. The models developed by Udagawa et al. [12, 13] strive to describe all significant processes affecting the performance of a unit cell. These authors proposed one-dimensional or pseudo-2D simulations based on a planar geometry and taking into account mass transport. Ni et al. [14, 15] have described the mass transport within the electrode along with the electrochemical kinetics. The principal results of these investigations [16, 17] lead to a parametric control of the SOEC operation. In addition, Jin and Xue [18, 19] have developed a 2D model for a planar SOEC. The simulation results lead to a better understanding of the internal mechanisms for regenerative SOFCs. This model has been used to study the delamination phenomenon at the oxygen electrode/electrolyte interface. Few micro-modeling investigations are applied to SOEC, and relevant electrochemical models are required to improve the micro-scale predictions [20, 21, 22, 23]. Nevertheless, the results of SOFC micro-modeling can be employed to appreciate, for example, particle size, graded or homogeneously distributed porosity, or composition influence on electrode performance.
In the subsection below a multiphysics model of SOEC has been built using the commercial software Comsol Multiphysics®. Electrochemical reactions within porous electrodes are described using the Butler-Volmer equation. Modeling is based on solving conservation equations of mass, momentum, and charge balance. Simulations allow the calculation of gas concentration, current density, and potential distribution within the electrodes (i.e., interconnects and electrodes). These simulations establish how porous electrode performance is affected by current collectors and electrode microstructural parameters. The model is then applied to engineer a design of the electrode structure or current collector configurations.
2. Model equations
In the present model, mass and charge transport phenomena coupled with chemical and electrochemical reactions have been investigated within the SOEC cell. This mathematical approach is based on classical SOFC assumptions, and thus the model should depend on operating conditions, intrinsic conductivities of materials, and geometric parameters such as porosity or grain size [21, 22]. Additionally, according to the 2D model of Kenney and Karan , the interconnects play a critical role. Thus, a 2D approach was performed in this work. In this study, a finite element method has been used to solve mass and charge balances including transport through porous media and electrochemical reactions within the porous electrodes. The set of resulting conservation equations has been solved using the commercial software Comsol Multiphysics®. In this computational approach, steady-state conditions have been imposed. This model of SOEC is based on the following assumptions:
Perfect current collectors (equipotential surface with perfect contact)
No contact resistances
Negligible convective flow through the electrodes
Ionic and electronic conductivities depending solely on temperature
Equations are detailed in this section one balance at a time.
2.1 Charge balance
The electrode material is a mixed electronic and ionic conductor. For modeling purposes, this electrode is considered as a porous gas diffusion electrode wherein the electrochemical reaction occurs at the triple phase boundary, i.e., at the interface between the electronic conductor, ionic conductor, and gas phase. The current in a porous electrode can be split into two parts: one part flowing through the ionic phase and the other through the electronic phase of the porous matrix. During electrochemical reactions, electrons are then transferred from the ionic phase to the electronic phase or vice versa. The transport of each kind of charges (e−, O2−) can be described using Ohm’s law. To account for charge transfer between electronic and ionic materials, a current source term (A m−3) is employed in the charge balance in Eq. (4):
The effective conductivity () depends on the material and the microstructure of each electrode. Their values can be computed using Eq. (5) , where S and M subscript are, respectively, ionic material and electronic material.
The bulk conductivities and molar fractions that have been used throughout this work are gathered in Table 1 . Nickel has been used as the electronic material at the cathode, whereas the anode was made of LSM-type perovskite [24, 25]. YSZ allows the transport of ions in both electrodes and the electrolyte. Within the electrolyte ceramic membrane, there are no current sources. Therefore, pure and dense YSZ electrolyte is considered () in Eq. (4). Similarly, current collectors are assumed to be ideal electronic conductors. The charge balance ensures that the current produced at the cathode is consumed at the anode. Additionally, in each electrode, the electronic current is the opposite of the ionic one. According to Costamagna et al. , the current source terms can be described by the classical Butler-Volmer expression Eqs. (7) and (8), and the current sources are expressed as follows:
|Nickel electronic conductivity, [S m−1] ||4.5 × 105|
|LSM electronic conductivity, [S m−1] ||1.6 × 105|
|YSZ ionic conductivity, [S m−1] |
|Cathodic volume fraction of nickel, ||0.4|
|Anodic volume fraction of LSM, ||0.5|
|Cathodic volume fraction of YSZ, |
Anodic volume fraction of YSZ, 
The selected parameters of Butler-Volmer equation remain valid at high fuel utilization and low value of hydrogen concentration corresponding to the SOEC mode. The current sources and therefore the expression of the electrochemical reactions are expressed by the following Butler-Volmer Eqs. (7) and (8), F being the Faraday’s constant (F = 95,485 C mol−1) and R the ideal gas constant (R = 8.314 J mol−1 K−1):
where are the overpotentials defined as the difference between the electronic potential on one hand and the ionic and equilibrium potentials on the other hand, displayed in Eq. (9). In the present model, hydrogen electrode is assumed to have a thermodynamic potential equal to 0 V:
Laurencin et al.  suggested that Butler-Volmer expression is suitable to describe electrochemical reaction involving the exchange of 2 electrons. The values for electrochemical geometric symmetric coefficients and are usually considered to be close to 0.5 in the literature [26, 28]. The exchange current densities are critical parameters to describe the current generated by the cell. Both parameters take into account the virtual specific surface area linked to the vicinity of the triple point boundaries (TPB) where the charge transfer occurs. These parameters have been deduced from previous work [16, 24] in order to observe usual SOEC performance, and are gathered in Table 2 .
Finally, the electrolyte material is YSZ, a suitable ionic conductor at high temperatures. The electrolyte potential is thus expressed by a classical Ohm’s law Eq. (4) without any current source term () within the electrolyte.
The Nernst law Eq. (12) is used to assess the open circuit voltage (OCV) at operating temperature T of 1173.15 with an E0 of 1.1 V which, in turn, is used to obtain the cell voltage through Eq. (11) with the parametric input of the polarization computation. Vpol has been defined in order to stabilize the numerical convergence near the OCV operating point. The partial pressures are according to the inlet feed.
2.2 Mass balance
As stated in the assumptions of the model, since electrodes are porous media with low permeability, mass transport is only due to diffusion process. Binary gas interaction and pore wall effect should be taken into account using the following general expression of the mass balance:
where is the mass source term, directly linked to the current sources by Faraday’s law as expressed below:
The binary diffusion coefficient of water in hydrogen was computed through Eq. (26), presented in appendix. Thus, at the cathode, only hydrogen and steam are taken into account. The effective diffusion coefficient takes into consideration the microstructure of the electrode according to Eq. (16) :
|Electrodes porosity, [−]||0.37|
|Electrodes grain diameter, [μm]||3|
|Electrodes tortuosity, [−]||4.8|
Usually, the dusty gas model considers Maxwell and Knudsen diffusion to describe the gas flow through a porous media. It was possible to suggest an approach equivalent to Fick’s law (Eq. (13)) by considering an effective equivalent diffusion coefficient , illustrated by Eqs. (19) and (20), with being the Fick diffusion form coefficient:
In the set of mass balance equations, the ideal gas law is applicable Eq. (21). This work has been done at a working temperature T of 1173.15 K and a pressure P of 105 Pa:
Similar modeling has been performed at the anode side to describe mass flow through the porous anode. O2 and N2 are substituted to H2O and H2, respectively. The binary diffusion coefficient is presented in Annex under Eq. (27). Eqs. (22)–(24) display the Knudsen diffusion coefficient, the effective diffusion coefficient, and the Fick diffusion form coefficient, respectively:
3. Simulation conditions
The commercial software Comsol Multiphysics® has been used to investigate the behavior of a simplified serial repeating unit (SRU) of SOEC. The modeled geometry will thus include both electrodes, the electrolyte, and the current collectors. The resulting set of conservation equations is solved using the commercial software. This section will focus on the description of how the model was implemented.
This study objective has been to investigate the influence of the electrode microstructure considering a realistic geometry of the SRU. Figure 1 displays the geometry that forms the basis of the model. For simplicity’s sake, a two-dimensional model of the SRU was used according to its cross section. The cell dimensions are gathered in Table 4 . A mapped mesh has been used to obtain a reasonable number of degrees of freedom allowing calculation convergence within an acceptable computation time. Its parameters are also listed in Table 4 . The studied geometry is a cross section of SRU which is perpendicular to the main direction of the gas flow. Considering that channels allow an ideal distribution of reactants on each electrode, it is possible to consider them as boundary conditions. This means that the produced hydrogen and oxygen are perfectly collected and exhausted from the SRU by the gas manifolds. Therefore, constant gas concentrations within the gas channels have been considered, in order to use this value as boundary conditions of mass balance.
|SRU dimensions||Mesh parameter|
|a [mm]||1||Degrees of freedom||63,982|
|b [mm]||2||Number of mesh point||13,757|
|c [mm]||1||Number of boundary elements||1352|
|d [μm]||40||Minimum element quality||0.5505|
3.2 Boundary conditions
This section presents the boundary conditions used to solve each balance. The symmetry existing in the geometry allows the solver to consider half of the mesh. The following Figure 2 summarizes the different boundary conditions defined. The present multiphysics problem uses partial derivative equations that need boundary conditions to be solved, such as electric and ionic charge balances Eq. (4) or the mass balance Eq. (13). Such conditions are expressed by the set of equations presented in Table 5 that are to be satisfied.
|Cathode inlet composition||[H2O]0 = 10.09 mol m−3 (90%molar)|
[H2]0 = 1.12 mol m−3 (10%molar)
|Anode inlet composition||[O2]0 = 2.35 mol m−3|
4. Results and discussion
The simulations obtained with the multiphysics model previously described have been developed to investigate the impact of diffusion on the SOEC performance and to quantify the location of current sources within functionally graded electrodes.
4.1 Control by diffusion
The diffusion phenomenon was investigated via four simulation cases, referred to with letters A to D. A is the reference case, based on the set of parameters and geometry as described in Tables 2 – 4 . This set was modified to give simulations B, C, and D. Case B considers a perfectly collected current without the collector pin: the boundary conditions of mass balance and electronic charge balance are merged. The effect of electrode thickness on gas diffusion and reaction distribution within the electrode was studied with case C considering an electrode thickness of 80 μm instead of the previous 40 μm (case A). Finally, case D displays a steam diffusion coefficient 10 times lower than in the standard case A. Figure 3 shows the results obtained for case A (left), taken as a reference for the discussion, and case B (right), where the collector pin has been removed.
The distribution of current sources (ibv ) appears to be neither continuous nor homogeneous. Figure 4 exhibits the steam concentration distribution along the electrode at varying abscissa. Whereas high current densities are reached in the vicinity of the interfaces electrode/electrolyte under the gas canals, a low fraction of the total current is produced under the current collector. The water concentration under the cathodic current collector is very low, and this shortage is close to depletion.
Consequently, most of the current source terms are located below the gas channel, and these sources are rather small below the current collector. Since most of the current is generated below the gas channels, the convergence of the electrons toward the collector pin causes hot spots where high current densities are observed. To separate the role of diffusion from other phenomena, such as conduction of charged species, additional simulations were performed. In order to complete the investigation of diffusion process, additional simulation cases B to D were performed. Consequently, as seen in Figure 3b , no water depletion is observed in case B when considering a perfectly collected current without the collector pin. The current–voltage characteristics for cases A, B, C, and D are presented in Figure 5 .
As expected, a homogeneous current collecting significantly improves the performance. Gaseous reactants and products are not impeded by the collector pin anymore. Decreasing the water diffusion coefficient induces a large shift of the current-tension characteristic toward lower currents and thus increases concentration overpotentials. According to Juhl et al. , increasing the electrode thickness improves the electrode performance for LSM/YSZ composite cathode of SOFC using thicknesses between 2 and 12 μm. Nevertheless, Virkar et al.  have computed that SOFC cathode overpotential increases or decreases with electrode thickness depending on the ionic conductivity of the composite. In the present simulations, both conductivities (electronic and ionic) and gas transport have been computed, and the competition between gas, ions, and electrons was investigated. The simulation results highlight a complex distribution of electrochemical active areas (see Figure 6 ). Due to the competition effect between gas (water) diffusion and ionic charge transport occurring in “countercurrent,” the current source terms are located close to both interfaces. In order to obtain better comprehensions, the current sources distributions have been analyzed. Figure 6 presents diagrams of Faradic currents (A m−3) along the cathode thickness for the studied cases. The electrode thickness has been divided into 10 layers from the electrolyte/cathode interface to the cathode/gas channel interface, and the ratio of the current sources to the whole current for each layer is plotted versus the electrode thickness. The results are expressed in percentage of the total current generation. Since the most interesting phenomena occurs under the edges of the collector pin, the current source terms are considered at abscissa x = 1 mm.
The competition between the transport of ions and the diffusion of gases is highlighted. Electron transport is not a limiting step as the electronic conductivity is about 104 times higher than the ionic one. When the water diffusion coefficient decreases, the competition between O2− and gas transport is distorted, and a relocation of current sources toward the cathode/gas channel interface occurs ( Figure 6 ). Therefore, gas access implies a higher reactivity of the electrode/gas channel interfaces than expected.
On all 40 μm thick electrodes, all layers exhibit a reactivity higher than 5% of the total generated current. However, the 80 μm electrode generates roughly 65% of the total current within the first 24 μm of electrode thickness and 30% within the last 16 μm ( Figure 6c ). In other terms, 50% of the electrode thickness is responsible for 95% of the current generation. The thinner the electrode, the more homogeneous the current sources seem to be.
Figure 7a shows the current source generation along the 80 μm electrode thickness located in the middle of the gas channel (abscissa x = 0 mm). It can be seen that more than 70% of the total current generation comes from the first 30% of the electrode thickness. That observation is coherent with Hussain et al.  who consider that the electrochemical reaction is occurring exclusively at the interface electrode/electrolyte. Figure 7a , b show that increasing electrode thickness leads to current being generated closer to the collector pin. This is compensated by a higher reactivity in the interface vicinity. Those edge effects, obvious in all cases, can be attributed to the electronic ohmic drop. Moreover, the ionic ohmic drop also controls the distribution of anodic current sources. Figure 7b displays the anodic equivalent to Figure 7a . Contrary to the cathode, no control of the faradaic currents by gas access exists at the anode, since the gas is being produced. Nevertheless, both electrodes display similar current sources distributions. Such observation remains valid for all simulations. Consequently the anodic current sources distribution is driven by electrical cathodic behavior.
4.2 Influence of graded grain diameter
It is generally accepted that the performances of composite electrodes as well as graded cathodes in SOFCs are largely governed by TPB number, mass transport, and ohmic drop. It is well known that the polarization resistance decreases [31, 32] when using graded electrodes. The improvement of the microstructure is one of the key parameters to reach high electrochemical performances.
In the context of current collecting considerations, this work has emphasized the effects of functionally graded electrodes on overall cell behavior (i.e., current collectors/electrodes/electrolyte). Thus, the influence of grain size distribution was investigated via the simulations gathered in Table 6 . A gradient of grain diameter distribution was introduced, leading to either the cathode or the anode presenting several layers of different grain diameters. The effects of such a change in microstructure on the currents obtained and their distribution were investigated. In this work, electrodes are constituted of two layers, each presenting a specific grain size. Their responses were modeled and the results analyzed. For each cathode and anode, five simulations were performed, referenced from I to V, with the subscripts “a” and “c” referring to anode and cathode, respectively. One goal of this work was to investigate the influence of the grain diameter on the reaction location and the cell performances. However, as Eq. (10) shows, the exchange current density is modified by the microstructure. To be able to compare the currents obtained for the different samples, the average exchange current density , given by Eq. (25) over the whole electrode thickness was kept constant:
|% thickness ω||[μm]||[μm]||Quote||% thickness ω||[μm]||[μm]||Quote|
To do so in multilayer samples, the increase in exchange current density caused by a layer composed of smaller particles was compensated by a larger particle layer. The different grain diameters were adapted according to the thickness of the layer so that the average exchange current density remains constant. The layer displaying the thinner grains is always closer to the electrolyte and will be referred to as layer ω and expressed as a percentage of the total thickness of the electrode.
The current–voltage characteristics obtained for the five cathodic and anodic cases are displayed in Figure 8a , b . On the cathode side, sample Vc seems to present an obvious performance increase if compared to the other modified cathodes. Additionally, case IVc shows higher current densities than case IIc, even if both electrodes are composed of grain of equal diameter in the layer closest to the electrolyte. On the other hand, the microstructural changes do not influence significantly the electrical anode behavior. The polarization curves can barely be differentiated.
Ni et al. [15, 33] have evaluated the potential of functionally graded materials for SOEC electrodes. These authors have compared  conventional nongraded electrodes with particle size-graded electrodes. For graded electrodes with a particle size decreasing by 50% from the gas/electrode interface (dg = 1 μm) to the electrode/electrolyte interface (dg = 0.5 μm), a negligible reduction in potential is observed in comparison with nongraded electrode. For a higher gradient of particle size decreasing by 70% (i.e., dg = 1 ➔ 0.3 μm) a significant saving of H2 electrode overvoltage has been observed. The present study exhibits similar results. In Figure 8 , the bilayer with a 50% decrease in particle size (Ic: dg = 4.048 ➔ 2 μm) showed poor improvement of performance. In addition, a 70% grain diameter decrease (Vc: dg = 4.95 ➔ 1.5 μm) led to a large enhancement of the current generated by the cell.
Thin particles at the electrode/electrolyte interface (cathode interlayer thickness) decrease the area-specific resistance (ASR). This increase in performance goes along with the relocation of the electrochemical reaction at the electrode/electrolyte interfaces. As shown in Figure 9 , 80% of the cathodic current arises from the first 10% of electrode thickness (case Vc), and on the anode side, more than 90% is generated by the first 20% (case IIa). The relocation in the volume observed for the smallest particles can easily be explained by the dependence of exchange current densities on grain diameter Eq. (10). Furthermore, forcing the current generation close to the electrolyte interface enables decreasing the ionic current path length.
The results show the influence of current collectors on gas access. Relevant control of material microstructures improves the diffusion of gaseous reactants and the current collecting. When diffusion is the limiting step, a relocation of the current sources within the volume of the electrodes is observed. On the contrary, if ionic ohmic drop becomes the rate determining step, current density sources are located close to the electrode/electrolyte interface. In some specific cases, the assumption that all electrochemical reactivity is located in the electrolyte interface vicinity cannot be made. That observation is emphasized when the competition between gas and ion transport is intentionally distorted, since a second reactive layer appears close to the cathodic gas channel.
It has been shown that it is possible to force the reaction to occur close to the electrolyte/electrode interfaces by layering the electrodes and introducing gradients of grain diameters. The obtained relocation is as high as 80% of the current being generated within the first 4 μm of the cathode thickness. The ohmic losses are reduced and gas access facilitated.
Contrary to the cathode, the changes of grain diameter gradient do not influence the electrochemical performance of the anode. However, its current production profile is consistently similar to the cathodic one. This means that specific attention should be given to the cathode microstructure.
Several aspects have been neglected in the present work and should be investigated to complete this approach and give global vision of the mechanisms that govern a SOEC response, e.g., the ohmic resistance due the dense ceramic membrane can be minimized using metal support technologies. In addition, the contact resistances shall be taken into consideration since they are a key parameter when optimizing the configuration of current collectors and will allow this model to be compared to experimental data of total SRU response.
Conflict of interest
No conflict of interest.
The diffusion coefficient for binary mixture of gases may be estimated from Fuller, Schettler, and Giddings relation with values coefficients for different molecules tabulated in :
|C||concentration [mol m−3]|
|D||diffusion coefficient [m2 s−1]|
|dg||mean grain diameter [m]|
|dp||mean pore diameter [m]|
|E||Nernst potential [V]|
|i||current density [A m−3]|
|L||electrode thickness [μm]|
|i0||exchange current density [A m−3]|
|M||molar mass [kg mol−1]|
|n||normal derivative vector [m]|
|Q||volumic current source [A m−3]|
|x||relative abscissa of electrode thickness [m]|
|Y||molar fraction or volume fraction [−]|
|α||factor of symmetry [−]|
|β||fick diffusion form coefficient [−]|
|Γ||mass source term [mol m−3 s−1]|
|σ||conductivity [S m−1]|
|ϕ||electrical potential [V]|
|Ω||ratio of electrode thickness [−]|
|i||type of conductor|
|k||type of species|
|a/c||anode or cathode|
|S/M||electronic or ionic conductor|
|LSM||strontium-doped lanthanum manganite La1-xSrxMnO3|