*2.1. Computational Model*

The fuel cell stack consists of multiple fuel cell units stacked together. The cooling plates are distributed at both ends of a single fuel cell and are in close contact with the bipolar plate. The heat generated during PEMFC operation enters the cooling plate through heat conduction in the bipolar plate, and then the heat is taken away by the coolant circulation in the cooling plate. Figure 1 shows the structure of the fuel cell stack.

**Figure 1.** Single cell structure of PEMFC.

Figure 2 shows the cooling plate model to be calculated. Heat is transported from both sides of the cooling plate during its actual working process. The cooling plate is divided from the central plane according to the cooling plate's symmetry for ease of calculation, and the half model of the overall cooling plate is analyzed to simplify the calculation. The heat flux acts on the bottom, and the value is a fixed value of 5000 W/m2, which is a common value encountered during typical PEMFC operation. The heat produced by PEMFCs is comparable to the output cell power (with PEMFCs with a rated power of 1 kW, around 1–1.5 kW of heat is produced) [36].

**Figure 2.** Calculation model of cooling plate.

Five different cooling channels are designed, as shown in Figure 3. Among them, model 1 is a multi-serpentine flow field, model 2 is a multi-turn flow field, model 3 is a multi-helical flow field, model 4 is a flow field with a uniform plate, and model 5 is a honeycomb structure flow field. The parameters of the geometric structure are shown in Table 1.

**Figure 3.** Design scheme of flow fields.

**Table 1.** Model parameters of flow fields.


The steady-state calculation formula for calculating the heat flow of a double cooling plate is as follows:

$$q = \frac{Q}{2A} \tag{1}$$

The cooling plate has two sides for heat transfer. For the *n*-cell stack with current *I*, when all the reaction enthalpies of the fuel cell are converted into electric energy and the aquatic product is water vapor,

$$Q = nI\left(-\Delta h\_{\rm f}^{0}/2F - V\right) \tag{2}$$

where *n* is the number of cells, *I* is the cell current, *V* is the output voltage of the cell and *A* is the total area of the cell, Δ*h*<sup>0</sup> <sup>f</sup> is the enthalpy of water formation, and *F* is the Faraday constant.

The regional uniformity index of the area-weighted variable *γ*<sup>a</sup> is calculated using the following formula:

$$\gamma\_{\mathbf{a}} = 1 - \frac{\sum\_{\mathbf{i}=1}^{\mathrm{n}} \left[ (\left| \phi\_{\mathbf{i}} - \overline{\phi}\_{\mathbf{a}} \right|) A\_{\mathbf{i}} \right]}{2 \left| \overline{\phi}\_{\mathbf{a}} \right| \sum\_{\mathbf{i}=1}^{\mathrm{n}} A\_{\mathbf{i}}} \tag{3}$$

*φ*<sup>a</sup> is the average of the variables across the surface:

$$\overline{\phi\_{\mathbf{a}}} = \frac{\sum\_{i=1}^{n} \phi\_{\mathbf{i}} A\_{\mathbf{i}}}{\sum\_{i=1}^{n} A\_{\mathbf{i}}} \tag{4}$$

where *γ*<sup>a</sup> is the uniformity index, *φ* is variable across the surface, *A* is a superficial area, *i* is the mesh face index with *n* mesh faces, and *n* is the number of grids.
