**3. Model Description**

The purpose of the model developed in this paper was firstly to compare the performance and efficiency of various configurations and then to develop a control strategy in subsequent research. To check the system response according to the control design, the model has to be dynamic and not too complex [1,23]. In addition, the needs of dynamic models have been addressed by several authors because of the high operating temperature of SOFCs [24,25]. Therefore, a dynamic and lumped component model was developed in the SIMULINK environment. First, the steady-state characteristics of each system were examined to check the effectiveness of the systems. Dynamic characteristics and control strategies will be analyzed in future papers. In the following sections, the descriptions of each component model are discussed.

#### *3.1. A SOFC Stack*

In this model, a 1 kW class DIR stack has been adopted for an efficient SOFC system, and it is a planar SOFC stack with an anode-supported cell. The conditions inside a SOFC stack are appropriate for the SMR reaction because of its proper operating temperature and electrochemically generated steam. Additionally, internal reforming has several advantages, such as reducing the size of the external reformer and air blower power consumption by means of decreasing the stack air flow rate for stack cooling [23,26,27].

The SOFC stack is composed of reactant channels and the positive-electrolyte-negative (PEN) structure. The model accounts for mass balances, thermal balances and electrochemical reactions. Specifications of the stack are presented in Table 2. Fuel and air mixtures are assumed to follow the ideal gas law.

**Table 2.** A SOFC stack specifications.


#### 3.1.1. Mass Balance Model

In the mass transfer model, the species in the anode channel are considered as CH4, H2O, CO, H2 and CO2, and those in the cathode channel are only N2 and O2. Four reactions inside a stack are considered. All the reactions are presented in Table 3.

**Table 3.** Reactions considered inside a stack.


Mass balance dynamics for *i* species in fuel and air channels are given in Equations (1) and (2). In the fuel channel, (R1)—(R3) reactions described in Table 3 take place. (R4) reaction occurs in the air channel. *ri* is the molar rate of formation for *i* species, and is obtained from the equations below.

$$\frac{\partial \mathbb{C}\_{i,f}}{\partial t} = -\mu\_f \frac{\partial \mathbb{C}\_{i,f}}{\partial x} + r\_{i,f} \quad \text{(i} = \text{CH}\_4, \text{ H}\_2\text{O, CO, } \text{H}\_2, \text{ CO}\_2\text{)}\tag{1}$$

$$\frac{\partial \mathcal{C}\_{i,a}}{\partial t} = -u\_a \frac{\partial \mathcal{C}\_{i,a}}{\partial x} + r\_{i,a} \quad (i = N\_{2\prime} \ O\_2) \tag{2}$$

For calculating the reaction rates of (R1) and (R2), there are various relations accounting for the reaction rates of SMR and WGS [1,20,28]. Among the equations, Chinda et al. [29] model presented in Equations (3)–(8) has been used. Chinda et al. derived the reaction rate from Arhenius' curve fits using the data by Lehnert et al. [30]. *Rk* and *kk* represent the reaction rate and forward reaction rate constant of reaction *k*, respectively. Chinda et al. considered that the SMR reaction occurs at the surface of the anode and the WGS reaction takes place inside the void volume of the anode.

$$R\_{\rm R1} = k\_{\rm R1} \left( P\_{CH\_4} P\_{H\_2O} - \frac{P\_{H\_2}^3 P\_{CO}}{K\_{\rm R1}} \right) \left[ \text{mol } \text{m}^{-2} \text{ s}^{-1} \right] \tag{3}$$

$$R\_{\rm R2} = k\_{\rm R2} \left( P\_{H\_2O} P\_{CO} - \frac{P\_{H\_2} P\_{CO\_2}}{K\_{\rm R2}} \right) \left[ \text{mol} \,\text{m}^{-3} \,\text{s}^{-1} \right] \tag{4}$$

$$k\_{\rm R1} = 2395 \exp\left(-\frac{231266}{RT}\right) \left[\text{mol m}^{-2} \text{ Pa}^{-2} \text{ s}^{-1}\right] \tag{5}$$

$$k\_{\rm R2} = 0.0171 \exp\left(-\frac{103191}{RT}\right) \left[\text{mol m}^{-3} \text{ Pa}^{-2} \text{s}^{-1}\right] \tag{6}$$

*Kk* indicates the equilibrium constant of reaction *k* and can be calculated by Equations (7) and (8) [2], *Z* is defined as (1000/T(K)-1), The unit of *K*R1 is Pa2, and *K*R2 is a dimensionless constant.

$$K\_{\rm R1} = 1.0267 \times 10^{10} \exp\left(-0.2513 Z^4 + 0.3665 Z^3 + 0.5810 Z^2 - 27.13 Z + 3.2770\right) \tag{7}$$

$$K\_{\rm R2} = \exp\left(-0.2935Z^3 + 0.6351Z^2 + 4.1788Z + 0.3169\right) \tag{8}$$

The reaction rates of (R3) and (R4) are equal to the electrochemical reaction rate of the fuel cell and can be obtained by Faraday's law as shown in Equation (9).

$$R\_{\rm R3} = R\_{\rm R4} = R\_{\rm R5} = \frac{J}{2F} \left[ \text{mol s}^{-1} \right] \tag{9}$$

From Equations (3)–(9), the molar rate of formation per volume (*ri*) can be inferred. Nitrogen does not react, so *rN*<sup>2</sup> is zero.

$$r\_{CH\_4} = -\frac{1}{t\_f} R\_{R1} \tag{10}$$

$$r\_{H\_2O} = -\frac{1}{t\_f} R\_{\mathbb{R}1} - \varepsilon\_a R\_{\mathbb{R}2} + \frac{1}{V\_a} R\_{\mathbb{R}3} \tag{11}$$

$$r\_{CO} = \frac{1}{t\_f} R\_{\mathbb{R}1} - \varepsilon\_a R\_{\mathbb{R}2} \tag{12}$$

$$r\_{H\_2} = 3\frac{1}{t\_f}R\_{\rm R1} + \varepsilon\_a R\_{\rm R2} - \frac{1}{V\_a}R\_{\rm R3} \tag{13}$$

$$r\_{\text{CO}\_2} = \varepsilon\_a R\_{\text{R2}} \tag{14}$$

$$
\sigma\_{O\_2} = -0.5 \frac{1}{V\_c} R\_{R3} \tag{15}
$$

The pressure decrease (Δ*p*) through the anode and cathode channels is calculated using Equation (16). Equations (17) and (18) represent the friction factor (*f*) for laminar and turbulent flow, respectively. *ϕ* is the aspect ratio of the channel, and *vg* is the gas velocity.

$$f = \frac{24}{\text{Re}[1 - 1.3553\rho + 1.9467\rho^2 - 1.7012\rho^3 + 0.9564\rho^4 - 0.2537\rho^5]} \text{ (Re} < 2000) \tag{16}$$

$$f = \frac{0.0791}{Re^{0.25}} \ (Re \ge 2000) \tag{17}$$

$$
\Delta p = f \frac{4L\_{ch}}{D\_h} \frac{1}{2} \rho\_\% u\_\%^2 \tag{18}
$$

3.1.2. Electrochemical Reaction Model

Cell potential can be computed by Equation (19). The thermodynamic reversible potential (*Vrev*) is determined by the Nernst equation in Equation (20) [31]. A voltage drop from the thermodynamic reversible potential exists because of overpotentials summarized in Table 4.

$$V\_{ccll} = V\_{rev} - \eta\_{act,a} - \eta\_{act,c} - \eta\_{ohm} - \eta\_{conc} \tag{19}$$

$$V\_{rev} = E\_0 - \frac{RT\_{PEN}}{2F} \ln \frac{P\_{H\_2O}}{P\_{H\_2} P\_{O\_2}^{0.5}} \tag{20}$$

$$E\_0 = 1.2723 - 2.7645 \times 10^{-4} T\_{PEN} \tag{21}$$



#### 3.1.3. Thermal Balance Model

For thermal energy calculations, the lumped capacitance method was used. Homogeneous temperature among the fuel cell components was assumed, so the temperature of the fuel cell was considered the same as *TPEN*. Heat absorption and release related to the reactions inside the fuel cell and convective heat transfer between the fuel cell and fuel/air bulk flows were considered in this model. Consequently, three temperatures were achieved from the thermal balance model: air, fuel bulk flow temperature and PEN temperature. The relevant equation for thermal balances in the PEN structure appears as follows.

$$\begin{split} \frac{\partial T\_{\rm PEN}}{\partial t} &= -\frac{q\_{\rm f\mathcal{E}}\boldsymbol{\mu}}{\rho\_{\rm PEN}\varepsilon\_{p,\rm PEN}L\_{\rm c}w\_{\rm c}t\_{\rm PEN}} - \frac{q\_{\rm f\mathcal{E}}}{\rho\_{\rm PEN}\varepsilon\_{p,\rm PEN}L\_{\rm c}w\_{\rm c}t\_{\rm PEN}}}{ + \frac{1}{\rho\_{\rm PEN}\varepsilon\_{p,\rm PEN}}\Delta H\_{\rm R1}R\_{\rm R1} - \Delta H\_{\rm R2}\varepsilon\_{\rm d}R\_{\rm R2} - \frac{1}{\rho\_{\rm PEN}\omega\_{\rm c}L\_{\rm c}}\Delta H\_{\rm R5}R\_{\rm R5}}} \\ \end{split} \tag{22}$$

For the calculation of heat transfer at gas channels, the Nusselt number (*Nu*) was considered to be a constant value of 3.39. Anode and cathode channels can act as fins of uniform rectangular cross-sectional areas, so relations for heat transfer at extended surfaces are used, as shown below, where *η<sup>f</sup>* is the fin efficiency and the convective heat transfer coefficient *hch* is obtained from *Nu*.

$$q\_{\mathcal{g}} = \eta\_f h\_{\rm ch} \left( N\_{\rm ch} A\_f + A\_b \right) \left( T\_{\rm PEN} - T\_{\rm g} \right) \tag{23}$$

#### *3.2. Steam Methane Reformer*

The methane steam reforming process requires a large amount of heat; thus, exhaust gas from the system is generally used as a heat source for the ESR, namely, an allothermal reformer. The reference and AOGR #1 systems suggested in this paper adopted this type of ESR. However, DIR stacks can mitigate the demand for high-performance ESR [4]. ESR can operate at a lower temperature range than the normal operating temperature

(973.15–1073.15 K) [32]. Hence, an adiabatic reformer was also examined in the AOGR #2 system. A shell and tube type catalytic steam reformer was selected in this paper, and the suggested specifications of the reformer are detailed in Table 5. The model consists of mass and thermal balance models. Gas mixtures are considered ideal gases, and the porosity of the bed is constant.

**Table 5.** Reformer specifications.


For catalytic steam reformers, it is important to find the most appropriate catalyst because it directly affects the performance of the reformer. Among various active metals, nickel (Ni) is widely used because of its high reactivity and long durability [33,34]. Xu and Froment examined the kinetics of the SMR process with a Ni/MgAl2*O*<sup>4</sup> catalyst under an operating range of approximately 675–1000 K. The results are the most widely used for SMR kinetcis [35]. The SMR model in this paper has been developed based on the Xu and Froment model. While various reactions take place in catalytic SMR, only the SMR reaction, the WGS reaction, and the direct steam reforming (DSR) reaction are considered. The SMR (R1) and WGS (R2) reactions are defined in Table 3, and the DSR (R6) reaction is presented below. It is assumed that the species in the reactant flow are CH4, H2O, CO, H2 and CO2. The SMR kinetics based on Xu and Froment model are organized in Table 6 [35].

Direct steam reforming reaction (DSR) *CH*<sup>4</sup> + 2*H*2*O* → *CO*<sup>2</sup> + 3*H*<sup>2</sup> (R6) (24)



Mass and thermal balances in both the gas and solid phases are presented below. Equations (25) and (26) represent the mass balance in the gas and solid phases for each species *i*, respectively [28]. The mass transfer coefficient (*kg*,*i*) is presented in Equations (27) and (28) [36].

$$\varepsilon\_{b}\frac{\partial \mathbb{C}\_{\mathbb{S},i}}{\partial t} = -v\_{\mathbb{S}}\frac{\partial \mathbb{C}\_{\mathbb{S},i}}{\partial \boldsymbol{x}} - k\_{\mathbb{S},i}a\_{\mathbb{v}}\left(\mathbb{C}\_{\mathbb{S},i} - \mathbb{C}\_{\mathbb{s},i}\right) \tag{25}$$

$$\frac{d\mathbb{C}\_{\rm s,i}}{dt} = -k\_{\rm g,i}a\_{\rm v}(\mathbb{C}\_{\rm g,i} - \mathbb{C}\_{\rm s,i}) + (1 - \varepsilon\_{\rm b})\rho\_{\rm cut}r\_{\rm SMR,i} \tag{26}$$

$$k\_{\mathbb{S},i} = j\_{D,i} \text{Re} \text{Sc}\_i^{1/3} \frac{D\_i}{d\_p} \tag{27}$$

$$
\varepsilon\_b j\_{D,i} = 0.765 \text{Re}^{-0.82} + 0.365 \cdot \text{Sc}\_i^{-0.398} \tag{28}
$$

Equations (29) and (30) are the thermal balances in the gas and solid phases. The convective heat transfer coefficient (*hg*,*SMR*) in Equation (30) is known from the Chilton– Colburn j-factor (*jH*). For forced convection through a packed bed, Yoshida et al. [37] suggested empirical correlations of *jH*. Ψ is an empirical coefficient depending on the particle shape, and its value is 1 for a sphere.

$$\varepsilon\_b \frac{\partial T\_\%}{\partial t} = -\upsilon\_\% \frac{\partial T\_\%}{\partial \mathcal{X}} + \frac{h\_\text{\mathcal{g},SMR} a\_\text{\mathcal{v}}}{\rho\_\text{\mathcal{g}} c\_{p,\text{\mathcal{g}}}} \left(T\_\text{\mathcal{v}} - T\_\text{\mathcal{S}}\right) \tag{29}$$

$$\frac{\partial T\_{\sf s}}{\partial t} = -\frac{h\_{\sf S,SMR}a\_{\sf v}}{\rho\_{\sf bcd}c\_{p,\sf b}} \left(T\_{\sf s} - T\_{\sf S}\right) + \frac{(1-\varepsilon\_{\sf b})\rho\_{\sf cat}}{\rho\_{\sf b}c\_{p,\sf b}} \sum\_{k} (-\Delta H\_{k}\eta\_{k}R\_{k}) \left(k = \text{R1, R2, R6}\right) \tag{30}$$

$$h\_{\mathbb{S},SMR} = j\_H \frac{c\_{p,\mathbb{S}} \rho\_{\mathbb{S}} v\_{\mathbb{S}}}{Pr^{2/3}} \tag{31}$$

$$j\_H = 0.91 \text{Re}^{-0.51} \Psi \ 0.01 < Re < 50 \tag{32}$$

$$j\_H = 0.61 \text{Re}^{-0.41} \Psi \text{ 50} < \text{Re} < 7000 \tag{33}$$

#### *3.3. Catalytic Combustor*

To compute the temperature and species in the CC, a mathematical model including mass and thermal balances was developed. A Pt-catalyzed monolithic combustor is analyzed. Specifications of CC are given in Table 7. The oxidation reactions over Pt considered in this model are CO, CH4 and H2 oxidation. The rate expressions and reaction rate per Pt surface area are organized in Table 8 based on Ref. [38]. Chemical reactions are assumed to occur only on the external surface of the catalytic wall.

**Table 7.** Catalytic combustor specifications.



**Table 8.** Equations for the reaction rate of CC based on Ref. [38].

Homogeneous temperature, concentration and velocity within the channel are assumed for mass and thermal balance computations. Equations (34) and (35) describe the mass balances for the gas and solid phases for each species *i*, where *i* refers to CH4, H2O, H2, CO, CO2, O2 and N2. The corresponding equations for thermal balances for the gas and solid phases are presented in Equations (36) and (37). For the heat transfer coefficient (*hg*,*CC*) calculation, Nu is considered to have a constant value of 3.39.

$$
\varepsilon\_m \frac{\partial y\_{\mathcal{G},i}}{\partial t} = -v\_{\mathcal{G}} \frac{\partial y\_{\mathcal{G},i}}{\partial x} - k\_{m,i} S \left( y\_{\mathcal{G},i} - y\_{s,i} \right) \tag{34}
$$

$$(1 - \varepsilon\_m)\frac{dy\_{s,i}}{dt} = k\_{m,i}S(y\_{\mathcal{S},i} - y\_{s,i}) + \frac{RT\_{\text{CC}}}{P\_{\text{tot}}}aR\_{\text{CC},i} \tag{35}$$

$$
\varepsilon\_{\rm m} \rho\_{\rm \mathcal{S}} c\_{p, \rm \mathcal{S}} \frac{\partial T\_{\rm \mathcal{S}}}{\partial t} = h\_{\rm \mathcal{S}, \rm CC} S \left( T\_{\rm \mathcal{S}} - T\_{\rm \mathcal{S}} \right) - v\_{\rm \mathcal{S}} \rho\_{\rm \mathcal{S}} c\_{p, \rm \mathcal{S}} \frac{\partial T\_{\rm \mathcal{S}}}{\partial x} \tag{36}
$$

$$(1 - \varepsilon\_{\rm ff})\rho\_{\rm s}c\_{p,\rm s}\frac{\partial T\_{\rm s}}{\partial t} = h\_{\rm g,\rm C\mathcal{C}}S\left(T\_{\mathcal{S}} - T\_{\rm s}\right) + a\left(\sum\_{i} (-\Delta H\_{i})R\_{i}\right) \tag{37}$$

#### *3.4. Balance of Plant*

3.4.1. Air and Recirculation Blower

A model of air and recirculation blowers was employed to examine the outlet temperature and power consumption. The stack air flow rate was determined by the air utilization factor, and the CC air flow rate was controlled to maintain a CC temperature below 1123.15 K. In terms of the recirculation blower, the target recirculation ratio determines the flow rate. The blower outlet temperature and power consumption are obtained from Equations (38) and (39). In the computation, the values of *ηisen*, *ηmotor* and *ηmech* were 0.8, 0.9 and 0.9, respectively.

$$T\_{blower,o} = T\_{blower,i} - T\_{blower,i} \frac{1 - \left(p\_o/p\_i\right)^{\kappa - 1/\kappa}}{\eta\_{isen}} \tag{38}$$

$$P\_{blower} = \dot{m}\_{blower}c\_{p,air} \left(T\_{blower,o} - T\_{blower,i}\right) \eta\_{mator} \eta\_{mech} \tag{39}$$

### 3.4.2. Heat Exchanger

Heat exchangers are employed to recover heat from AOG, COG and exhaust gas. They act as fuel/air preheaters, a steam generator and a HR-HE. The outlet temperature of each gas and heat transfer rate are defined by the effectiveness-NTU method, as shown in

Equation (40) since the information is insufficient to use the LMTD method. *εHE* indicates the effectiveness of a heat exchanger and was set to 0.75.

$$q\_{HE} = \varepsilon\_{HE} \mathbb{C}\_{\min} (T\_{h,i} - T\_{c,i}) \quad \mathbb{C}\_{\min} = \min \{ \mathbb{C}\_{hot}, \mathbb{C}\_{cool} \} \tag{40}$$

### *3.5. Performance Factor*

*Uf uel* and *Uair* indicate the fuel and air utilization factors, respectively. Methane is only supplied as fuel in this model. If there are other hydrocarbons in fuel, they should be added to the denominator in Equation (41). The fuel and air flow rates of the stack are determined from the target utilization factor. Recirculated hydrogen, however, is not considered when the fuel flow rate is calculated.

$$
\Delta I\_{fuel} = \frac{I}{8F\dot{n}\_{CH\_4}}\tag{41}
$$

$$\mathcal{U}\_{air} = \frac{J}{4Fy\_{O2}\dot{n}\_{air}}\tag{42}$$

The electrical, thermal and total efficiencies are shown in Equations (43)–(45). The electrical efficiency of the system is the ratio of the net generated energy of the system to the chemical energy of the supplied fuel. To estimate the thermal efficiency, the temperature of the final exhaust gas from the CC (*Tvent*) is assumed to become 393.15 K after the HR-HE.

$$
\eta\_{elc} = \frac{P\_{stack} - P\_{FCRlower} - P\_{CCBllower}}{\dot{m}\_{CH\_4} LHV\_{CH\_4}} \tag{43}
$$

$$\eta\_{\rm th} = \frac{\mathcal{C}\_{p, \rm COOG} \dot{m}\_{\rm COOG} (T\_{\rm COOG,O} - T\_{\rm rent})}{\dot{m}\_{\rm CH\_4} LHV\_{\rm CH\_4}} \tag{44}$$

$$\eta\_{\rm tot} = \frac{\eta\_{\rm elle} + \eta\_{\rm th}}{\dot{m}\_{\rm CH\_4} LHV\_{\rm CH\_4}} \tag{45}$$

### **4. Results and Discussion**
