## 2 Modeling of Fuel Cell Systems
### Fuel Cell Single Cell Modeling
#### 2.1.2 Bipolar Plate Temperature Model  
One end of the bipolar plate is in contact with the end plate layer, and the heat flow is by heat conduction [21]. The other end is in contact with the MEA, because the contact part of the bipolar plate has a negative/anode gas flow channel, so the heat flow mode has both heat conduction and heat convection. Some fuel cells are designed with extra channels in the plate, which are cooled by cooling water. Bipolar plates have high electrical conductivity and produce a small amount of joule heat during operation. In addition, the enthalpy changes of gas and liquid in and out of the bipolar plate will also affect the plate temperature. The energy conservation equation of the bipolar plate is:  
\[ \begin{array}{c}\left[c_{\rm p,avg}(m_{\rm gas}+m_{\rm liq})+(\rho_{\rm L2} A_{\rm L2}h_{\rm L2}c_{\rm p,L2})\right]\frac{\rm dT_{\rm L2}}{\rm dt}\\ =Q_{\rm L1}+Q_{\rm int1}+Q_{\rm L3}+Q_{\rm L3R}+\Delta Q_{\rm liq}+\Delta Q_{ \rm gas}+\Delta H_{\rm gas}+\Delta H_{\rm liq}\end{array}  \tag{4}\]  
In Formula (4), \(c_{p,arg}\) is the mixing specific heat capacity of gas and liquid; \(m_{gas}\) and \(m_{liq}\) are the flow rates of gas and liquid, respectively. \(Q_{\rm L1}\) is the heat transfer from the end plate layer to the bipolar plate; \(Q_{\rm int1}\) is the joule heat generated by resistance; \(Q_{\rm L3}\) and \(Q_{\rm L3R}\) are the heat transfer and heat convection from the MEA to bipolar plate, respectively; \(\Delta Q_{\rm liq}\) and \(\Delta Q_{\rm gas}\) are, respectively, the heat change of cooling water and gas flowing through the electric reactor. \(\Delta H_{\rm gas}\) and \(\Delta H_{\rm liq}\) are the enthalpy changes of gas and the liquid in the reactor. The calculation formula of \(Q_{\rm L1}\) is similar to that of Formula (3). The calculation formula of \(Q_{\rm L3}\), \(Q_{\rm int1}\) and \(Q_{\rm L3R}\) is:  
\[Q_{\rm L3}=U_{\rm L3}A^{\prime}(T_{\rm L3}-T_{\rm L2}) \tag{5}\]  
\[Q_{\rm int1}=(iA)^{2}\frac{\rho_{\rm R,L2}h_{\rm L2}}{A} \tag{6}\]  
\[Q_{\rm L3R}=U_{\rm L3R}A^{\prime}_{\rm R}(T_{\rm L3}-T_{\rm L2}) \tag{7}\]  
In Formulas (5)-(7), \(U_{\rm L3}\) and \(U_{\rm L3R}\) are the equivalent heat transfer and heat convection coefficient from the MEA layer to the plate layer, respectively. \(T_{\rm L3}\) is the temperature of theMEA layer; I is current density; \(\rho_{\rm{R,L2}}\) is the resistivity of the plate layer; \(W_{\rm{cl}}\) is the cooling water flow rate; and \(c_{\rm{p,cl}}\) is the specific heat capacity of water. The calculation method of heat change caused by water and gas entering and leaving the reactor is similar. Here, the heat change of gas is taken as an example:  
\[Q_{\rm{gas}}=L_{\rm{P,gas}}P_{\rm{Cs,gas}}h_{\rm{gas}}\big{(}T_{\rm{gas.in}}-T _{\rm{gas.out}}\big{)} \tag{8}\]  
In Formula (8), \(L_{\rm{P,gas}}\) is the length of gas passage; \(P_{\rm{Cs,gas}}\) is the cross-section circumference of the gas passage; \(h_{\rm{gas}}\) is the gas convection coefficient; and \(T_{\rm{gas.in}}\) and \(T_{\rm{gas.out}}\) are the temperatures at which the gas enters and leaves the reactor, respectively. In the above equation, it is necessary to calculate the Reynolds number and Nusselt coefficient of the fluid, to obtain the convection coefficient of the gas. The Reynolds number is calculated as follows:  
\[D_{\rm{h,gas}}=\frac{4A_{\rm{gas}}}{P_{\rm{cs,gas}}} \tag{9}\]  
\[Re_{\rm{gas}}=\frac{\rho_{\rm{gas}}Q_{\rm{c,gas}}D_{\rm{h,gas}}}{\mu_{\rm{gas }}A_{\rm{gas}}} \tag{10}\]  
In Equations (6) and (7), \(D_{\rm{h,gas}}\) is the hydraulic diameter of the flow passage; \(P_{\rm{cs,gas}}\) is the cross-sectional area of the flow passage; \(\rho_{\rm{gas}}\) is the gas density; \(Q_{\rm{c,gas}}\) is the gas flow rate; and \(\mu_{\rm{gas}}\) is the dynamic viscosity of gas. In the simulation model, the air gas flow rate is 1-15 L/min, and the hydrogen gas flow rate is 1-5 L/min.