Novel porous electrode designs for reversible solid oxide hydrogen planar cell through multi‐physics modeling

A comprehensive multiphysics 3D model of an anode‐supported planar reversible solid oxide cell (rSOC) with a half‐channel‐unit‐cell geometry is created and validated. The physical phenomena that are modeled include reversible electrochemistry/charge transport, coupled with momentum/mass/heat transport. Several electrode microstructures comprising the homogeneous and functionally graded porosity distributions are applied to the validated model, to evaluate and compare the current‐voltage (j‐V) performance in both fuel cell mode and electrolysis mode. The results indicate that increasing the porosity in a homogeneous porous electrode does not always promote the cell's j‐V performance. An optimal porosity emerges where the effect of porosity on the mass transport is maximized, which ranges between 0.5 and 0.7 in the working conditions of the present study. Compared with homogeneous porous electrodes, the heterogeneous porous electrode design with a functionally graded porosity distribution is found to be a potential option to better the overall j‐V performance of the rSOC. Furthermore, it is discovered that theoretically grading the porosity in the width direction (i.e., increasing porosity from the center of each gas channel to the center of each adjacent rib) brings an outsize benefit on the cell's performance, compared to the traditional way of improving the porosity along the cell thickness direction.

INTRODUCTION [3]. The capability of operating sequentially between fuel cell mode and electrolysis mode makes rSOCs a promising technology for electrical energy storage systems [4].
Compared with other types of electrolyzers and fuel cells, such as low-temperature proton exchange membranes and alkaline electrolyzers, rSOCs have a higher tolerance for fuel impurity and are able to work at a higher reaction rate with a lower electrical power requirement (i.e., theoretically lower applied potentials for the same current density) [2]. However, key challenges persist when utilizing the conventional state-of-the-art SOFC (Ni-YSZ/YSZ/LSM-YSZ), ranging from practical cell performance (particularly with the electrolysis reactions). Other issues regarding degradation and durability as well as system integration must be addressed before the widespread adoption of rSOC [5]. The critical issues that are required to be addressed comprise high-temperature related challenges, handling of the dynamics of the rSOC and the associated switching issues, cell/stack design, system design, and so on. Amongst these, the design of the electrodes, particularly with respect to the structure and porosity distribution, to improve the overall performance of the reversible cell merits special consideration.
Surveying the literature shows that the performance of rSOCs is strongly dependent on porous microstructure distributions of both anodic and cathodic electrodes [6][7][8][9][10][11]. A number of studies focus on this aspect through either experimental or modeling work. Jung et al. [9] and Jung et al. [10] have studied the electrode materials through experiments, for low-temperature SOFCs and rSOCs, respectively. Jung et al. [9] found that Pt is the best material under 400 • C for both electrodes in terms of fuel cell performance among the three SOFC electrode materials being studied, porous Pt, Ni, and lanthanum strontium cobaltite (LSC), while LSC is a better cathode material over 450 • C. Jung et al. [10] have tested two different oxygen electrode materials, La 0.8 Sr 0.2 MnO 3 (LSM) and LSM/yttria-stabilized zirconia (LSM/YSZ), under both SOFC mode and SOEC mode for a hydrogen electrode supported cell, finding that the performance of both SOEC and SOFC modes are improving with the operating temperature. Furthermore, the LSM oxygen electrode is a better option for SOFC performance, while LSM/YSZ oxygen electrode is more durable for alternative SOFC/SOEC operating cycling.
Mathematical modeling has been proven to be a costeffective and reliable method for understanding fundamental mechanisms and optimizing designs of solid oxide cells at different levels [12]. Through accurate prediction and numerical analysis of the impact of the material properties on the cell performance, an optimal design may be achieved. Shi and Xue [11] presented a twodimensional (2D) SOFC model for the porous electrode microstructure design and optimization, which indicated that among the various porosity distributions investigated, an inverse parabolic distribution was shown to be especially promising in terms of the positive effects on cell voltage-current performance. Fashalameh et al. [7] presented a planar multi-layer anode-supported SOFC fabricated through slurry-based 3D printing. From their study, it was shown that the hierarchically macro mesopores can create higher power density for the same planar area. Furthermore, they believed that the next generation of SOFCs can be flexibly designed with various performance goals in mind by different microstructures, architecture, and features with thinner elements through 3D printing methods. Yan et al. [6] developed a graded-porosity cathode model, in which the influence of the microstructure on the activation, concentration, and ohmic overpotentials was investigated. They found that this numerical modeling is helpful for searching for an optimal porosity gradient profile for La 0.6 Sr 0.4 Co 0.2 Fe 0.8 O 3+δ (LSCF) cathodes. Zhang et al. [8] developed a single-channel multiphysics computational fluid dynamics (CFD) model for a pore phase tortuosity study. By investigating the relationship between the SOFC performance and microstructure parameters within the porous electrode, they found that the electronic tortuosity has a big impact on the SOFC performance, as well as establishing that there is an optimum electronic phase volume fraction. Following the above, it is observed that most of the present literature on the subject of electrode structure and porosity were predominantly focusing on the performance of SOFCs. For rSOCs, which have the same structure as SOFCs but are required to operate in SOEC and SOFC mode sequentially, the performance in the SOEC mode needs to be considered as well. If the optimal design for its SOFC working mode is unsuitable for the SOEC mode, then a compromise should be made in order to improve the overall performance in both modes. Furthermore, most of the studies on the microstructure distribution of the porous electrode were only considering the thickness direction; this is reasonable from both a fabrication and performance standpoint, as from their results it is shown that the microstructure distribution in this direction could make a significant impact in terms of the cell's j-V (current density-voltage) performance. However, it has not been established if the microstructure distribution(s) on the cell length and width directions would have a similar impact as in the thickness direction. If so, it follows that the optimal design of the cell's electrodes regarding the microstructure distribution could have more options, other than just along the cell thickness direction, especially if it can be realized using modern fabrication/manufacturing techniques.
This paper aims to investigate and predict the overall performance of the two working modes of a reversible solid cell, for the variation of the geometric distribution of porosity. To facilitate this study, a comprehensive 3D multiphysics model has been built using COMSOL (COMSOL Multiphysics 6.0, COMSOL Ltd.) and detailed in this study to simulate the two working modes (fuel cell mode and electrolyzer mode) of a single reversible solid oxide planar cell. Coupled with reversible electrochemistry, the present multiphysics model involves phenomena comprising mass and energy conservation through momentum equations, fluid transport through gas channels and porous media (working electrodes), and heat transfer.
To reduce the computation time and allow a better focus on the microstructure features of the cell, this model was built with a half-channel-unit-cell geometry based on the symmetry conditions of the planar cell. This model was validated through open literature [13], which has similar working conditions as the present study. The validated model then was used for the parametric study on the porosity distributions of electrode materials in both the Z direction (thickness) and XY direction (planar) across the cell, in both SOFC and SOEC operational modes. The multi-directional porosity study and its effect in reversible modes are, to the authors' knowledge, a first-in-literature. Functionally graded porosity distributions in three directions (x, y, z) and two directions (x, y) were trialed to find the optimal distribution for the improved j-V performance of the cell across different operating regimes. Subsequently, analysis and discussions of the simulation results are presented through the model properties, including the j-V characteristics, contributions from the kinetic, ohmic, and mass transport overpotentials as well as thermal and heat transfer properties. It is discovered that compared with the homogeneous porous electrode, the heterogeneous porous electrode with a functionally graded porosity distribution could be a potential option to optimize the overall j-V performance of the rSOC. Furthermore, theoretically grading the porosity in the width direction brings a larger benefit to the cell's performance than the traditional way of improving the porosity along the thickness direction.

METHODOLOGY
A multiphysics 3D model was built and validated first, on which several electrode microstructures by using different functionally graded porosity distributions were applied, to evaluate and compare the cell's j-V performance in both the SOFC mode and SOEC mode.
The overall methodology of the study is structured into two parts. First is the model description and numerical solution approach as well as a general/reference model val-idation. Subsequently, once validated, the model will be replicated and amended by adding two different porosity distribution cases to the electrodes, while the specifications and bulk geometry are exactly the same as the first model.

2.1
Model description

Model geometry
In the existing research, many 2D or even 0D models have been presented to describe the mass transport and electrochemical process(es) of the cell, obtaining generally good agreement with experimental results. However, for a study that is to predict the performance for a variation in layer thicknesses or microstructural features (such as porosity), a comprehensive 3D model is needed. This research presents a 3D model of an anode-supported planar cell with a new half-channel-unit-cell geometry based on the symmetry conditions of the planar cell. As shown in Figure 1, the electrolyte is sandwiched between the two porous gas diffusion electrodes (GDEs). The gas feed in the two electrodes is arranged in a counterflow configuration. Table 1 lists the geometry parameters used in this single-unit cell. When the cell is operated in fuel cell mode, the Ni-YSZ electrode of the cell is the anode, while in electrolyzer mode, it should be conventionally designated as a cathode. To avoid this confusion, hydrogen electrode (Ni-YSZ) and oxygen electrode (e.LSCF) are used to denote the respective electrodes instead, as only the solid oxide hydrogen cell has been studied in this research. The composite electrodes are assumed such that the electrochemical reaction sites are uniformly distributed in both hydrogen and oxygen electrodes.

Governing equations
The physical phenomena that are selected in this multiphysics model include reversible electrochemistry and charge transport, coupled with momentum/mass/heat transport. The electrochemical reactions act as the source terms for heat, mass, momentum, species, and charge. In order to predict and compare the cell's performance from an electrochemical viewpoint, solutions that combine the reversible electrochemistry with CFD and heat transfer module are used to evaluate the j-V performance, flow, and species distributions, as well as temperature across the cell. The physics parameters used in the model and their corresponding source references are given in Tables 2-4.
Oxygen electrode (air side) For SOFC mode, hydrogen diffuses from the fuel side channel into the porous hydrogen electrode where it combines with the oxygen ions coming from the oxygen electrode, to form steam and generate electrons. At the air side, oxygen diffuses from the air side channel into the porous oxygen electrode, where the oxygen molecules combine with electrons to form oxygen ions. The electrical current and power are generated through this process. For SOEC mode, the electrochemical reactions are reversed as shown in Equations (1) and (2)

Charge transport and electrochemical reactions
Both electronic and ionic transport is allowed in the two composite electrodes, whilst the electrolyte only allows ions to migrate through. Secondary current distribution is applied to define the transport of charged ions in an electrolyte of uniform composition using Ohm's law in combination with a charge balance. Butler-Volmer-type equations are applied to describe the relationship between charge transfer and overpotential. According to Ohm's law, the governing equations below are applied to describe the charge balance [17]: where and are the effective conductivities for the electrolyte and electrode phases, respectively, Φ and Φ are the corresponding potentials. Q and Q represent the current source or sink in the domain equations. In the present model, the electrolyte domain is considered to only conduct ionic current in the ion-conducting phase, while the TA B L E 4 The parameters presently used for the heat transfer model
The conductivity values used in the model are shown in Table 3. The conductivities of the material in each phase of the electrodes are taken from Ferguson et al. [20]. The electrode volume fraction is used to calculate the effective conductivity of the porous matrix [17], whereas Bruggeman's model is used for calculating the effective conductivity for both porous electrodes as shown below [21,22]: where σ eff, i , V i , σ i , and ε i represent the effective conductivity, volume fraction, material conductivity, and porosity respectively, for each phase of each electrode as shown in Table 3.
The operating cell voltage could be written as [17,23,24] cell = eq, − eq, − where E eq denotes the equilibrium potential, subscript O and H stand for the oxygen electrode and hydrogen electrode respectively, and η represents the overpotentials. The activation overpotential in the electrodes is defined as [21,22]: where stands for the open-circuit potential, which is equal to the electromotive force, given by the Nernst equation when the hydrogen-steam mixture is only considered on the fuel side [17,22]: where 0 is the temperature-dependent open-circuit potential at standard pressure (1 atm), p i is the partial pressure at triple phase boundary (TPB) in atm and T is the process temperature. Butler-Volmer type equations are used to describe the oxidation and reduction reactions and their voltagecurrent density relationships [13,14]: , where ,loc and ,loc represent local current density in hydrogen electrode and oxygen electrode respectively, α H and α O represent the charge transfer coefficients for hydrogen and oxygen electrode given in Table 2. F is Faraday's constant and R is the universal gas constant. 0 is the exchange current density, defined below (from Njodzefon et al. [13] and Leonide et al. [14]): where stands for the partial pressure at the TPB. E act, and E act, are the activation energies for the respective oxidation/reduction reactions, along with the prefactors and . The exponent terms , b, and m are empirical constants given in Table 2.
Mass transport H 2 and H 2 O (in steam/vapor form) are assumed to be transferred in the fuel side channel and the porous elec-trode, and O 2 , N 2 , and H 2 O are in the air side channel and the porous electrode in the present model.
The basic equation for the conservation of mass of each species i is given [11,25,26]: where ρ is the density of the mixture, u is the mass average velocity of the mixture, ω i the weight fraction of species i, j i is the mass flux relative to the mass average velocity, and R i is the rate expression of the source term.
Maxwell-Stefan diffusion model is more computationally expensive but also a more detailed diffusion model, compared with the simpler Mixture-averaged model. Considering the already simplified half-channel-unit cell geometry used in our model, the Maxwell-Stefan diffusion model is still applicable in terms of the computational cost, but Knudsen Diffusion is not accounted for (i.e., molecules colliding with walls). Hence the following equations are employed to describe the mass transfer for each constituent i [19,23]: where d k is the diffusional driving force, D ik is the multicomponent Maxwell-Stefan diffusivity defined in Equation (18), T is the temperature, is the thermal diffusion coefficient, and x k is the mole fraction.
where k d , T, p, M, and v are the reference diffusivity, temperature, pressure, the molecule molar mass, and kinetic volume respectively, of which values are given in Table 2.
The effective diffusivities are calculated by using Bruggeman's model as below:

Momentum equations
Momentum equations are solved to determine the fluid velocity and pressure. Brinkman equations are applied in the present model to describe the flows in the flow channels and the porous electrodes [12,27,28]: For the flow channels, )] For the porous electrodes, where ρ is the density, ε is the porosity, u is the velocity vector, p is the pressure, T is the absolute temperature, I stands for unit matrix, and F is the volume force vector. µ and k represent respectively the dynamic viscosity and permeability given in Table 2. Q m is the mass source or sink.

Heat transfer
Temperature is a critical factor that could have a significant impact on the rSOCs performance. Therefore, heat transfer phenomena are included in the present model, which are heat conduction and convection, coupled with electrochemical heating to define the domain and boundary heat sources based on the sum of irreversible (Joule heating and activation losses) as well as reversible heat. The energy conservation equation used in the present heat model [28][29][30]: where ρ is the density, C p is the heat capacity at constant pressure, k is the thermal conductivity, u is the velocity field, and Q is the heat source or sink. The thermal conductivities, densities, and heat capacities of the solid phase, which include the porous electrode matrix, electrolyte, and interconnected parts and ribs, are shown in Table 4. The fluid in the flow channels and porous electrodes is the mixture of species H 2 O and H 2 on the hydrogen electrode side and O 2 , N 2 , and H 2 O on the oxygen electrode. Their thermal conductivities and specific heat values are sourced from Celik and Akhtar et al. [17,19]. For the porous medium, local thermal equilibrium is chosen in the present model. The effective thermal conductivities of the solid-fluid system are calculated assuming that plane layers are parallel to the heat flow.

Boundary conditions
Two symmetry boundaries are set at the two sides of the present cell model, as this model is a unit half-channel model based on an axial symmetry condition. The electrical potential at the hydrogen electrode current collector, which is between the hydrogen electrode and the rib, is set to zero. The potential at the oxygen electrode current collector, which is between the rib and the oxygen electrode, is set as the cell working potential, . The electrical insulation is assumed around the outside walls of the two electrodes and the electrolyte.
The pressure difference between the inlets and outlets of the two flow channels (hydrogen and air) is defined as the boundary conditions for momentum transport, whose values are given in Table 2 for the model validation. Two pressure drop conditions for the hydrogen side channel are also applied in the porosity study, which are 40 and 70 Pa. At the hydrogen electrode side, two inflow volume fractions of hydrogen and steam are applied for each working mode, which are 80 % H 2 and 20% H 2 O for SOFC mode, 20% H 2 and 80% H 2 O for SOEC mode, and 50% H 2 50% H 2 O for both working modes. At the oxygen electrode side, the air is supplied. For the mass insulation boundary, it is assumed that no mass flows across the walls.
For the heat transfer model, thermal insulation is assumed on the top and bottom sides of the cell, and both ends except the inlets and outlets. The temperature boundary condition for the inlets is set to the operating temperature. The outlet boundaries are set as outflow for a convective-dominated heat transfer condition. closer to the electrolyte-electrodes interface and both of the inlets and outlets to improve the accuracy and convergence, while the coarser mesh for the remaining parts to save the computation cost, shown in Figure 1. The nonlinear PARDISO solver has been used, and the relative tolerance is set to 0.1% for the main variables and 0.01% for the current density value.

Model validation
This model is validated in comparison with experimental data from Njodzefon et al. [13]. As shown in Figure 2, the model is in generally good agreement with experimental data in the electrolyzer mode. In fuel cell mode, there is a slightly lower performance compared to the experimental data. The error in SOFC mode is about 30 mV, corresponding to 4.3% at 2.0 A/cm 2 , which is acceptable since a linearly increased error in the SOFC mode was observed in the compared source paper as well, which seems to be caused by specific temperature effects because of hydrogen generation consuming electrical and thermal energy in SOEC mode as discussed in the source paper. [13] It is also worth mentioning here that the validation is performed purely in order to establish that the fundamental physical parameters in the model are sound in general, prior to being applied in the dedicated microstructural studies.

Model with microstructure modification and porosity distribution
Two main porous electrode cases, homogenous and heterogeneous porous electrode designs in terms of different porosity distributions, have been studied in the present TA B L E 5 The parameters used in the porosity study The prefactor in the function (z) k2 The grading coefficient in the functions (y) and (y) k3  Table 5. Considering SOFCs are usually operated at temperatures of about 700-900 • C [31], the modeling temperature is set at 800 • C.

Homogenous porous electrode case
For a homogenous electrode case, an evenly distributed microstructure electrode is implemented, which means the porosity distribution is set to a constant number. Having a structure of a hydrogen electrode (anode) supported planar cell, the modeled cell has an oxygen electrode (30 µm) which is much thinner than its hydrogen electrode (400 µm). For a thin oxygen electrode, the porosity change does not have much impact on the cell performance [6]. For this reason, the parametric sweep of porosity values (0.1-0.9) was only set for the hydrogen electrode, while the oxygen electrode porosity was set at 0.8 constantly.
It is worth mentioning here that in reality, an rSOC is unlikely to have an electrode with porosity as high as 0.8. Also, some effective medium approximation models (i.e., Bruggeman's model) may have different forms or parameters in those extreme cases. However, such a big porosity range might be helpful for identifying the current vs. porosity trend in general. The pressure condition is not the focus of this study. Nevertheless, two different pressure drop conditions at the hydrogen electrode flow channel, 70 and 40 Pa, are modeled, considering that for a given porosity of the electrode, cell performance is also dependent on the operating conditions, typically the pressure drop across the electrode [11]. Two groups of mole fractions of H 2 /H 2 O are applied at the hydrogen electrode flow channel, 0.8/0.2, 0.5/0.5 (for the SOFC mode), and 0.2/0.8, 0.5/0.5 (for SOEC mode).

Heterogeneous porous electrode case
To facilitate a comparison between the homogenous case and different porosity distributions in the heterogeneous electrode case, functionally graded porosity distributions in three single directions (x, y, z) and two (x, y) directions were modeled to investigate the current performance under the two working modes of rSOC. The porosity functions used in heterogeneous electrode cases are shown in Figure 3 and described by Equations (25)- (31). The related parameters in the functions are shown in Table 5.
For the thickness(z) direction, Shi and Xue [11] found that among the linear and some nonlinear functionally graded porosity distributions, the inverse parabolic graded porosity distribution shows highly promising performance. which could be explained by that the inverse parabolic distribution allows the porosity decreases more rapidly near the electrode/electrolyte interface compared to the linear and parabolic distributions, and consequently facilitate the increase of the volumetric electrochemical reactive area, and thus promote the electrochemical reaction rate of the cell [11]. So, in the present model, the inverse parabolic function (z) is employed in comparison with other distributions, as shown in Figure 3a and Equations (25) and (26). The corresponding porosity distribution along the z direction is shown in Figure 4a.
where subscripts H and O stand for hydrogen and oxygen electrode, respectively. In the y direction, which is along the flow channel width, the linear function (y) (Figure 3b and Equation (27)) was set for both electrodes such that the porosity increases lin-early from the center of the flow channel to both sides of the unit cell as shown in Figure 4b. Grading the porosity in the width (y) direction could enable the porosity increase across the electrodes especially at the region under the ribs, where the mass flow resistance is believed to be higher than the region right under the gas flow channels, thus the overall mass transport could be improved.
Along the cell length in the x direction, the porosity distribution function is set as shown in Figure 3c and described by Equations (28) and (29). As shown in Figure 4c, the corresponding porosity distributions on the hydrogen electrode and oxygen electrode present a linear increase in the cell length from the inlet side to the outlet side to enhance the mass transport along the flow channel.
( ) = 4 ⋅ + 4 (28) In addition to the three unidirectional functionally graded porosity cases described above, electrode structures were also created whereby the porosity was graded simultaneously in two directions (x and y direction). The distribution is shown in Equations (30)

Examination of homogenous porous electrode case
The modified model with homogenous porous distribution is first utilized to investigate the variation of average current with the variation of porosity value. As can be seen from the current vs. porosity curves in Figure 5, the current densities increase with the increment of porosity value until reaching the maximum current densities which are obtained around the porosity values of 0.5 to 0.7 in both working modes; thereafter the current performance starts to drop with increasing porosity. This result for the SOFC F I G U R E 5 Variation of current density with homogeneous porosity when V cell = 0.75 V at solid oxide fuel cell (SOFC) mode and V cell = 1.25 V at solid oxide electrolyzer cell (SOEC) mode mode, which is denoted by the upper part of Figure 5 is in agreement with the result from Shi and Xue [11]; interestingly, it also shows the porosity has a larger impact on the cell's current on the SOEC mode, which is the lower part of Figure 5, especially for the case of higher pressure drop value (70 Pa).
In order to explain the presented current variation with the porosity value in Figure 5, the distributions of mass concentration, the variation of conductivity, and the average temperature in the homogeneous porous electrodes are investigated in both the SOFC (0.75 V) and the SOEC mode (1.25 V), shown in Figures 6-8. Four typical cases with four porosity values (0.1, 0.2, 0.5, and 0.9) are picked in each working mode according to the current variation with porosity value in Figure 5, to examine the mass concentration distribution in the electrodes. It can be observed that the respective concentration distribution of H 2 and H 2 O in SOFC and SOEC mode is improved from the porosity of 0.1-0.5, whereas increasing the porosity from 0.5 to 0.9 does not show any significant enhancement. Compared with the SOFC mode, this improvement is much more apparent in the SOEC mode, shown by the larger impact of porosity on the current in the lower part than the upper part in Figure 5. It is reasoned that generally speaking, the mass concentration of reactant (H 2 in SOFC and H 2 O in SOEC) could be the dominant factor leading to the current improvement in both cases within the porosity range of 0.1-0.5. When the porosity is around 0.5, where the mass concentration reaches its maximum, the optimal current is obtained as well. Beyond this value, mass concentration is no longer the dominant factor contributing to the current change. In addition, it is observed even in this As the oxygen electrode(the lower thin part) was hidden under the electrolyte and the hydrogen electrode(the upper thick part) originally, the image of the oxygen electrode is deformed from its actual position and moved downwards to get fully exposed. The same deformation is used in Figures 10 and 11. homogenous electrode scenario that there is a large asymmetrical impact of the conditions, in that the pressure and concentration increases had a much more outsize effect on the SOEC case. Other factors such as the variations of the conductivities and the average temperatures may keep affecting the cell's current performance.
Besides the well-known fact that the effective electrochemical reactive area decreases with the increase of electrode porosity, the electrolyte and electric conductivities will also decrease due to the increase of porosity (shown in Figure 7). This reduction of conductivities is proposed as the dominant factor that leads to the attenuation of the current after the improvement in reactant mass transport reaches its maximum. Furthermore, other F I G U R E 9 j-V performance with heterogeneous porosity case factors such as the average working temperature across the electrodes will vary as well with the different porosity values (shown in Figure 8), which might aggravate the current change, but is not of sufficient magnitude to offset the effect of the lower conductivity towards overall current density.

Examination of heterogeneous porous electrode case
The electrodes with functionally graded porosity distribution separately on the cell's length (ε(x)), width(ε(y)), and thickness (ε(z)), as well as simultaneously on the x and y direction are simulated (Figures 3 and 4). The cell's j-V performance in its two working modes is examined under two running conditions for each mode (50% and 80% H 2 for SOFC mode, 50% and 80% H 2 O for SOEC mode), shown in Figure 9.
According to Shi and Xue [11], the inverse parabolic function on the thickness direction demonstrated better performance over linear and parabolic functions generally. For this reason, in the present study, the inverse parabolic function is also applied to the thickness direction to compare with other three functionally graded porosity distributions on the cell length and flow channel width direction, ε(x), ε(y), and ε(x, y).
Interestingly, at SOFC mode working voltages (0.8-0.7 V), the linear function ε(y) surpasses the other three variations, showing a constant better performance (Figure 9, right). Moreover, this superiority of ε(y) over the other three in terms of current performance tends to be more obvious with the bigger current density range (above about 1.5 A/cm 2 ). In SOEC mode (Figure 9, left), within the current density of around -2.5 A/cm 2 , ε(y) is still ranked first in terms of the cell's current density from high to low, followed by ε(z), ε(x, y), and ε(x) within the current range of about -2.5 to -2.38 A/cm 2 and ε(x, y), ε(z), and ε(x) when current is below -2.38 A/cm 2 . When the current density is beyond -2.5 A/cm 2 , however, ε(z) yields a higher electrolysis current for the same voltage.
In order to understand this performance change with different porosity distributions of the electrodes, the corresponding mass concentration distributions in the SOFC mode ( Figure 10) and SOEC mode ( Figure 11) are employed. Generally, the ε(y) and ε(z) distributions show a more uniform mass concentration gradient in both hydrogen and oxygen electrodes than the ε(x, y) and ε(x) distributions, and ε(x, y) shows a better concentration distribution than the ε(x) for both working modes. In SOFC mode, a slightly better hydrogen consumption could be observed on the ε(y) distribution with a minimal hydrogen mole concentration of 0.16 (as shown in the legend), compared with the value of 0.21 for ε(z) distribution. At SOEC mode, however, ε(z) surpasses ε(y) showing a slightly better hydrogen concentration gradient. Accordingly, it is reasonable to estimate that the average mass flow resistance in these four porosity distributions could be ranked from high to low as ε(x), ε(x, y), ε(z), and ε(y) for SOFC mode, and ε(x), ε(x, y), ε(y), and ε(z) for SOEC mode. As a result, the current generated in these four cases could be in the same order from low to high, which exactly matches the cell's j-V performance shown in Figure 9.
A further comparison of current density between homogeneous and heterogenous cases is carried out under the same working conditions of 70 Pa for the pressure and 50% H 2 and 50% H 2 O for the inlet gas fraction. As shown in Figure 12, all the functionally graded porosity F I G U R E 1 2 Comparison of current density between homogeneous and heterogenous cases at V cell = 0.75 V for solid oxide fuel cell (SOFC) mode and V cell = 1.25 V for solid oxide electrolyzer cell (SOEC) mode distributions present an overall better performance over the homogeneous cases in terms of the current density, and ε(y) gives the largest current in the SOFC mode and the second largest in the SOEC mode. ε(y), the functionally graded porosity distribution on the flow channel width direction enables the porosity to increase rapidly across the electrodes, especially at the region under the ribs, where the mass flow resistance is believed to be higher than the region right under the gas flow channels, thus the overall mass transport is improved, and consequently benefits the average current generated.

CONCLUSION
A comprehensive multiphysics 3D model of an anodesupported planar rSOC with a new half-channel-unitcell geometry was built and validated, on which several electrode microstructures, including functionally graded porosity distributions as well as homogeneous porosity distributions, were applied, to evaluate and compare cell's j-V performance on both the SOFC mode and SOEC mode. From the simulation results, it is discovered that compared with homogeneous porous electrodes, the heterogeneous porous electrode design with a functionally graded porosity distribution could be a potential option to better the overall j-V performance of the rSOC. Furthermore, theoretically grading the porosity in the width direction brings an outsize benefit to the cell's performance, compared to the traditional way of improving the porosity along the cell's thickness direction.
From the view of electrode design, particularly with respect to the structure and porosity distribution, this study extends the previous studies which were predominantly focusing on the performance of SOFC mode and considering changing the porosity merely on the cell's thickness direction. It offers a multi-directional porosity study and its effect in reversible modes theoretically. From a manufacturing point of view, it might not be realistic yet to grade the porosity functionally no matter which direction. However, grading the porosity approximately linear is possible through modern fabrication and manufacturing techniques, such as 3D printing. Apart from the traditional way of changing the porosity in laminated layers along the thickness direction, grading or increasing the porosity in any way through the channel width direction, especially under the ribs may bring some surprising improvement in the cell's j-V performance. Of course, it needs to be further tested experimentally, which we encourage the rSOC community to do so.
Furthermore, although the functionally graded porosity on the cell length direction presents the least ideal performance among the four functionally graded distributions in terms of current density based on the simulation results, it is more likely to be realized and worthy of further testing if it can bring an extra improvement to the cell performance, considering its size, which is normally from several centimeters for a single planar cell, compared with the channel width or cell thickness which are only in an order of micrometers.
Lastly, from the electrode design point of view, in addition to optimizing the cell's j-V performance, other factors such as degradation, mechanical strength, etc. should also be considered. It is worth further testing if the functionally graded porous electrode has other benefits other than the current-voltage performance, which were outside the initial scope of this study, such as the efficiency (voltage as well as fuel utilization/H 2 yield), and dynamic behavior, among others.