Dynamic Modeling of a Parallel-Connected Solid Oxide Fuel Cell Stack System

: This study proposes novel simulation methods to model the power delivery function of a parallel-connected solid-oxide-fuel-cell stack system. The proposed methods are then used to investigate the possible thermal runaway induced by the performance mismatch between the employed stacks. A challenge in this modeling study is to achieve the same output voltage but di ﬀ erent output current for each employed stack. Conventional fuel-cell models cannot be used, because they employ fuel ﬂow rates and stack currents as the input variables. These two variables are unknown in the parallel-connected stack systems. The proposed method solves the aforementioned problems by integrating the fuel supply dynamics with the conventional stack models and then arranging them in a multiple-feedback-loop conﬁguration for conducting simulations. The simulation results indicate that the proposed methods can model the transient response of the parallel-connected stack system. Moreover, for the dynamics of the power distribution, there exists an unstable positive feedback loop between employed stacks when the stack temperatures are low, and a stable negative feedback loop when the stack temperatures are high. A thermal runaway could be initiated when the dynamics of the stack temperature is slower than that of the current distribution.


Introduction
Fuel cells (FCs) have become a focus for power generation systems because of their environmentally friendly qualities and high energy-conversion efficiency [1,2]. An FC power plant has to aggregate several FC modules (approximately 10-50 kW), using power electronics such as DC-DC converters, to achieve the kilowatt or megawatt power capacity. Each FC module comprises several FC stacks (approximately 1-5 kW). Each stack contains several hundreds of cells, and each cell outputs a power of 20-30 W [3]. Because thousands of FC cells are within one FC module, each FC stack must have different characteristics due to manufacturing nonuniformity. Furthermore, inside each FC module, these FC stacks are often connected in parallel via mechanical manifolds for receiving the fuel supply and connected in parallel via electrical interconnection for obtaining the electricity output [4,5]. When connecting stacks with different characteristics in parallel, both the fuel supply and load distribution to each stack cannot be even. This nonuniformity problem can lead to severe stack operation and aging problems [6]. Therefore, currently, the uniformity requirement of the FC stack is strict, and the stack is expensive.
The stack resistance variation is one of the major concerns for the solid-oxide-fuel-cell (SOFC) system. The SOFC cell consists of porous electrodes with the thickness of tens of micrometers. Because of the fabrication technology, each cell is different from the other [7], such as electrode thickness and material uniformity. This also makes the initial resistance of each cell different [8]. Under SOFC operation conditions, the cell resistance would increase due to the electrode delamination, formation of that those physic dynamics are not only coupled inside FC stacks but also across FC stacks.
The DC-DC converter is used to isolate the power module from the external load and to maintain a stable voltage at the load. Moreover, the converter prevents the current from flowing back into power modules when power modules with different characteristics are connected together. A conventional DC-DC converter mainly comprises a pulse width modulation (PWM) regulator, capacitor, inductor, switch, diode, and resistor. In this configuration, when the power module does not supply sufficient power to the load, the load voltage is less than its designated value. The feedback signal requests the PWM regulator to increase the duty ratio of the switch control signal, and thus the inductor drains more current out of the power module. In other words, when adapting the DC-DC converter, the power request from the load is converted to the output current request of the power module. Figure 1. Schematic of the two fuel-cell (FC) stacks connected in parallel. The stack system is followed by a DC-DC converter to deliver power to the load. The green arrows represent the fuel supplying and exiting the stacks; the blue arrows represent electricity output; the red arrows represent the heat exchange between stacks.

Mathematics Models of the Parallel-Connected SOFC Stacks
According to the system configuration displayed in Figure 1, the system model should include the dynamics of the electrochemical reactions, thermal fluid, and power electronics (such as DC-DC converter). However, because the dynamics of the DC-DC converter is much faster than the others, the dynamics of the DC-DC converter is neglected in the simulations for simplicity. Figure 1. Schematic of the two fuel-cell (FC) stacks connected in parallel. The stack system is followed by a DC-DC converter to deliver power to the load. The green arrows represent the fuel supplying and exiting the stacks; the blue arrows represent electricity output; the red arrows represent the heat exchange between stacks.
The power generation of the FC system depends on electrochemical reactions inside the FC stacks. These reactions are determined by factors such as the fluid mass flow rates, gas pressures, stack temperature, and power loading. The electrochemical reaction is complex because the aforementioned physics dynamics are all coupled together. Moreover, the reaction becomes more complicated when these stacks are connected in parallel. For example, the fuel distribution in a parallel-connected stack system varies with the channel pressures of each stack. Furthermore, the channel pressure and fuel composition in each stack are affected by electrochemical reactions, which are determined by the power distribution. However, the power distribution cannot be determined without knowing the operation conditions of each stack first. The aforementioned statements indicate that those physic dynamics are not only coupled inside FC stacks but also across FC stacks.
The DC-DC converter is used to isolate the power module from the external load and to maintain a stable voltage at the load. Moreover, the converter prevents the current from flowing back into power modules when power modules with different characteristics are connected together. A conventional DC-DC converter mainly comprises a pulse width modulation (PWM) regulator, capacitor, inductor, switch, diode, and resistor. In this configuration, when the power module does not supply sufficient power to the load, the load voltage is less than its designated value. The feedback signal requests the PWM regulator to increase the duty ratio of the switch control signal, and thus the inductor drains more current out of the power module. In other words, when adapting the DC-DC converter, the power request from the load is converted to the output current request of the power module.

Mathematics Models of the Parallel-Connected SOFC Stacks
According to the system configuration displayed in Figure 1, the system model should include the dynamics of the electrochemical reactions, thermal fluid, and power electronics (such as DC-DC converter). However, because the dynamics of the DC-DC converter is much faster than the others, the dynamics of the DC-DC converter is neglected in the simulations for simplicity.

Fluid Dynamics Modeling
The mass of each gas inside the FC stack can be modeled by the law of conservation of mass. By using the anode channel as an example, the mass of gas i in the FC stack j is described as follows: . n react S j ,an,i = . n out S j ,an,i = y S j ,an,i · . n out S j ,an , j = 1, where . n in S j ,an,i , . n react S j ,an,i , and . n out S j ,an,i are the inlet, reactive, and outlet mass flow rates of gas i, respectively; N 0 represents the number of cells in an FC stack; I S j is the current output of the stack j, x S j ,an,i , and y S j ,an,i are the mole fractions of gas i of the stack j; . n in S j ,an denotes the fuel entering rate; and . n out S j ,an represents the discharging rate. According to [27], the mass flow rate can be modeled by linearizing the flow equation of the nozzle when the pressure difference between upstream and downstream is small and the flow is operating in the subcritical region. Therefore, . n in S j ,an and . n out S j ,an can be modeled as follows: . n in S j ,an = K in S j ,an · (P up S,an − P S j ,an ) . n out S j ,an = K out S j ,an · (P S j ,an − P down S,an ) (6) where K in S j ,an and K out S j ,an are the two constants of the inlet and outlet flow rates; P up S,an and P down S,an are the upstream and downstream pressures of the stack; and P S j ,an is the anode channel pressure of the stack j. Because the inlets (outlets) of the two stacks are connected, the inlet and outlet pressures of the two stacks are the same. Similarly, the mass of the gases in the cathode channel (n S j ,ca,O 2 and n S j ,ca,N 2 ) can be calculated analogously by using Equations (1)- (6).
The gas pressure in the anode channel can be modeled by using the ideal gas law as follows [20]: where R is the gas constant, V an is the volume of the anode channel, and T S j represents the stack temperature. The FC stack temperature can be modeled by using the energy conservation law as follows: . . .
Energies 2020, 13, 501 5 of 20 where A rad is the effective area of the heat radiation, ε is the emissivity of air, and σ denotes the Boltzman constant [28]. Here, m cell and Cv cell are the equivalent mass and specific heat capacity of the FC cell, respectively. Moreover, m int and Cv int are the equivalent mass and specific heat capacity of the interconnector, respectively, and m ste and Cv ste are the equivalent mass and specific heat capacity of the stack fixture. P FC is the electricity power of the FC stack, and T in S j ,an,i and T in S j ,ca,i represent the temperatures of the gases entering the anode channel and the cathode channel, respectively. Cp i is the specific heat capacity of gas, i; ∆Ĥ o r is the reaction heat of the electrochemical reactions; and . Q rad is the thermal radiation between stacks.

Electrochemical Reaction Modeling
In an SOFC stack, oxygen dissociates into oxygen ions when it receives electrons from cathode electrodes. The oxygen ions pass through the electrolyte and react with hydrogen in the anode channel. These reactions generate electrons to source the external load. The ideal cell voltage in the stack j (E cell j ) can be described by using the Nernst equation [19,20].
where −∆g j represents the Gibbs free energy, R is the ideal gas constant, and F is the Faraday constant. Several voltage loss (polarization) mechanisms are present when the FC cell is delivering currents. The output voltage of the FC cell is modeled as follows: where η act j is the activation polarization, η ohm j is the ohmic polarization, η conc j represents the concentration polarization, i cell j is the output current density of a cell, and D j is the ohmic deterioration factor. A a j , A c j , A e , B a , B c , and B e are the ohmic-loss coefficients of the electrodes; δ a , δ c , and δ e are the thicknesses of the electrodes; i 0,a j and i 0,c j represent the respective exchange current density of the anode and cathode for calculating the activation polarization; and i L,a j and i L,c j are the respective limiting current densities of the anode and cathode for calculating the concentration polarizations. The aforementioned current density can be obtained by the following equation [29,30]: where γ a , γ c , E act, a , and E act, c are the activation polarization coefficients, and D e f f , a and D e f f , c represents the effective gaseous diffusivities through the electrode. According to Equations (13)- (21), the cell voltage depends on the gas pressure, output current, and stack temperature. Moreover, both the Nernst voltage and ohmic polarization have negative temperature coefficients. Both the concentration and activation polarizations are complicated functions of the stack temperature.
Equations (1)- (12) can be used to obtain the gas mass, gas pressure, and stack temperature. The values of these parameters are fed into Equations (13)- (21), to calculate the output voltage of a single cell. Finally, the output current (I FC j ), stack voltage (V FC j ), and output power (P FC j ) of the stack j can be calculated as follows: where A cell is the cell area. The stack model displayed above is very similar to the lumped stack model used in previous studies [19,20]. It differs from the lumped stack model in the following aspects: (1) a flow rate model is included to distribute the fuel to each stack, and (2) the heat radiation term, . Q rad , should be included to model the heat exchange between stacks. Note that Equations (13)- (21) indicate that the stack voltage can be calculated when the stack current and chamber pressure are specified beforehand.

Proposed Algorithms for Modeling the Parallel-Connected Stacks
To comply with physics specifications, each power module in a parallel-connected network should have the same output voltage. Moreover, when power electronics are adapted to connect the power modules and the load, the output characteristics (current and voltage) of the power module are uncertain. The aforementioned statements are particularly true for an FC system because the system is not a constant-voltage source. Thus, conventional stack models that use the stack current as input variables to determine the stack voltage must be applied in different manners to model the parallel-connected FC stack system. Figure 2 presents a block diagram of the proposed simulation method that can model the power-delivering function of the FC system displayed in Figure 1. The proposed method arranges two FC stack models in a multiple-feedback-loop system. Each stack model comprises the expressions presented in Equations (13)- (21), which use the stack current as the input to determine the stack voltage and the output power. Each stack voltage and output power values are used to calculate the averaged stack voltage and total power output of the multiple-stack system. Because the response from the power input to the voltage output of the DC-DC converter is much faster than that of the electrochemical reactions, the DC-DC converter is modeled as an energy-conversion factor (C p−p ), as shown in Figure 2.
In a parallel-connected power system, the output current of each module should adhere to two requirements: (1) the current distribution should ensure that the voltage output of each module is the same, and (2) the total output power should be sufficient to maintain the designated load voltage. These two requirements are ensured by two sets of feedback loops that are shown in Figure 2. The first set of loops feedback each stack voltage and the averaged stack voltage. The difference between these Energies 2020, 13, 501 7 of 20 two voltages is processed by the integral operations to determine a portion of the stack current (I v,1 and I v,2 , noted in the plot). The second set of loops feedback the total power delivered. The difference between the delivered power and the designated load power (P ref , in Figure 2) is also processed by using the integral operations to determine another portion of the stack current (, in Figure 2). If the parameters of these three integral operations (K I1 , K I2 , and K I3 ) are designed properly, the two stack voltages should be the same. Moreover, the power delivery would be the same as the designated load power.
Energies 2020, 13, x FOR PEER REVIEW 7 of 20 These two requirements are ensured by two sets of feedback loops that are shown in Figure 2. The first set of loops feedback each stack voltage and the averaged stack voltage. The difference between these two voltages is processed by the integral operations to determine a portion of the stack current ( ,1 and ,2 , noted in the plot). The second set of loops feedback the total power delivered. The difference between the delivered power and the designated load power ( , in Figure 2) is also processed by using the integral operations to determine another portion of the stack current ( , in Figure 2). If the parameters of these three integral operations ( 1 , 2 , and 3 ) are designed properly, the two stack voltages should be the same. Moreover, the power delivery would be the same as the designated load power.

Parameter Design for Integral Operations
The design of the parameters 1 , 2 , and 3 must satisfy the following constraints: (1) every feedback loop of the system shown in Figure 2 should be stable, (2) the outcome of the current distribution should not be influenced by the part of the algorithms in which the same voltage is achieved for the two stacks, and (3) the model displayed in Figure 2 should present the dynamics of the power-delivering function of the FC system. According to the block diagram presented in Figure  2, the total output current of the system is as follows: To ensure that the current distribution is not disturbed by the history of the searching process, 1 = 2 and ,1 (0) = ,2 (0) = 0. Moreover, 1 and 2 must be designed in a manner such that all the feedback loops are stable and the convergence of ( ) is fast. Therefore, ( ( ), 2 ( )) could represent the transient response of the stack voltage and current of the parallel-connected system. In this manner, the output power of the system can be described as follows.

Parameter Design for Integral Operations
The design of the parameters K I1 , K I2 , and K I3 must satisfy the following constraints: (1) every feedback loop of the system shown in Figure 2 should be stable, (2) the outcome of the current distribution should not be influenced by the part of the algorithms in which the same voltage is achieved for the two stacks, and (3) the model displayed in Figure 2 should present the dynamics of the power-delivering function of the FC system. According to the block diagram presented in Figure 2, the total output current of the system is as follows: To ensure that the current distribution is not disturbed by the history of the searching process, K I1 = K I2 and I v,1 (0) = I v,2 (0) = 0. Moreover, K I1 and K I2 must be designed in a manner such that all the feedback loops are stable and the convergence of V avg (t) is fast. Therefore, (V avg (t), 2I P (t)) could represent the transient response of the stack voltage and current of the parallel-connected system. In this manner, the output power of the system can be described as follows.
Energies 2020, 13, 501 8 of 20 After taking the Laplace transformation of the aforementioned equation, the dynamics of the power-delivering function can be approximated as follows: where V avg is the stack voltage at the operation point. According to Equation (27), the value of K I3 is determined by the response time of the FC system. Moreover, the values of K I1 and K I2 should be set based on the stability and bandwidth of the feedback system displayed in Figure 2. However, the FC system is highly nonlinear, and the stability analysis of the system is not the focus of this study. We designed the values of K I1 and K I2 by using empirical methods. Moreover, a suitable selection of the values of these parameters is as follows: K I1 = K I2 = 2 and K I3 = 1.

Simulation Results
In this study, we used the MATLAB/Simulink (MathWorks, Massachusetts, United States) to construct the model and investigate the performance of the FC system. We first validated the proposed models and simulation methods and then checked the performance of the system when the employed stacks have different characteristics. The parameters of the SOFC stack model are listed in Tables 1  and 2. Table 1 lists the parameters of the thermal-fluid dynamics of the stack. All of these parameters, except for the mass of the stack fixture, are excerpted from a previous study [30]. Table 2 lists the parameters of electrochemical reactions. All of these parameters, except for the cell area, A cell , and the number of cells in a stack, N o , are excerpted from a study [30].

Model Validation
Ideally, the proposed model and simulation methods should be verified by using the experimental results of the stack transient response. Unfortunately, we could not find those data in the literature. Instead, we compared our results with the experimental results of the static stack response presented in [31] for different operating temperatures, as shown in Figure 3. The static stack response (I-V curve) was tested under the following conditions: 97% of hydrogen and 3% of steam in the anode channel, 21% of oxygen and 79% of nitrogen in the cathode channel, 101325 Pa of the anode and cathode pressures, and stack temperature ranging from 873 to 1073 K. Experimental results excerpted from the aforementioned study [31] were denoted by solid-marks (triangles, circles, and diamonds); the I-V curve predicted by the static stack model (Equations (13)- (21) and parameters listed in Table 2) was represented by hollow-marks, and predicted by the steady-state value of the dynamic models (Equations (1)- (21) and parameters listed in Tables 1 and 2) were depicted by a red line. According to the plot, the proposed model predicts less polarization loss than the experimental data. The averaged relative-error between the static model and experimental data is 11.6%. Besides, the results from the proposed dynamic model are very close to the static stack model at the temperature of 873 K. Ideally, the proposed model and simulation methods should be verified by using the experimental results of the stack transient response. Unfortunately, we could not find those data in the literature. Instead, we compared our results with the experimental results of the static stack response presented in [31] for different operating temperatures, as shown in Figure 3. The static stack response (I-V curve) was tested under the following conditions: 97% of hydrogen and 3% of steam in the anode channel, 21% of oxygen and 79% of nitrogen in the cathode channel, 101325 Pa of the anode and cathode pressures, and stack temperature ranging from 873 to 1073 K. Experimental results excerpted from the aforementioned study [31] were denoted by solid-marks (triangles, circles, and diamonds); the I-V curve predicted by the static stack model (Equations (13)- (21) and parameters listed in Table 2) was represented by hollow-marks, and predicted by the steady-state value of the dynamic models (Equations (1)- (21) and parameters listed in Tables 1 and 2) were depicted by a red line. According to the plot, the proposed model predicts less polarization loss than the experimental data. The averaged relative-error between the static model and experimental data is 11.6%. Besides, the results from the proposed dynamic model are very close to the static stack model at the temperature of 873 K. Comparisons between the I-V curve data of a SOFC system obtained from three sources: experimental data excerpted from [31], simulation data obtained from the static electrochemical model presented in this study, and simulation data obtained from the dynamic model stated in this study.

Dynamic Response of the Parallel-Connected SOFC Stacks
In the following simulations, the initial temperature of the two FC stacks was 873 K. Moreover, the initial pressures were 1.287 × 10 5 and 1.236 × 10 5 Pa for the anode and cathode channels, respectively. From 5 × 10 2 to 3.5 × 10 3 s, the load power requirement ramped up from 0 to 1.5 × 10 3 W at a rate of 0.5 W • s −1 . From 3.5 × 10 3 to 10 4 s, the load power was maintained at 1500 W. For comparison purposes, the next few figures reveal the performance of each stack under two operating conditions. The responses of stack 1 in case 1, stack 2 in case 1, stack 1 in case 2, and stack 2 in case 2 were represented by blue line, green dotted line, red dotted line, and purple dashed line, respectively.

Same Stacks in the Parallel Connection
In this simulation, stacks with the same performance were connected in parallel for the power delivery. The ohmic deterioration factors, presented in Equation (16), of the two stacks were ( 1 , 2 ) = (1, 1) in case 1 and ( 1 , 2 ) = (3.5, 3.5) in case 2. Figure 4 displays the dynamic response of the output current, voltage, output power, and temperature of each stack. According to the simulation results, Figure 3. Comparisons between the I-V curve data of a SOFC system obtained from three sources: experimental data excerpted from [31], simulation data obtained from the static electrochemical model presented in this study, and simulation data obtained from the dynamic model stated in this study.

Dynamic Response of the Parallel-Connected SOFC Stacks
In the following simulations, the initial temperature of the two FC stacks was 873 K. Moreover, the initial pressures were 1.287 × 10 5 and 1.236 × 10 5 Pa for the anode and cathode channels, respectively. From 5 × 10 2 to 3.5 × 10 3 s, the load power requirement ramped up from 0 to 1.5 × 10 3 W at a rate of 0.5 W · s −1 . From 3.5 × 10 3 to 10 4 s, the load power was maintained at 1500 W. For comparison purposes, the next few figures reveal the performance of each stack under two operating conditions. The responses of stack 1 in case 1, stack 2 in case 1, stack 1 in case 2, and stack 2 in case 2 were represented by blue line, green dotted line, red dotted line, and purple dashed line, respectively.

Same Stacks in the Parallel Connection
In this simulation, stacks with the same performance were connected in parallel for the power delivery. The ohmic deterioration factors, presented in Equation (16), of the two stacks were (D 1 , D 2 ) = (1, 1) in case 1 and (D 1 , D 2 ) = (3.5, 3.5) in case 2. Figure 4 displays the dynamic response of the output current, voltage, output power, and temperature of each stack. According to the simulation results, both stacks had the same response because they were exactly the same. At the end of the simulation time, both stacks exhibited responses of 62.7 A, 12 V, 750 W, and 1305 K in case 1 and 67.6 A, 11.1 V, 750 W, and 1392 K in case 2. Both systems output power of 1500 W, as required. The dynamics of the stack temperature was slower than that of power delivery. Moreover, an FC system employing small ohmic-loss stacks delivered smaller currents, provided higher voltages, and had a lower stack temperature compared with a system employing large ohmic-loss stacks.
Energies 2020, 13, x FOR PEER REVIEW 10 of 20 both stacks had the same response because they were exactly the same. At the end of the simulation time, both stacks exhibited responses of 62.7 A, 12 V, 750 W, and 1305 K in case 1 and 67.6 A, 11.1 V, 750 W, and 1392 K in case 2. Both systems output power of 1500 W, as required. The dynamics of the stack temperature was slower than that of power delivery. Moreover, an FC system employing small ohmic-loss stacks delivered smaller currents, provided higher voltages, and had a lower stack temperature compared with a system employing large ohmic-loss stacks.  Figure 5 displays the response of the stack voltage composition, which includes Nernst voltage, activation polarization, ohmic polarization, and concentration polarization. According to simulation results, the Nernst voltage decreased as time increased, because of its negative temperature coefficient. Both activation and ohmic polarizations increased rapidly at the beginning and then decreased slowly. This phenomenon occurs because the output current increased rapidly at the beginning, and then the stack temperature increased slowly and dominated the response. The concentration polarization increased at the beginning, due to the increase in the output current, and kept increasing due to its positive temperature coefficient. Figure 6 presents the anode channel pressure, cathode channel pressure, hydrogen utilization, and oxygen utilization in this example. According to simulation results, the anode and cathode pressures were 1.287 × 10 5 and 1.235 × 10 5 Pa , respectively, in both cases. The hydrogen utilization was 61.4% in case 1 and 66.2% in case 2, whereas the oxygen utilization was 47.7% in case 1 and 51.1% in case 2. These simulation results indicated that the channel pressure did not vary significantly when the current (electrochemical reactions) in two cases varied from 67.6 to 62.7 A.  Figure 5 displays the response of the stack voltage composition, which includes Nernst voltage, activation polarization, ohmic polarization, and concentration polarization. According to simulation results, the Nernst voltage decreased as time increased, because of its negative temperature coefficient. Both activation and ohmic polarizations increased rapidly at the beginning and then decreased slowly. This phenomenon occurs because the output current increased rapidly at the beginning, and then the stack temperature increased slowly and dominated the response. The concentration polarization increased at the beginning, due to the increase in the output current, and kept increasing due to its positive temperature coefficient. Figure 6 presents the anode channel pressure, cathode channel pressure, hydrogen utilization, and oxygen utilization in this example. According to simulation results, the anode and cathode pressures were 1.287 × 10 5 and 1.235 × 10 5 Pa, respectively, in both cases. The hydrogen utilization was 61.4% in case 1 and 66.2% in case 2, whereas the oxygen utilization was 47.7% in case 1 and 51.1% in case 2. These simulation results indicated that the channel pressure did not vary significantly when the current (electrochemical reactions) in two cases varied from 67.6 to 62.7 A.

Different Stacks in the Parallel Connection
In this simulation, stacks that have different ohmic properties were connected in parallel for the power delivery function. The ohmic deterioration factors of the two stacks were (D 1 , D 2 ) = (1, 3.5) in both cases. Case 1 includes the heat exchange between two stacks . Q rad , and case 2 excluded . Q rad . Figure 7 shows the dynamic response of the stack current, voltage, power, and temperature of each stack. In case 1, the output current of stack 1 (small ohmic loss) was much larger than that of stack 2 (large ohmic loss) during the power ramp up. The current of stack 2 slowly became equivalent to stack 1 as time increased. At 10 4 s, the currents of stacks 1 and 2 were 67 and 63.2 A, which correspond to the power output of 771.9 and 728.1 W, respectively. Both stacks had the same voltage even when they had different currents. The temperature responses of the two stacks were almost the same due to the heat radiation. In case 2, stack 1 still retained a larger current output than stack 2. Both currents increased and the simulation terminated when the current of stack 1 attained a value of 98.4 A, which was the current limit (i L,a j in Equation (20)) of the concentration polarization. The temperature responses of the two stacks differed largely from each other because there was no heat exchange between them. According to these simulation results, the parallel-connected stacks that had different ohmic properties encountered a thermal runaway when the heat exchange between stacks was prohibited. Moreover, compared with the temperature results presented in the previous simulation, the output current in this case attained its current limit, which was not resulted from the stack temperature. This suggests that the temperature control of the power module may not be useful to prevent this thermal runaway from occurring.

Different Stacks in the Parallel Connection
In this simulation, stacks that have different ohmic properties were connected in parallel for the power delivery function. The ohmic deterioration factors of the two stacks were ( 1 , 2 ) = (1, 3.5) in both cases. Case 1 includes the heat exchange between two stacks ̇, and case 2 excluded ̇. Figure 7 shows the dynamic response of the stack current, voltage, power, and temperature of each stack. In case 1, the output current of stack 1 (small ohmic loss) was much larger than that of stack 2 (large ohmic loss) during the power ramp up. The current of stack 2 slowly became equivalent to stack 1 as time increased. At 10 4 s, the currents of stacks 1 and 2 were 67 and 63.2 A, which correspond to the power output of 771.9 and 728.1 W, respectively. Both stacks had the same voltage even when they had different currents. The temperature responses of the two stacks were almost the same due to the heat radiation. In case 2, stack 1 still retained a larger current output than stack 2. Both currents increased and the simulation terminated when the current of stack 1 attained a value of 98.4 A, which was the current limit ( , in Equation (20)) of the concentration polarization. The temperature responses of the two stacks differed largely from each other because there was no heat exchange between them. According to these simulation results, the parallel-connected stacks that had different ohmic properties encountered a thermal runaway when the heat exchange between stacks was prohibited. Moreover, compared with the temperature results presented in the previous simulation, the output current in this case attained its current limit, which was not resulted from the stack temperature. This suggests that the temperature control of the power module may not be useful to prevent this thermal runaway from occurring.   13 of 20 different current distributions. In case 2, the Nernst voltages of the two stacks differed largely due to the large temperature difference between the two stacks. Moreover, the concentration polarization increased sharply because the output current reached its maximum current limit. This led to the zero stack voltage shown in Figure 7. increased sharply because the output current reached its maximum current limit. This led to the zero stack voltage shown in Figure 7. Figure 9 displays the channel pressure and fuel utilization in this example. According to the simulation results, in both cases, the anode channel pressures of the two stacks were stable at 1.287 × 10 5 Pa, while the cathode channel pressures varied from 1.241 × 10 5 to 1.231 × 10 5 Pa. The pressure variation was less than 0.8% when the current output varied from 0 to 98.4 A. In case 1, the hydrogen utilization rates of the two stacks were 65.6% and 61.6%, and the oxygen utilization rates were 50.7% and 48.0%. In case 2, the hydrogen utilization rates of stack 1 reached 100%. Therefore, it can be concluded that the channel pressure and, thus, the inlet fuel distribution are not the major reasons that caused the instability.   Figure 9 displays the channel pressure and fuel utilization in this example. According to the simulation results, in both cases, the anode channel pressures of the two stacks were stable at 1.287 × 10 5 Pa, while the cathode channel pressures varied from 1.241 × 10 5 to 1.231 × 10 5 Pa. The pressure variation was less than 0.8% when the current output varied from 0 to 98.4 A. In case 1, the hydrogen utilization rates of the two stacks were 65.6% and 61.6%, and the oxygen utilization rates were 50.7% and 48.0%. In case 2, the hydrogen utilization rates of stack 1 reached 100%. Therefore, it can be concluded that the channel pressure and, thus, the inlet fuel distribution are not the major reasons that caused the instability.

Stacks Have Less Heat Capacity
In this simulation, the stack parameters were the same as those used in Section 5.2.2, except that the equivalent fixture masses ( in Equation (8)) of the two stacks were reduced to 0.01 • . According to the results presented in Figure 10, in case 1 (allows heat exchange between stacks), stack 1 had an output of 66.7 A, 11.5 V, 767.1 W, and 1368 K, and stack 2 had an output of 63.5 A, 11.5 V, 730.3 W, and 1368 K. The stack with the small ohmic loss had a larger output power than the stack with a large ohmic loss. In case 2 (without heat exchange between stacks), the transient responses of stacks were similar to those in case 1. However, the performance difference between stacks in case 2 was larger than that in case 1. Moreover, stack 2 had almost the same steady-state performance as that in stack 1, although the employed stacks are different, and no heat was exchanged between them.

Stacks Have Less Heat Capacity
In this simulation, the stack parameters were the same as those used in Section 5.2.2, except that the equivalent fixture masses (m ste in Equation (8)) of the two stacks were reduced to 0.01 m ste . According to the results presented in Figure 10, in case 1 (allows heat exchange between stacks), stack 1 had an output of 66.7 A. According to the results presented in Figure 10, in case 1 (allows heat exchange between stacks), stack 1 had an output of 63.5 A, 11.5 V, 767.1 W, and 1368 K, and stack Energies 2020, 13, x FOR PEER REVIEW 15 of 20

Discussion
This study proposed a novel method to model the current distribution of parallel-connected SOFC stacks. According to the simulation results presented in Figures 4, 7, and 10, each parallelconnected stack has the same stack voltage, even when its output current differs from the others or is unbounded in some cases. Therefore, the proposed method is effective to solve this "same voltage but different current" problem. Moreover, the proposed method is better than the existing approaches that employ numerical searches to find the current distribution due to its lower amount of effort on the algorithm implementation. For example, a 50 kW SOFC power module normally comprises eight stacks in parallel [4]. To find the current distribution by conducting a numerical search, one has to derive an 8 × 8 Jacobian matrix for those highly nonlinear dynamics presented in Equations (13)- (21) and may also have to perform the matrix inversion. By contrast, the proposed method only requires eight feedback loops, which largely reduce the simulation effort required.
As aforementioned, the power-delivering behavior of the parallel-connected FC stacks is complicated because of the coupled physics inside the stacks and across stacks. Those coupled physics can be classified as fluid, electric, and thermal dynamics. According to the simulation results illustrated in Figures 5 and 8, the channel pressures varied by less than 0.8% when the output current varied from 0 to 98.4 A. The channel pressures do not vary much, because the same amount of steam is generated when the hydrogen is consumed in the anode channel. Moreover, oxygen only accounts for approximately 10% of the total volume in the cathode channel. After excluding the effect of channel pressure (fluid dynamics), the power distribution is dominated by the coupled dynamics between electricity and temperature.

Discussion
This study proposed a novel method to model the current distribution of parallel-connected SOFC stacks. According to the simulation results presented in Figure 4, Figure 7, and Figure 10, each parallel-connected stack has the same stack voltage, even when its output current differs from the others or is unbounded in some cases. Therefore, the proposed method is effective to solve this "same voltage but different current" problem. Moreover, the proposed method is better than the existing approaches that employ numerical searches to find the current distribution due to its lower amount of effort on the algorithm implementation. For example, a 50 kW SOFC power module normally comprises eight stacks in parallel [4]. To find the current distribution by conducting a numerical search, one has to derive an 8 × 8 Jacobian matrix for those highly nonlinear dynamics presented in Equations (13)- (21) and may also have to perform the matrix inversion. By contrast, the proposed method only requires eight feedback loops, which largely reduce the simulation effort required.
As aforementioned, the power-delivering behavior of the parallel-connected FC stacks is complicated because of the coupled physics inside the stacks and across stacks. Those coupled physics can be classified as fluid, electric, and thermal dynamics. According to the simulation results illustrated in Figures 5 and 8, the channel pressures varied by less than 0.8% when the output current varied from 0 to 98.4 A. The channel pressures do not vary much, because the same amount of steam is generated when the hydrogen is consumed in the anode channel. Moreover, oxygen only accounts for approximately 10% of the total volume in the cathode channel. After excluding the effect of channel pressure (fluid dynamics), the power distribution is dominated by the coupled dynamics between electricity and temperature. Case 1 in Section 5.2.2 presents the system response when the parallel-connected stacks have different ohmic properties. Because the two stack temperatures are almost the same (Figure 4d), the temperature effect can be excluded from the power (current) distribution. Moreover, if we only focus on the current difference in the two stacks, the portion of the current that is responsible for the load power delivery (I p in Figure 2) can be ignored. Consequently, the current distribution is determined by the current that makes the two stack voltages the same (I v,1 and I v,2 in Figure 2). In that case, I v,1 is larger than I v,2 at the beginning because stack 1 has a smaller ohmic loss than stack 2. Due to the increase in the power demand, I p increases. Thus, the difference between I v,1 and I v,2 has to be increased to compensate for the enlarged difference of the ohmic polarization. The stack temperature increases with the stack currents and then dominates all the polarization losses. Thus, the polarization differences between stacks and current distribution decreases. These mechanisms explain why the current is almost equally distributed in the steady state, albeit the employed stacks have different ohmic properties.
When the temperature effect is included, the power-sharing function in the parallel-connected system is very complicated. Figure 11 shows the plot of the cell voltage versus temperature by using Equations (13)- (21), when the pressure effect is excluded from the calculations. According to the plot, the stack voltage decreases as the stack current increases. Moreover, the stack voltages have positive temperature coefficients in the low-temperature and large-current region. The lower the temperature and the larger the current are, the more positive the temperature coefficients are. When the stack temperature increases, the stack voltages are less affected by the current, and their temperature coefficients decrease gradually and become negative.
Energies 2020, 13, x FOR PEER REVIEW 16 of 20 Case 1 in Section 5.2.2 presents the system response when the parallel-connected stacks have different ohmic properties. Because the two stack temperatures are almost the same (Figure 4d), the temperature effect can be excluded from the power (current) distribution. Moreover, if we only focus on the current difference in the two stacks, the portion of the current that is responsible for the load power delivery ( in Figure 2) can be ignored. Consequently, the current distribution is determined by the current that makes the two stack voltages the same ( ,1 and ,2 in Figure 2). In that case, ,1 is larger than ,2 at the beginning because stack 1 has a smaller ohmic loss than stack 2. Due to the increase in the power demand, increases. Thus, the difference between ,1 and ,2 has to be increased to compensate for the enlarged difference of the ohmic polarization. The stack temperature increases with the stack currents and then dominates all the polarization losses. Thus, the polarization differences between stacks and current distribution decreases. These mechanisms explain why the current is almost equally distributed in the steady state, albeit the employed stacks have different ohmic properties.
When the temperature effect is included, the power-sharing function in the parallel-connected system is very complicated. Figure 11 shows the plot of the cell voltage versus temperature by using Equations (13)-(21), when the pressure effect is excluded from the calculations. According to the plot, the stack voltage decreases as the stack current increases. Moreover, the stack voltages have positive temperature coefficients in the low-temperature and large-current region. The lower the temperature and the larger the current are, the more positive the temperature coefficients are. When the stack temperature increases, the stack voltages are less affected by the current, and their temperature coefficients decrease gradually and become negative. Because the electrochemical reaction is exothermic, the stack that outputs a higher amount of current would have a higher stack temperature. In the parallel-connected stack system, the stacks that have the more positive temperature coefficients have to output a higher amount of current than the others to establish the same stack voltage. In the example presented in case 2 in Section 5.2.2, the initial current of stack 1 is larger than that of stack 2 due to the smaller ohmic loss of stack 1. The temperature coefficient of the stack 1 (large current) is more positive than that of stack 2 (small current) when both stack temperatures are low. Consequently, these two parallel-connected stacks enter a positive feedback loop, and the current difference between stacks increases sharply. Simultaneously, both stack currents increase due to the load demand. Before the temperatures of both stacks increase, which would minimize the polarization difference between stacks, the current of stack 1 attains its maximum current limit. These mechanisms explain the thermal runaway induced by the mismatch between the stack characteristics. In case 2, presented in Section 5.2.3, the temperature response is rapid. The increase in the temperature of stack 1 decreases the temperature coefficient of its stack voltage. Thus, the value of the coefficient is close to that of stack 2. The unstable current distribution may either exist only for a short duration or may not exist at all. Since both stack Because the electrochemical reaction is exothermic, the stack that outputs a higher amount of current would have a higher stack temperature. In the parallel-connected stack system, the stacks that have the more positive temperature coefficients have to output a higher amount of current than the others to establish the same stack voltage. In the example presented in case 2 in Section 5.2.2, the initial current of stack 1 is larger than that of stack 2 due to the smaller ohmic loss of stack 1. The temperature coefficient of the stack 1 (large current) is more positive than that of stack 2 (small current) when both stack temperatures are low. Consequently, these two parallel-connected stacks enter a positive feedback loop, and the current difference between stacks increases sharply. Simultaneously, both stack currents increase due to the load demand. Before the temperatures of both stacks increase, which would minimize the polarization difference between stacks, the current of stack 1 attains its maximum current limit. These mechanisms explain the thermal runaway induced by the mismatch between the stack characteristics. In case 2, presented in Section 5.2.3, the temperature response is rapid. The increase in the temperature of stack 1 decreases the temperature coefficient of its stack voltage. Thus, the value of the coefficient is close to that of stack 2. The unstable current distribution may either exist only for a short duration or may not exist at all. Since both stack currents increase stably, the stack temperatures increase and the polarization differences between stacks decrease. Consequently, both the stack currents and stack temperatures are similar to each other, albeit there is no heat exchange between stacks.
According to the aforementioned analysis, the mismatches of the temperature coefficients and stack currents establish an unstable positive feedback loop across the stacks in the low-temperature region. Moreover, there is a stable negative feedback loop for minimizing the performance difference of the employed stacks when the stack temperature is high. Due to the interactions between the two aforementioned feedback loops, thermal runaway may occur in the following cases: (1) a large power requirement leads to a large-current output and higher stack temperatures, which in turns lower the maximum current limit, and (2) slow temperature dynamics compares to the current distribution dynamics, which enlarges the current differences between stacks. Thus, one of the employed stacks attain its maximum current limit. According to the system configuration shown in Figure 1, load power shortage always leads to an increase in the stack current. This is because of the associated power electronics and the associated control methods presented in Figure 1. However, increase in the power output of the FC systems can be realized by either increasing or decreasing the output current, depending on the operation point. The proposed analysis suggests that it is beneficial to develop power electronics exclusively for FC systems. Thus, the employed FC stacks may have larger performance differences and the parallel-connected system can operate closer to its maximum power operation point. Further study of this problem will be conducted in the future.

Conclusions
This study proposed a novel method to model the power delivering operation of a parallel-connected SOFC system. This model comprises a parallel-connected fuel supply architecture, parallel-connected electrical output, and heat transfer between stacks. One of the challenges in this task is to obtain the same voltage but different current output for each employed stack. The proposed method solves this problem by using multiple feedback loops with the integral actions in the computation algorithms. The simulation results of this study indicated that when the proposed methods are used, the parallel-connected stacks have the same voltage output when each stack current is different or unbounded.
The proposed model and simulation methods were then used to investigate the possible thermal runaway induced by the mismatch between the stack performance. In these simulations, the required load power ramped up from 0 to 1500 W, at a rate of 0.5 W · s −1 . The simulation results indicate that, when the two stacks are the same, even they have degenerated (the ohmic polarization increases 3.5 times), the parallel-connected stack system can deliver the power that is requested. However, when the employed stacks are different (one stack's ohmic polarization is 3.5 times higher than that of the other), the parallel-connected stack system may become unstable when the temperature response is slow and the heat transfer between stacks is insufficient.
This study indicates that the mismatches between the temperature coefficients and stack currents establish an unstable positive feedback loop when the stack temperatures are low. Moreover, a stable negative feedback loop is present to minimize the performance difference between the employed stacks when the stack temperatures are high. Due to the interactions between these two feedback loops, thermal runaway may occur in the following cases: (1) a large power requirement leads to a large-current output and higher stack temperatures, which in turns lower the maximum current limit, and (2) slow temperature dynamics compares to the current distribution dynamics, which enlarge the current differences between stacks. Thus, one of the employed stacks reaches its maximum current limit.
Author Contributions: Methodology preparation, software usage, data curation, and writing of the original draft was done by C.-C.W. Investigation, formal analysis, supervision, and review and editing of the draft was performed by T.-L.C. All authors have read and agreed to the published version of the manuscript.
Funding: This research received no external funding.
Acknowledgments: This work was supported by the Ministry of Science and Technology (MOST 108-2221-E-009-108).

Conflicts of Interest:
The authors declare no conflicts of interest. constant pressure specific heat, J · kg −1 · K −1 Cv constant volume specific heat, J · kg −1 · K −1 σ Boltzman constant, m 2 · kg · s −2 · K −1 ε emissivity of air A area, m 2 ∆Ĥ o r reaction heat of the electrochemical reactions, J · mol −1 i 0 exchange current density, A · m −2 i L limiting current densities, A · m −2 δ a , δ c , δ e thicknesses of the electrodes, m A a j , A c j , A e ohmic-loss coefficients of the electrodes, Ω · m B a , B c , B e ohmic-loss coefficients of the electrodes, K γ a , γ c activation polarization coefficients, A · m −2 E act, a , E act, c activation polarization coefficients, J · mol −1 D e f f , a , D e f f , c effective diffusivities, m 2 · s −1 Superscript in inlet out outlet react reaction up upstream down downstream