Effect of Electrolyte Thickness on Electrochemical Reactions and Thermo-Fluidic Characteristics inside a SOFC Unit Cell

: We investigated the effect of electrolyte thickness and operating temperature on the heat and mass transfer characteristics of solid oxide fuel cells. We conducted extensive numerical simulations to analyze single cell performance of a planar solid oxide fuel cell (SOFC) with electrolyte thicknesses from 80 to 100 µ m and operating temperatures between 700 ◦ C and 800 ◦ C. The commercial computational ﬂuid dynamics (CFD) code was utilized to simulate the transport behavior and electrochemical reactions. As expected, the maximum power density increased with decreasing electrolyte thickness, and the difference became signiﬁcant when the current density increased among different electrolyte thicknesses at a ﬁxed temperature. Thinner electrolytes are beneﬁcial for volumetric power density due to lower ohmic loss. Moreover, the SOFC performance enhanced with increasing operating temperature, which substantially changed the reaction rate along the channel direction. This study can be used to help design SOFC stacks to achieve enhanced heat and mass transfer during operation.


Introduction
Recently, the demand for a solid oxide fuel cell (SOFC) as a power source has increased with the development of portable small electronic devices [1,2]. Various studies have been conducted to miniaturize SOFCs while maintaining their promising high-power output [3]. Small fuel cells are projected to replace conventional batteries for portable applications due to improved power density [4]. In general, many researchers are interested in changing the configuration of the membrane electrode assembly (MEA) to find out the interrelationship with SOFC performances [5]. Electrodes and electrolyte geometries were basic design constraints to control the ohmic resistance inside SOFC [6]. In particular, making electrolytes thinner significantly reduces the ohmic loss of SOFCs and results in good fuel cell performance [7]. The interconnect configuration was also changed to evaluate cell performance to determine the geometrical effects [8]. Some researchers are trying to develop the high energy density power source which could be changed through various approaches [9]. Several models for operating conditions were proposed to reach the maximal efficiency of the SOFC system [10]. The design and optimization of SOFCs were conducted to understand the relationship between the parameters and performance when using the numerical model [11]. However, it is a time-consuming Energies 2018, 11, 473 2 of 15 process to experimentally identify the influence of geometry and operating conditions on each parameter. Moreover, it is difficult to accurately measure the thermo-fluidic characteristics inside the SOFCs due to the sealing problem and the high operating temperature [12].
As an alternative way to study the thermo-fluidic characteristics, numerical studies have been widely utilized to predict the performance of SOFCs by calculating the heat and mass transfer through electrochemical reactions [13]. Numerical modeling with a simplified configuration of SOFC stacks was conducted for comparing the performance curve with experimental data [14], and optimization of channel design was performed to obtain more efficient SOFC cells through a computational fluid dynamics (CFD) approach [15]. The performance of a SOFC can be affected by the channel configuration, which may change the electrode/electrolyte interfacial area. Lee et al. developed a high-fidelity physical model to investigate the effect of fuel utilization on the thermo-fluidic characteristics inside SOFC stacks [16]. They selected advanced models reported in the literature for numerical analysis, so we employed these models for the CFD simulations in this study. Shi et al. conducted isothermal mathematical modeling of anode-supported SOFCs to investigate the electrochemical processes [17]. They assumed simplified models, which must be improved. Jeon provided a comprehensive CFD model to predict SOFC performance with different operating conditions [18]. Although they obtained useful results, it is necessary to extend the dimensions because the two-dimensional model has limitations in explaining the effects in the span-wise direction. Shi et al. evaluated the SOFC performance with different micro-flow channel designs and compared fuel/gas pressure losses by using CFD analysis [19]. A comparison of the performance changes with porosity (based on the manufacturing process) was provided as a sensitivity test in this study. He et al. studied the correlation factors and effective parameters in their computational model [20]. Although various numerical analyses have been performed to estimate SOFC performance until now, an integrated database established using a parametric study of various geometries and operating conditions is still needed.
There are various reasons why such a parametric study should be conducted as follows: the configuration of each SOFC is different from the others, several types of software exist for numerical analysis, and many numerical approaches have employed different models even for the same physics. Therefore, we chose a basic rectangular type channel and SOFC unit cell validated with experimental data from the literature. Commercial finite element analysis software (COMSOL Multiphysics v.5.3, COMSOL, Inc., Burlington, MA 01803, USA) for CFD simulation was utilized to focus on the electrochemical reactions inside the SOFC stacks. High fidelity was achieved by using a higher-order model to simulate the physical phenomena based on a summary from the literature. Thus, our goal in this study is to predict thermo-fluidic characteristics and electrochemical reactions inside the SOFC stacks with the use of CFD simulations. In addition, we investigated the effect of the electrolyte thickness and operating temperature on the output power and air-fuel reaction to carry out a quantitative evaluation.

Model Geometry and Physics
In the present study, a three-dimensional non-isothermal model of a planar SOFC was developed. The model geometry is based on the literature and is shown in Figure 1a. The unit cell is composed of a strontium-doped lanthanum manganite/yttria-stabilized zirconia (LSM-YSZ) cathode, yttria-stabilized zirconia (YSZ) electrolyte, nickel/yttria-stabilized zirconia (Ni-YSZ) anode, air/fuel channels and a metal interconnect. A counter flow configuration was adopted for the air and fuel in the gas channel. Figure 1b shows the computational domain and boundary conditions with the dimensions summarized in Table 1. The model includes charge, mass, momentum, species, energy conservation equations, and electrode kinetics used to consider the electrochemical reaction. Some assumptions used in this study are as follows. The SOFC unit cell considered in this study operates at steady state. The porous electrode is made of isotropic and homogeneous materials [21]. The only gas species participating in the chemical reaction on the anode side is hydrogen (i.e., the reforming Energies 2018, 11, 473 3 of 15 reaction was not considered) and all gas inside the SOFC is treated as an ideal gas [21]. The flow regime in the gas channel is laminar due to a low Reynolds number based on the hydraulic diameter of the channel as a characteristic length (Re D = 168) [21]. Material properties are calculated based on operating temperature. The electrodes can be regarded as transparent optical materials, and the radiation heat transfer effect is negligible [22,23]. The local temperature equilibrium hypothesis for a porous medium (T f = T s = T) is applied. This means that the local temperature of the solid and gas phase in the porous electrode is the same [24]. Activation energy is constant in the temperature range from 700 to 800 • C. The electrolyte is a fully dense and solid material, and the electrode is a porous material [24]. channels and a metal interconnect. A counter flow configuration was adopted for the air and fuel in the gas channel. Figure 1b shows the computational domain and boundary conditions with the dimensions summarized in Table 1.The model includes charge, mass, momentum, species, energy conservation equations, and electrode kinetics used to consider the electrochemical reaction. Some assumptions used in this study are as follows. The SOFC unit cell considered in this study operates at steady state. The porous electrode is made of isotropic and homogeneous materials [21]. The only gas species participating in the chemical reaction on the anode side is hydrogen (i.e., the reforming reaction was not considered) and all gas inside the SOFC is treated as an ideal gas [21]. The flow regime in the gas channel is laminar due to a low Reynolds number based on the hydraulic diameter of the channel as a characteristic length ( Re D = 168) [21]. Material properties are calculated based on operating temperature. The electrodes can be regarded as transparent optical materials, and the radiation heat transfer effect is negligible [22,23]. The local temperature equilibrium hypothesis for a porous medium (Tf = Ts = T) is applied. This means that the local temperature of the solid and gas phase in the porous electrode is the same [24]. Activation energy is constant in the temperature range from 700 to 800 °C . The electrolyte is a fully dense and solid material, and the electrode is a porous material [24].

Charge Conservation Equation
There are two charge carriers, electrons, and ions in the SOFC. The conservation equations for each charge carrier are as follows: where φ and σ are the ionic/electronic potential and the ionic/electronic conductivity, respectively. The ionic/electronic conductivities of the constituent charge conductors can be obtained from [25].   where ϕ and σ are the ionic/electronic potential and the ionic/electronic conductivity, respectively. The ionic/electronic conductivities of the constituent charge conductors can be obtained from [25].
where T is the gas temperature. We used the effective properties to consider the effect of pore domains in porous media. Effective ionic/electronic conductivity was estimated by using the volume fraction of an electron conductor and the porosity of the electrode with various values depending on the manufacturing process [16]. σ where θ is the volume fraction of electronic conductor and ε is the porosity. Q ion and Q elec represent the charge carrier source terms. These terms only appear in the porous electrode domain where the electrochemical reaction occurred.
where S a is activation area per unit volume, and i means the local current density generated by electrochemical reactions.

Mass and Momentum Conservation Equation
In this study, the computational domain included the fluid flow channels and porous media. The general Navier-Stokes equation was applied to the fluid channels, and the Brinkman equation was employed for the porous media. The Brinkman equation was utilized to compute the fluid velocity and pressure fields of porous media for a single-phase flow in the laminar flow regime: where µ is the viscosity of the fluid, and ψ is the viscous stress tensor. Small variations in the density caused by the temperature variations of the fluid were considered in the momentum equation. The density of gas mixtures was estimated by the ideal gas law as follows: The dynamic viscosity of each gas species was calculated as a function of temperature, and the molar average was obtained to account for the mixture properties. where a i is the species dependent parameter listed in [20], and x i is mole fraction of species i. The equation shown below represents the momentum loss of the fluid passing through the porous medium.
where k is permeability. This source term was only applied to the electrode domain, which has a porous structure.

Species Conservation Equation
There are several models which can describe the diffusion phenomena. The Maxwell-Stefan equation can accurately account for multi-component diffusion, but it cannot consider the effect of Knudsen diffusion within the porous electrodes, which have micron-sized pores. The Dusty gas model can consider the effect of Knudsen diffusion of the multi-component gas mixture, but it has an implicit form and is very difficult to use for the three-dimensional numerical analysis. Therefore, we used the modified Fick's law, which has a relatively simple form compared to the Dusty gas model and can consider multi-component diffusion and Knudsen diffusion at the same time. There are four kinds of gas species (oxygen and nitrogen on the cathode side, hydrogen and water vapor on the anode side) in the computational domain. Conservation equations of each species i are as follows: where ω i is the mass fraction of species i, and D is the diffusivity. The diffusion flux, J, is evaluated using modified Fick's law. In the gas channels, the binary diffusivity of species i and j is calculated using the following formula [16]: where M i is the molar mass, and ν i is diffusion volume. The averaged diffusivity of the mixture was calculated by [16]: There are two types of diffusion mechanisms in porous material: molecular diffusion and Knudsen diffusion. Molecular diffusion is related to collisions between different types of gas molecules. Knudsen diffusion is related to a collision between the gas molecules and the pore wall, which occurs when the scale length of a system is comparable or smaller than the mean free path of the particle involved. Since the pore size of the porous electrode was the same as the order of the mean free path of a gas molecule, Knudsen diffusion was considered in this study.
where r pore is the pore radius. The diffusivity of each gas species in the porous electrode having an irregular pore distribution was estimated by calculating the effective diffusivity as follows [16]: where ε is porosity, and τ is tortuosity of a porous material. R is the species source term regarding electrochemical reactions, as shown in Equations (22)- (24): where M is molar mass, F is the Faraday constant, S is an activation area per unit volume, and subscripts a and c mean the anode and cathode, respectively.

Energy Conservation Equation
The general form of the energy conservation equation can be written as: where c p is the specific heat capacity and k is thermal conductivity. In the solid domain, the advection term may be neglected. In the fluid domain, whole terms should be considered. The energy conservation equation defined in porous media domains correspond to the convection-diffusion equation with averaged thermodynamic properties to account for both solid matrix and fluid properties. It is valid when the temperatures of the porous matrices and the fluid are in equilibrium. The specific heat capacity and thermal conductivity for each gas species are as shown in Equations (26) and (27), and the coefficients b i and c i refer to [26].
Molar averaged specific heat capacities, and thermal conductivities for gas mixtures can be calculated as follows: The effective thermal conductivity of porous electrodes can be calculated by following the mixture rule [27]: This rule applies by multiplying the first term of right side by the solid volume fraction, 1 − ε, multiplying the second one by the porosity, ε, and summing these resulted terms. The heat source term, which can be represented by the sum of Joule heating, polarization overpotential losses, and electrochemical reactions are as shown below [27]: where η act is the activation overpotential, E eq is the equilibrium potential of the electrodes, and σ is the ionic/electronic conductivity of the charge carriers. Due to the irreversible heat dissipation by transport of ionic/electronic charge carriers, Q joule is applied to the electrode/electrolyte and electrode/interconnect interface. Q act is a heat source term due to the activation overpotential, which is irreversible, and Q chem represents the heat sourced by exothermic/endothermic chemical reactions on the three-phase boundary in the electrode. These two terms, Q act and Q chem , are applied at the anode and cathode. The compression and viscous dissipation are neglected due to low flow velocity.
Since the radiation heat transfer is thought to be very small compared to other thermal energy sources, the radiation effect is ignored in this study.

Electrochemistry
The rate of the electrochemical reactions can be described as the reaction rate to the activation overpotential. For an electrode reaction, the activation overpotential can be defined as follows [27]: where ϕ elec is electronic potential, ϕ ion is ionic potential, and E eq is equilibrium potential. The mathematical relationship between the activation overpotential and current density in porous electrodes is described by electrode kinetics using a Butler-Volmer model, which can be written by: where i 0 is the exchange current density, S a is the active area per unit volume in a porous electrode and α is a charge-transfer coefficient. The exchange current density was calculated using the following Equation [28]: where E act and k e are activation energy (137 kJ/mol for the cathode, 140 kJ/mol for the anode) and pre-exponential factor (2.35 × 10 11 for the cathode, 6.54 × 10 11 for the anode), respectively. The exchange current density was calculated using the above equation with validated data from [29].

Numerical Details
All the operating and boundary conditions were determined from the literature [29]. Figure 1b presents the grid system and boundary conditions. The boundary conditions for solving the governing equations are as follows. The condition of the top and bottom interconnects was estimated as an electrical potential from 0.5 to 0.9 V and electrical ground, respectively. A pressure inlet condition was employed at the air channel inlet with velocity of 3 m/s, air (O 2 :N 2 = 0.15:0.85) at 800 • C. For the fuel channel inlet, the pressure inlet condition was a velocity of 0.4 m/s, fuel (H 2 :H 2 O = 0.4:0.6) at 800 • C. The outlet of both channels was at a pressure of 0 Pa. Span-wise side walls were considered as insulation. Initially, all the systems are regarded as being filled with a gas at atmospheric pressure (0 Pa) and operating temperature (800 • C). Detailed boundary conditions are summarized in Table 2. The electrochemical properties of materials used to predict the air-fuel reaction are organized in Table 3 [25][26][27]29,30]. Electronic and ionic conductivities were considered as a function of temperature. Diffusion parameters for the transport phenomena were estimated at 800 • C and are shown in Table 4 [16,29,31]. Material and thermal properties are summarized in Tables 5 and 6 [29,30,34]. Five simulation cases are presented in Table 7. Case 1 (t cell = 100 µm, T op = 800 • C) was used as a reference case for performance comparison. Other cases with different electrolyte thicknesses and operating temperatures were employed to investigate the effects of these factors on the overall performance.

Results and Discussion
The three-dimensional numerical simulation was conducted to evaluate the effect of parameters on the thermo-fluidic characteristics inside the SOFC unit cell stack. The grid independent test was performed to efficiently obtain accurate results by comparing the performance curve generated by using a finer grid system that has five times more elements in this study. Considering accuracy and computational efficiency, 14,960 elements were employed, as shown in Figure 2a. The cell voltage and the associated current density were first investigated to compare the numerical results with previous literature. As shown in Figure 2b, experimental and numerical data from the literature were obtained to quantitatively validate the accuracy of the numerical model for the present study. The maximum deviation was estimated to be 9% as the current density increased, and the numerical results were similar to those of [29]. The difference may be because the numerical analysis was assumed to be an ideal process for the transport or sealing among multi-layers, unlike the experiments. In addition, the shape parameters for a porous electrode such as tortuosity, porosity, pore diameter, and active specific surface area have various ranges of values depending on the manufacturing process, which determines the internal pore structures even with the same material. Moreover, we chose the shape parameters for the electrode to be constant values obtained from several papers [26,28,30] that were not specified in [29]. So, some numerical errors relative to the experimental data were observed. In general, the increased error can be found in the region of high current density and low voltage. The maximum error was also estimated to be about 10% in the high current density region, and the discrepancy in the high current density region was due to temperature differences caused by ohmic heating in the experiment [29]. Therefore, this does not reflect the accuracy of our results but is a general tendency. There are still other variables determining fuel cell performance, particularly the electrode morphology, which dictates the electrochemical reaction difference at the triple phase boundaries. This is very difficult to model numerically since the morphology changes with increasing temperature and affects the difference between simulation results and literature. Nevertheless, the results of the present study were more consistent with the experimental data than the numerical data found in the literature, especially in the lower current density region. Therefore, the presented numerical modeling may be more reliable for predicting fuel cell performance.
The performances of SOFC unit cells operating under various conditions were evaluated to determine the effect of each parameter. Figure 3 shows the output voltage with current density for all cases. An inverse proportional tendency was observed between the output voltage and current density, which is similar to previous literature [16][17][18][19][20]29,32]. The output voltage decreased with operating temperature because of a reduction in reaction rate and ionic conductivity, which is very sensitive to the ambient temperature. In addition, the difference in output voltage became significant when the current density increased. The maximum relative difference was estimated to be 32.57% based on the results of case 1. To investigate the effect of electrolyte thickness, five total simulation cases were evaluated to compare the performance of SOFCs operated under different conditions as shown in Figure 3. The power density increased as the electrolyte thickness decreased. This is mainly due to the reduction in ohmic loss, which is proportional to the length of ion transport pathway in the electrolyte [18]. Some authors reported that the performance decreased when the electrolyte thickness was near the sub-micron scale [35]. However, it is not essential to consider this criterion compared to the electrolyte thickness scale in this study.
Chemical reactions must be analyzed in terms of the heat and mass transfer to investigate the effect of various parameters. When the air and fuel are injected, chemical reactions simultaneously take place at the interface between the electrolyte and electrodes. Air and fuel channels were designed in a cross-flow configuration in this study. Both oxygen and hydrogen consumption rate showed identical tendencies because the ionic and electronic transfer were coupled in the physical model during the computation. Oxygen and hydrogen distributions along the channel are illustrated at a 0.9 V operating voltage in Figure 4. Case 3 had the lowest operating temperature and showed the highest mass fraction of hydrogen. The profile of hydrogen mass fraction on the anode side is presented in Figure 5. Hydrogen was consumed as the flow passed by, and its reaction rate declined. This decreasing reaction rate along the channel could result in a non-uniform species distribution inside the cell. Moreover, the amount of consumed hydrogen decreased with operating temperature, having Energies 2018, 11, 473 11 of 15 around 33.21% of the maximum difference between case 1 and case 3 at the outlet position. It is clear that the reaction rate profile along the channel length was highly affected by the ambient temperature exerted as energy source determining initial state. On the other hand, the amount of consumed hydrogen slightly increased with decreasing electrolyte thickness. This is because the geometry and material properties resulted in ohmic and concentration losses [36]. Therefore, the SOFC with a thinner electrolyte can provide better performance regarding energy efficiency as well as stacking efficiency for the entire system.
Temperature variation in the unit cell was one of the important factors governing the electrochemical reactions and transport phenomena inside SOFC stacks. In this study, various material properties were considered as a function of temperature to generate accurate results. The three-dimensional temperature contour of case 1 is illustrated in Figure 6. The maximum temperature was estimated to be 802.1 • C in the midst of the cell length at 800 • C. This means that the temperature variation was influenced by the interaction between the transferred heat and the heat of reaction along the channel length. SOFCs are exothermic reaction devices that include oxidation reactions, which can cause temperature differences along the channel length. Temperature profiles along the fuel channel are presented at an operating voltage of 0.9 V, as shown in Figure 7a. The inlet temperature of the fuel channel was fixed at a boundary condition of 800 • C at the zero position. The peak temperature was observed at the normalized channel length of 0.4. An increased temperature scale and variation tendency were reasonable compared to the literature data [16]. Temperature variation concerning the operating temperature was observed for case 1 as shown in Figure 7b. From these results, the maximum temperature difference was estimated to be around 52.11 • C at 0.5 V operating voltage, causing significant thermal stress along the channel. This high-temperature difference occurred under low voltage and was probably due to high current conditions promoting secondary chemical reactions [18]. Moreover, this temperature difference could be controlled by changing the thermal conductivity of interconnects for effective thermal management [24]. The peak temperature increased with decreasing operating voltage, resulting in variations in operating temperature and current densities, hence changing the efficiencies [37]. Moreover, the heat generated by the chemical reaction was much smaller than the order of operating temperature.
The decrease in the operating temperature resulted in a lower reaction rate, which decreased the peak temperature [35]. The thermal expansion of the cell was estimated to be higher in the middle of the channel than other locations due to the maximum temperature difference. Non-uniform ionic/electronic transfer and thermal stress exerted on the stack may cause a loss in SOFC performance and the inability to seal. In addition, differential thermal expansion along the cell can lead to buckling, or even rupture [38]. Thus, fluid-solid interfacial (FSI) simulation involving thermal expansion and fracture should be conducted to design mechanically robust SOFCs. inside the cell. Moreover, the amount of consumed hydrogen decreased with operating temperature, having around 33.21% of the maximum difference between case 1 and case 3 at the outlet position. It is clear that the reaction rate profile along the channel length was highly affected by the ambient temperature exerted as energy source determining initial state. On the other hand, the amount of consumed hydrogen slightly increased with decreasing electrolyte thickness. This is because the geometry and material properties resulted in ohmic and concentration losses [36]. Therefore, the SOFC with a thinner electrolyte can provide better performance regarding energy efficiency as well as stacking efficiency for the entire system. Temperature variation in the unit cell was one of the important factors governing the electrochemical reactions and transport phenomena inside SOFC stacks. In this study, various material properties were considered as a function of temperature to generate accurate results. The three-dimensional temperature contour of case 1 is illustrated in Figure 6. The maximum temperature was estimated to be 802.1 °C in the midst of the cell length at 800 °C . This means that the temperature variation was influenced by the interaction between the transferred heat and the heat of reaction along the channel length. SOFCs are exothermic reaction devices that include oxidation reactions, which can cause temperature differences along the channel length. Temperature profiles along the fuel channel are presented at an operating voltage of 0.9 V, as shown in Figure 7a. The inlet temperature of the fuel channel was fixed at a boundary condition of 800 °C at the zero position. The peak temperature was observed at the normalized channel length of 0.4. An increased temperature scale and variation tendency were reasonable compared to the literature data [16]. Temperature variation concerning the operating temperature was observed for case 1 as shown in Figure 7b. From these results, the maximum temperature difference was estimated to be around 52.11 °C at 0.5 V operating voltage, causing significant thermal stress along the channel. This high-temperature difference occurred under low voltage and was probably due to high current conditions promoting secondary chemical reactions [18]. Moreover, this temperature difference could be controlled by changing the thermal conductivity of interconnects for effective thermal management [24]. The peak temperature increased with decreasing operating voltage, resulting in variations in operating temperature and current densities, hence changing the efficiencies [37]. Moreover, the heat generated by the chemical reaction was much smaller than the order of operating temperature.
The decrease in the operating temperature resulted in a lower reaction rate, which decreased the peak temperature [35]. The thermal expansion of the cell was estimated to be higher in the middle of the channel than other locations due to the maximum temperature difference. Non-uniform ionic/electronic transfer and thermal stress exerted on the stack may cause a loss in SOFC performance and the inability to seal. In addition, differential thermal expansion along the cell can lead to buckling, or even rupture [38]. Thus, fluid-solid interfacial (FSI) simulation involving thermal expansion and fracture should be conducted to design mechanically robust SOFCs.

Conclusions
A numerical simulation of a three-dimensional unit cell stack was conducted to investigate the effect of electrolyte thickness and operating temperature on SOFC performance. We studied unit cell stacks, but the effect from one channel to another was also considered through periodic boundary conditions. We should fully understand how the reactions and heat transfer characteristics inside

Conclusions
A numerical simulation of a three-dimensional unit cell stack was conducted to investigate the effect of electrolyte thickness and operating temperature on SOFC performance. We studied unit cell stacks, but the effect from one channel to another was also considered through periodic boundary conditions. We should fully understand how the reactions and heat transfer characteristics inside

Conclusions
A numerical simulation of a three-dimensional unit cell stack was conducted to investigate the effect of electrolyte thickness and operating temperature on SOFC performance. We studied unit cell stacks, but the effect from one channel to another was also considered through periodic boundary conditions. We should fully understand how the reactions and heat transfer characteristics inside SOFC stacks depend on modeling parameters in the design stage. Therefore, a numerical analysis of a single channel of SOFC stacks having the same geometry from literature was conducted to estimate the variation of mass, momentum, energy transport and electrochemical reactions. Parametric study results showed that an electrolyte as thin as 20% had a higher power density, resulting in a maximum difference of 5% in our simulation case. We also confirmed that the thickness of the electrolyte was highly related to the ohmic loss, which could be used to determine the performance of SOFCs. Moreover, at an operating temperature between 700 • C and 800 • C, the SOFC performance could be substantially enhanced up to 33% due to active chemical reactions by hydrogen utilization. It is essential to evaluate the performance in terms of the overall efficiency of SOFC while considering the operating temperature regarded as input energy. Thus, this study will help predict the interactions between parameters and in designing optimal SOFCs. In future work, the presented numerical model will be integrated with multi-physics modeling to more reliably design a practical SOFC stack.