Multiphysics Coupling Simulation and Parameter Study of Planar Solid Oxide Fuel Cell

In this paper, a numerical model of gas flow, heat transfer, mass transfer and electrochemical reaction multi-physics field coupling of a planar SOFC is established and solved. According to the calculation results, the distribution of velocity, temperature and concentration inside the SOFC cell is analyzed. The influence of cathode inlet flow rate, porosity, rib width and other parameters on the performance of SOFC is also discussed. The results show that within a certain range, increasing the cathode inlet flow rate can significantly increase the average current density of the cell. Increasing the porosity of the electrode can improve the gas diffusion of the porous electrode, thereby increasing the rate of the electrochemical reaction. Increasing the width of the ribs will result in a significant decrease in cell performance. Therefore, the rib width should be reduced as much as possible within the allowable range to optimize the working performance of the cell.


INTRODUCTION
As a new energy technology, anode supported planar solid oxide fuel cell (SOFC) has attracted much attention in recent years (Kong et al., 2016;Wang et al., 2018a;Wang et al., 2018b). Compared with other types of fuel cells, SOFC has the advantages of diversified fuel selection, high power generation efficiency, and no need for precious metal catalysis Presto et al., 2018;Zhang et al., 2018;Futamura et al., 2019;Hagen et al., 2019;Sahli et al., 2019). In order to further improve the performance of SOFC, a detailed understanding of its specific working process is required (Kong et al., 2020). The operation of SOFC involves a variety of physical processes such as heat transfer, material transport, fluid flow and chemical and electrochemical reactions. These physical processes are highly coupled and interact with each other, and their internal velocity field, temperature field and electric field are extremely complex. In general, experimental methods for SOFC are costly and time-consuming because of its multiscale feature of the porous electrodes. In contrast, the numerical simulation based upon computational fluid dynamics (CFD) method is much cheaper and more efficient, which can easily change various parameters of SOFC and analyze their impact on SOFC performance. Therefore, numerical simulation method for SOFC is widely used. Ferguson et al. (1996) established three-dimensional mathematical models of SOFC of flat and tubular structures, and studied and analyzed the effects of three different flow modes, namely, codirectional flow, reverse flow and cross flow, on electric field, temperature field and concentration field. Iwata et al. (2000) designed a quasi-two (co-flow and counter-flow) and quasi-three (crossflow) dimensional simulation program considering mass, charge and heat balance, and obtained the temperature and current density distribution along the flow direction and perpendicular to the electrolyte membrane. Jeon et al. (2006) established a two-dimensional reaction interval model containing a complete single channel, and optimized the width ratio of the flow channel and the interconnector through the research results. Stygar et al. (2012) established a 3dimensional corrugated SOFC model, and studied the influence of the unit geometry on the temperature distribution and heat transfer rate in the interconnector. The influence of co-current and counter-current on temperature distribution is considered in the simulation. The results show that the temperature gradient is relatively low and the temperature distribution is relatively uniform during counter-flow, which has a significant impact on the electrical conductivity, corrosion performance and durability of the connector. Yuan et al. (Andersson et al., 2012) comprehensively considered the physical processes such as heat transfer, mass transfer, fluid flow, charge transport, internal reforming reaction of fuel and electrochemical reaction, and coupled multiple physical fields to build a numerical model of SOFC single channel at medium temperature. The effects of ion and electron transport resistance in the electrode as well as operating temperature and excessive airflow on the cell performance were studied. Zhang et al. (2015) established an anode-supported threedimensional SOFC single channel model based on computational fluid dynamics (CFD), and compared the calculated results with the experimental results, which were well verified. The distribution of temperature, reactant velocity, current density and reactant concentration in different flow modes is explained in detail by comparing three different flow modes (co-, counter-and cross-flow). It is found that the output power of the current and counter-flow mode is much higher than that of the cross-flow mode when other operating parameters are kept constant. Dong (2019) established a complete SOFC multi-physics coupling model based on the open-source fuel cell code, and used the model to predict the performance of the cell, and discussed the transmission phenomenon and electrochemical characteristics inside the cell. Zeng et al. (2018) developed a three-dimensional model of the solid oxide fuel together with the optimized interconnection design. The results show that the current density and thermal stress are related to the shape of the interconnection tip and the depth of the cathode. Compared with rectangular and quadrilateral tip interconnects, triangular tip interconnects have the best electrochemical performance. With the increase of the tip depth, the current density will increase accordingly, and only the current density of the trapezoidal ribs decreases with the increase of the tip depth. Wang et al. (2018a), Wang et al. (2018b) established a threedimensional model of a planar anode supported solid oxide fuel cell, and the influence of working pressure on the cell performance was studied. The results indicate that increasing the operating pressure can improve the cell performance by increasing the open circuit voltage and reducing the activation overpotential, and enhance the electrochemical reaction near the electrolyte. Lee et al. (2017) presented a three-dimensional model of a reversible solid oxide fuel cell, and the effects of different geometric structures and operating parameters (electrode support layer thickness; interconnector rib size; fuel gas composition) on current-potential characteristics and round-trip efficiency were studied. The results show that due to the uneven distribution of the reactants, the size of the ribs has a significant impact on the cell performance, especially when the support layer is thin. Li et al. (2018) proposed a three-dimensional model that couples mass transfer, electron transport, and electrochemical reactions based on a solid oxide fuel cell with porous electrodes and gas channels, the results show that the current density distribution is significantly related to the gas composition distribution. Meanwhile, the simulation results also indicate that, the ridgeline of the interconnection line will limit the mass transfer inside the cell and negatively affect the electrochemical reaction.
However, the current research rarely involves the influence of the rib width parameter of the cell on the cell performance. Only the research by Li et al. (2018) considered the rib width parameter without considering the effect of temperature field in the simulation. The present work discussed the influence of rib width with models coupled with temperature field. In this research, a three-dimensional multi-physics model of planar SOFC is established to study the physical fields distribution in the cell. And the effect of some operation parameters such as cathode inlet flow rate, inlet temperature on the cell is discussed based on the model. Meanwhile, the cell performance under different structural parameters (porosity and rib width) was studied, the rib width is an important cell structural parameter.   (Iman, 2017

MODEL GEOMETRY
A three-dimensional planar SOFC cell is shown in Figure 1. It contains 10 fuel and air channels. The flow direction of air and fuel in this cell are counter-flow. This model was established for result verification and parameter study and its geometry parameters are selected based on the work by Khazaee et al. (Iman, 2017), as defined in Table 1.

Electrochemical Model
For hydrogen fuel cells, hydrogen and oxygen are fed into the anode and cathode channels of the cell respectively. When hydrogen and oxygen pass through the electrode surface, they diffuse to the electrode/electrolyte interface through the porous electrode structure, and the oxygen is at the cathode and electrolyte interface. Obtain electrons to form oxygen ions. The oxygen ions are transferred to the anode/electrolyte interface through the solid electrolyte layer, and chemically react with the hydrogen on it, and chemically react with the hydrogen on it to generate water vapor, while the oxygen ions lose electrons; The external circuit conducts from the anode to the cathode, forming a closed loop and generating current at the same time. The electrochemical reactions occurring at the anode and cathode of the cell are as follows: Anode : Cathode : The electrical energy output by the solid oxide fuel cell is converted from the chemical energy of the fuel. It can be known from the thermodynamic formula that when the temperature and pressure remain unchanged, the maximum output electrical work of the SOFC is equal to the decrease in Gibbs free energy. Thus, the open circuit voltage (EOCV) of SOFC can be obtained, that is, the Nernst potential is as the following: where E Nernst denotes Nernst potential(V), ΔG(T, P) denotes Gibbs free energy under the condition of temperature T and pressure P, P H 2 denotes hydrogen partial pressure at the threephase interface (Pa), P O 2 denotes oxygen partial pressure at the three-phase interface, P H 2 O denotes the partial pressure of water vapor at the three-phase interface.
In the actual working process of SOFC, there are some resistances that make the actual cell output voltage less than the ideal voltage E Nernst . This voltage loss phenomenon is called polarization. The main polarization losses are divided into three types, namely activation polarization, concentration polarization and ohmic polarization, defined as follows: where η act denotes activation polarization, η conc,a denotes concentration polarization of anode, η conc,c denotes concentration polarization of cathode, η ohm denotes ohmic polarization, and i denotes the local current density (A·m −2 ), i 0 denotes the reference exchange current density (A·m −2 ). i 0 can be expressed by Equation (8), where k e is the pre-exponential factor, E a is the activation energy.
The actual working potential of SOFC is defined as follows: where V op denotes the actual working potential. During the operation of SOFC, the charge is transferred in the electrodes, electrolytes and interconnectors in the form of electronic current and ionic current. The transport equations of electron current and ion current are expressed as follows: where S current denotes a general source term, k denotes an index that is el for the electrode or io for the electrolyte, A v denotes surface area to volume ratio/m 2 ·m −3 , φ k the potential (SI unit: V) and σ eff k denotes the effective conductivity (SI unit: S/m) which can be calculated as:  (Hajimolana et al., 2011).

Input parameters
Anode Cathode k e 6.54 × 10 11 2.35 × 10 11 E a /10 3 × J·mol −1 140 137 Frontiers in Chemistry | www.frontiersin.org January 2021 | Volume 8 | Article 609338 3 σ eff s,a 9.5 × 10 7 where V is the volume fraction for the specific materials in the electrodes and t is the structure-dependent tortuosity factors. Table 2 shows the input parameters of the electrochemical model.

Mass Transport
The material transport process is mainly composed of diffusion and convection. Factors such as concentration gradients, temperature gradients, pressure gradients, etc. will drive the transfer of gas components. At present, there are mainly three models used to describe the mass transfer process, namely the Fick model, the Maxwell-Stefen model and the Dusty gas model. The mass transfer model used in this paper is the Maxwell-Stefen model, which is expressed as follows: where ω i denotes the mass fraction of material i, D T i the thermal diffusion coefficient (assumed to be zero at this work), S i the source term, M i the molar mass of material i, and D eff ,ij is the effective diffusion coefficients, which can be calculated as follows: where D ij denotes the binary diffusion coefficient, and can be expressed by Eq. (16):

Momentum Transport
The Navier-Stokes equation is usually used to describe the flow process in the flow channel in SOFC. The combination of the N-S equation and the continuity equation can be expressed as: where ρ denotes density (kg m −3 ), μ the viscosity coefficient of fluid (kg m −1 s −1 ).
The momentum transfer in porous media can be described by the Darcy-Brinkman equation. The combination of Brinkman equation and continuity equation can be expressed as: where κ is the permeability of the porous medium (m 2 ), S i the source term.

Heat Transport
To accurately simulate the actual working process of the SOFC, thereby optimizing cell performance and improving thermal stability, the internal energy transfer process of the SOFC must be considered. For the solid area (electrolyte and connecting body), only heat conduction needs to be considered, the thermal energy equation can be expressed as: For the flow channel, heat conduction and convection heat transfer need to be considered at the same time, the thermal energy equation can be expressed as: For the porous electrode area, the local thermal equilibrium assumption is generally adopted, that is, the temperature of the porous medium material in the local area is equal to the temperature of the fluid, so the energy governing equation is: where C f is the molar concentration of fluid, C p,g the specific heat of gas mixture, ρ g the density of gas mixture, S the source term, and k eff effective thermal conductivity. The source term S for the thermal energy equation due to the electrochemical reactions can be expressed as:

Input Parameters and Boundary Conditions
The operating conditions and input parameters for the developed model are listed in Table 3. The potential at the anode current collector is set to zero and at of the cathode current collector is as the cell operating voltage. All other external boundaries and interfaces are insulated.

RESULT AND DISCUSSION
In order to prove the correctness of the model built in this paper, the simulated calculation results are compared with the experimental data provided by Yakabe et al. (2001). in the literature. The comparison results of the volt-ampere characteristic curves are shown in Figure 2. The results show that the calculation results of the planar SOFC model established in this paper are in good agreement with the experimental data given in the literature, and the maximum error is less than 10%. Therefore, the model can be used to study the performance of SOFC.

Temperature Distribution
The distribution of the temperature field inside the cell in the case of counter-flow is shown in Figure 3. It can be seen that the cell has a low temperature (1073 K) at both ends and a high temperature (around 1110 K) in the middle. This is because the fuel gas and air in the model are in a reverse flow mode, and the inflow temperature is both 1073 K. With the continuous inflow of fuel gas and air Inside the cell, the heat generated by the electrochemical reaction continues to accumulate, and the cell temperature rises accordingly, reaching a maximum near the middle of the cell. At both ends of the cell, the high-temperature gas flowing out of the cell exchanges heat with the gas flowing into the other side of the cell, and the temperature gradually drops.

Velocity Field Distribution
It can be seen from Figures 4A, B that the inlet velocity of the anode fuel gas is 0.4 m s −1 . After entering the cell, the boundary layer is formed first, and the flow velocity rapidly increases to about 0.58 m s −1 . At a relatively stable level, along the fuel flow direction, the hydrogen flow velocity in the anode flow channel hardly changes; while the cathode air inlet velocity is 3 m s −1 at the inlet. After entering the cell, the boundary layer is also formed. The flow velocity rapidly increased to about 4.5 m s −1 , and then the gas flow velocity in the cathode flow channel decreased continuously along the air flow direction. There is a difference in the velocity field distribution of the cathode and anode flow channels. This is because the anode consumes 1 mol of hydrogen and generates 1 mol of water vapor, while the total molar flow rate of the mixed gas inside the anode remains unchanged, so the velocity is almost unchanged; The continuous consumption of the flow direction causes the total molar flow of the cathode gas to be continuously reduced, and therefore the flow rate of cathode fluid is continuously reduced.

Concentration Field Distribution
The oxygen concentration distribution on the cross section of the cell are shown in Figure 5. It can be seen from Figure 5 that the oxygen concentration in the flow channel has been maintained at a high level (around 0.1). The oxygen concentration inside the porous electrode directly in contact with the flow channel is also relatively high, but the oxygen concentration in the cell ribs (the part of the porous electrode directly in contact with the interconnector) is very low (less than 0.01). There are three main reasons for this phenomenon: First, It is because the oxygen concentration in the flow channel is low  Frontiers in Chemistry | www.frontiersin.org January 2021 | Volume 8 | Article 609338 5 (the oxygen inlet mole fraction is 0.14); secondly, the rib area does not directly contact the flow channel, so its oxygen transportation mainly depends on diffusion, and the binary diffusion coefficient of oxygen is lower than other gas in SOFC, which is not conducive to the diffusion of oxygen, and for common anode-supported SOFCs, the thickness of the cathode is generally thin, which will also make the diffusion of oxygen more difficult. The imbalance of the oxygen concentration distribution will lead to the imbalance of the electrochemical reaction, resulting in uneven temperature distribution inside the cell. Figure 6 shows the current density distribution on the cross section of the flow channel. The arrow direction represents the direction of current flow. It can be known that the current density distribution on the entire cross-section is relatively uniform and the values are relatively small (less than 1 × 10 5 A m −2 ), but near the contact surface of the flow channel and the interconnector, the current density has a maximum value around 7 × 10 5 A m −2 . This is because the electronic current in the electrode can only pass through the interconnector part, and always takes the shortest path, but the flow channel is not conductive, therefore, the electronic current generated by the porous electrode part will be transmitted laterally to both sides of the

Parameter Study
The change of SOFC's working parameters and structural parameters will affect the distribution of physical fields inside the cell. This part mainly studies the influence of the several parameters on the performance of SOFC. The research results are calculated under the condition of output voltage V op 0.7 V.

Influence of Cathode Inlet Flow Rate
In order to study the influence of the cathode inlet flow rate on the performance of SOFC, the multi-physical field distribution of the cathode flow rate from 1 to 5 m s −1 was calculated. Figure 7 shows the change in current density when the cathode inlet flow rate changes (the standard operating condition is the cathode inlet flow rate u c 3 m s −1 ). It can be seen that when the cathode flow rate increases from 1 to 2 m s −1 , the current density increases rapidly from around 2700 A m −2 to around 3200 A m −2 . When the cathode inlet flow rate continues to increase, the increase in current density gradually slows down. This is because the initial air inlet flow is too small, so the oxygen in the back of the cell is exhausted. The main reason for limiting the cell performance is insufficient supply of reactants. At this time, increasing the cathode air flow will increase the current density significantly. As the cathode inlet flow rate continues to increase, the air on the cathode side gradually becomes excessive, but the electrochemical reaction in the cell has already obtained enough reactants, so the improvement of cell performance is not obvious, When the inlet flow velocity of Cathode increases from 2 to 5 m s −1 , the current density only increases by less than 100 A m −2 . Figure 8 shows the distribution of oxygen concentration in the flow channel under different cathode inlet flow rates. It can be seen that the lower the cathode inlet flow rate, the longer the oxygen stays in the cell, and the more oxygen is consumed. When the cathode flow rate is 1 m s −1 , the downward trend of the oxygen mole fraction at the outlet has a significant slowdown, and is close to 0, indicating that the flow rate of oxygen is not enough under this flow rate condition. On the other hand, with the continuous increase of the cathode inlet flow rate, the consumption of oxygen in the flow channel is getting lower and lower, and Figure 7 shows that the current density is continuously increasing, which shows that after sufficient reactants are provided, increasing the oxygen flow rate mainly improves the diffusion process of oxygen in the porous medium, thereby increasing the current density and improving the cell performance.

Influence of Inlet Temperature
The electrochemical reaction process in the SOFC working process is greatly affected by temperature. In order to study the influence of the inlet temperature on the performance of the SOFC, In the range of 873-1273 K, every 100 K is selected as a   working condition point, a total of five inlet temperature conditions are solved separately. Figure 9 shows the corresponding relationship between inlet temperature and current density. It can be seen that the current density increases significantly as the inlet temperature increases (from around 500 A m −2 to around 6500 A m −2 ). This is consistent with the formula given in the electrochemical model, that is, as the temperature rises, the local current density will increase, thereby increasing the average current density of the cell.

Influence of Porosity
Porosity is the microstructure parameter of the porous electrode, which reflects the difficulty of gas transfer in the porous medium. Figure 10 describes the corresponding relationship between current density and porosity. It can be seen that as the porosity continues to increase, the current density also increases (from around 2900 A m −2 to around 3200 A m −2 ). This is because the porosity improves the diffusion characteristics of the electrode, making the diffusion of oxygen and hydrogen inside the porous electrode easier, thereby improving cell performance.

Influence of Rib Width
The different rib width parameters selected for numerical calculation are 0.5, 0.75, 1.25, 1.5 mm. Figure 11 shows the corresponding relationship between rib width and current density. It can be seen that as the rib width increases, the current density gradually decreases (from around 3800-2600 A m −2 ). In order to understand the gas transport situation of the cell ribs, a straight line is selected as the calculation area on the center plane of the cathode ribs, and the distance between the straight line and one side of the rib is 0.1 times of the whole rib width. For different rib widths, the oxygen concentration distribution on the straight line is shown in Figure 12. Generally, as the rib width increases, the average oxygen concentration in the rib area becomes lower. When the rib width is 0.5 mm, the oxygen mole fraction at the cathode exit is about 3%, and when the rib width is other values, the oxygen mole fraction at the cathode exit is very close, both being about 1%. This is because when the rib width is greater than 0.75 mm, the oxygen concentration near the exit is too low, which is close to the lower limit of concentration preset in the simulation. It can be considered that the oxygen has been completely consumed at this time and has no contribution to the electrochemical reaction. Meanwhile, the larger the rib width, the lower the oxygen concentration at the entrance, that is, the smaller the total contribution of the entire rib area to the cell performance. Therefore, increasing the rib width is equivalent to increasing the area of the adverse reaction zone, on the other hand, increasing the rib width makes the oxygen transmission path in the porous electrode longer, and it is more difficult for oxygen to transport to the middle of the rib. In addition, increasing the rib width will also lengthen the current transmission length in the cathode, which increases the ohmic polarization. Combining the above three factors, increasing the rib width will reduce the performance of the cell. Correspondingly, the performance of the cell can be improved by reducing the rib width. However, the rib width is restricted by the manufacturing process and cost. Considering all factors, reducing the rib width as much as possible can improve the working performance of the cell.

CONCLUSION
In this paper, a numerical model of multi-physics coupling of a planar SOFC is established and solved. According to the calculation results, the distribution of various physical parameters inside the SOFC cell is analyzed. The influence of cathode inlet flow rate, porosity, rib width and other parameters on the performance of SOFC is discussed with a fully coupled models. The following conclusions are mainly drawn: (1) Due to the thinner cathode porous electrode and the smaller oxygen binary diffusion coefficient, the oxygen concentration in the electrode area directly in contact with the flow channel is much greater than the oxygen concentration in the electrode area (rib) directly in contact with the interconnector. The concentration varies greatly between the electrode area covered by channel and the electrode area covered by the ribs. This uneven concentration distribution will affect the electrochemical reaction rate; the current density distribution inside the SOFC is not uniform, and its maximum value is located near the interface between the cathode ribs and the flow channel (2) When the cathode inlet flow rate is less than a certain value, the oxygen flowing into the cathode will be exhausted before it reaches the cell outlet. In this range, increasing the cathode inlet flow rate can significantly improve the working performance of the cell. When the cathode inlet flow rate is greater than the critical value, the inlet air can meet the needs of the electrochemical reaction, and it is difficult to increase the inlet flow rate to significantly improve the cell performance.
(3) Increasing the porosity can improve the gas diffusion inside the porous electrode and reduce the concentration polarization, thereby increasing the output voltage of SOFC; Increasing the rib width will extend the gas diffusion path and the electron transmission path in the porous electrode. In order to improve the working performance of the cell, the rib width should be reduced as much as possible within the range allowed by objective conditions such as the manufacturing process.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.

AUTHOR CONTRIBUTION
ZD responsible for the overall organization, conception and writing of the paper. XS responsible for model construction, data analysis, and paper writing. JM responsible for data processing and assist in writing of the paper, ZJ responsible for data processing and assist in writing of the paper. GX responsible for paper conception and method guidance.

SYMBOL DESCRIPTION
A v Surface area to volume ratio/m 2 ·m −3