Three Dimensional Structure Optimization of Proton Exchange Membrane Fuel Cell with Radial Flow Field Based on Topology Optimization

Abstract

This article uses the topological optimization method to optimize the flow field of radial proton exchange membrane fuel cells. Due to the high computational cost of this method, the model was simplified to a two-dimensional model, and the cathode flow field was studied under isothermal and steady-state conditions. Using the gradient-based algorithm SNOPT, the flow field evolves freely in the sector design domain to maximize battery power and minimize reaction energy loss. Perform 3D modeling of the obtained optimized model and optimize the longitudinal depth by applying different tilt angles and obstacles. The results indicate that the application of obstacles has a significant optimization effect on speed. Setting the gradient has a better effect on the oxygen concentration distribution. With the increase of the gradient, the uniformity of the oxygen concentration distribution is also increased. As the inclination increases, the pressure difference of the reaction gas gradually increases, and the pressure uniformly decreases. After setting obstacles, the pressure of the reaction gas decreases step by step after passing through each obstacle. The effect of increasing the inclination to increase the average current density is better than setting obstacles, and the comprehensive effect of setting an inclination of 2° is the best at working voltages of 0.5 V and 0.6 V.

1. Introduction

The fuel cell is an important conversion medium in hydrogen energy, and in new energy, hydrogen energy has attracted people’s attention as a clean secondary energy, known as the “ultimate energy” of the 21st century. Proton exchange membrane fuel cell (PEMFC) have the characteristics of high power density, low operating noise, and pollution-free products, making them highly competitive in future transportation power systems.^1,2 The electrode plate is an important component of PEMFC, which can connect adjacent monomers and support the stack structure. Its surface design will affect the distribution of reaction gases and hydrothermal management, affecting the performance of the battery.³ Therefore, the study of bipolar plates is particularly important for improving the performance of PEMFC. In this paper, the topology optimization method is used to find the optimal transverse channel layout in the radial flow field design domain of PEMFC, and the gradient and obstacles are applied to the longitudinal depth for optimal design, and compared with the traditional flow field.

Topology optimization (TOM) is a structural optimization method based on load conditions, constraints (such as stress, displacement, frequency, weight, etc.), and performance indicators (stiffness, weight, etc.), using finite element analysis to make the design domain reach the optimal material layout for simulation calculation.⁴ TOM is suitable for solid, fluid, and fluid-solid coupling devices, and can meet multi-objective and multi constraint conditions. It can be flexibly applied to various complex working conditions. Because the model after topology optimization is difficult to produce directly, it is usually used as a preliminary design in design, and further artificial optimization design is carried out according to production and actual conditions.

Andreas⁵ et al. use the gradient penalty of elements to obtain stability and prevent chessboard, while the very fine structure can solve the sharp problems in topology optimization through the adaptive mesh division of the material void interface, and obtain the same qualitative structure and quantitative stiffness as the uniform mesh. Zhang⁶ et al. proposed an effective improved evolutionary topology optimization (M-ETO) method to design hyperelastic structures with large deformation. This improvement enables nonlinear topology optimization problems to be solved at a relatively large evolution rate. Wang⁷ et al. designed a novel porous structure optimization method based on local volume constraints for TOM, which specifies the topological skeleton of the stress field as a density field as a solid simulation element. These elements guide the deposition of materials around them and can also be reshaped or removed during the optimization process, improving the convergence and optimization speed of information optimization. Pedro⁸ et al. proposed a two-step solution method that decouples the overall aeroelastic performance target from the solid-air topological structure on the search structure, solving the problem of reference fully coupled fluid-structure, and providing numerical examples of compliant wings with performance targets at two Mach numbers. Pollini⁹ designed a density-based topology optimization method to solve the problem of wind farm layout optimization and maximized its energy production under the condition of limiting the minimum and maximum number of wind turbines placed and the minimum spacing between wind turbines.

Topology optimization is also widely used in redox batteries. Mathieu-Potvin¹⁰ et al. used the gradient-based algorithm for topology optimization of the platinum density in the cathode catalyst layer to improve the utilization rate of the catalyst, to maximize the current density of the total fixed platinum. It was found that under the high stoichiometric ratio, platinum is concentrated in the near film layer with a high reaction rate. Xun¹¹ et al. developed a single-stage power conversion topology that can directly couple energy storage devices, applying a fluctuating voltage topology, thereby reducing costs and making it more compact without a DC/DC converter. The feasibility and effectiveness of this topology structure were demonstrated through simulation results of fuel cell supercapacitor powertrain modeling. Kim¹² et al. described the derivation and solution of a new mathematical equation for topology optimization, which combines density-based algorithms, moving asymptote interpolation (MMA), and incompressible Navier Stokes equations to design the topology of parallel flow fields, and studied the influence of key design parameters and pressure drop on topology optimization. Zhang¹³ et al. analyzed the stress and displacement distribution of the end plate through numerical simulation. Three optimization topologies were developed to minimize flexibility, uniform stress distribution, and two objective couplings. The inlet end plate and blind end plate were reconstructed respectively, reducing the mass of the optimized inlet end plate by 35 % and the weight of the optimized blind end plate by 46 %. Reza¹⁴ et al. used the average depth model, designed 3DPEMFC into a 2D model, and carried out topology optimization for 2D models of the parallel flow field. Compared with the traditional flow field, the topology optimization design significantly enhanced the power generation performance, improved the reactant distribution, and reduced the pressure drop.

The radial flow field can quickly and uniformly distribute the reaction gas, and timely discharge the water generated by the reaction, thereby improving the output power of PEMFC and ensuring its reliability under various working conditions, thus having great development space. Nakanishi,¹⁵ Diethelm,¹⁶ and others first proposed a kind of radial flow field structure in the case of building a solid oxide fuel cell (SOFC) stack. They pointed out in their research that this structure can effectively solve the problem of brittle bipolar plates in the process of fuel cell stack assembly. Batawi¹⁷ adopted a radial structure design in the high-temperature fuel cell stack and for the first time integrated the exhaust gas treatment device of the post-combustion process onto the bipolar plate, greatly reducing the volume and manufacturing cost of the stack. To improve the water and gas transmission capacity on the cathode side, Schuler¹⁸ used various circular bipolar plates with labyrinth flow channel structures in his case and confirmed the superior performance of such structures in gas transportation, distribution, exhaust emissions, and other aspects through subsequent research. Doggwiller^19,20 used a lattice structure flow field in his case and compared it with a maze structure flow field. The results indicate that the flow field using a lattice structure also exhibits excellent material transport performance. Hernandez Guerrero²¹ selected a quarter of the radial flow field area and simulated the current density distribution characteristics and the mass transfer speed in GDL of PEMFC with 4-channel, 8-channel, and 12-channel structures respectively. The results show that the three structures show significant performance improvement compared with the traditional parallel flow field, serpentine flow field, and other structures. Vladimir²² integrates the characteristics of the radial flow field and interdigital flow field, forcing gas to pass through the diffusion layer under pressure to participate in the reaction through discontinuous flow channel settings, thereby improving gas utilization and battery power density, ultimately achieving the goal of improving battery performance. The radial flow field structure proposed by Perez²³ can avoid large pressure drops, but this type of structure is too complex. When one of the outlets is blocked by water, it often leads to problems such as water flooding and uneven pressure distribution. Friess²⁴ et al. designed a novel cathode radial flow field that uses a flow control ring to distribute gas more evenly. This design can effectively improve battery performance and reduce the pressure drop by nearly 50 % compared to the serpentine flow field. Zhang²⁵ et al. designed a new rim-shaped radial flow field structure and studied the effects of layer number, number of branches, and waveform amplitude of branches on PEMFC. It was found that the battery performance was best when there were 8 branches, 4 layers, and 15° branch curvature.

This study is based on the average depth model, simplifying the 3D model to a 2D model, and using topology logic optimization to design the model. Design a TOM model to minimize the power loss and energy dissipation of the battery. The resulting TOM model is modeled in three dimensions and optimized in terms of longitudinal depth by applying different inclinations and obstacles.

2. Numerical Simulation

In PEMFC composed of the cathode side and anode side, oxygen is reduced and hydrogen is oxidized. Topology optimization can reduce computer research time and computing costs. In this study, we will take the cathode reaction zone of radial PEMFC composed of cathode plate and cathode gas diffusion layer (GDL) as the research object, and other influencing factors of the fuel cell layer as the boundary conditions and feature settings of the model. COMSOL Multiphysics^{® 26} is commonly used for finite element analysis in the field of fuel cells, and can also be used for simulation design of topology optimization. However, the topology optimization problem of 3D PEMFC is more complex, so the average depth model can be used to simplify the complex problem of 3D PEMFC into a two-dimensional plane structure model. Using a fully coupled solver to solve the finite element method, the convergence of the nonlinear equation is improved through multiple iterations and gradual approximation, and the problem of coupling multiple physical fields is solved. The optimization solver optimizes the flow field structure of radial PEMFC and terminates the optimization iteration calculation when the maximum number of iterations is reached.

2.1 Model assumptions

This article adopts the following assumptions in constructing the PEMFC model

1. 1.    Steady state, laminar flow;
2. 2.    The reaction gas is compressible at the inlet;
3. 3.    Maintaining a constant model temperature is controlled by an external controller
4. 4.    The single-phase weighted hybrid model is assumed to reduce the calculation cost of the topology optimization model. Considering the small height of the channel and GDL, water management can be simplified no matter how the flow field is designed,²⁷
5. 5.    Porous materials are isotropic with uniform porosity;
6. 6.    Assuming that the membrane is impermeable;
7. 7.    Neglecting water transport through membranes;
8. 8.    Set reactant gas and water vapor as an ideal gas.

2.2 Simplified model

To simplify the 3D model, the same transmission conditions need to be used when using the average depth model. The three-dimensional PEMFC model and the two-dimensional model using the average depth model are shown in Fig. 1. Based on the material characteristics of the channel or rib plate in the three-dimensional model and GDL, the equivalent porosity ($\epsilon $), permeability (κ), and diffusion coefficient (D) of the porous channel (ch) or porous rib (rib) in the two-dimensional model can be calculated as follows:   

\begin{equation} \bar{P}_{\text{a}} = \frac{P_{\text{gdl}}H_{\text{gdl}} + P_{\text{a}}H_{\text{a}}}{H_{\text{gdl}} + H_{\text{ch}}},\ P = \{\epsilon,\kappa,D\},\ \text{a} = \{\text{ch},\text{rib}\} \end{equation} (1)

Figure 1.

3D radial PEMFC and 2D model using average-depth model.

Where, H is the height of each layer of PEMFC, $\bar{P}$ is the average property. Equivalent porosity ($\epsilon (\gamma )$), Brinkman constant (α(γ)), and diffusion coefficient (D(γ)) will be taken as the design parameters for topology optimization.

2.2.1 Conservation of mass and momentum

The Brinkman equation, which combines the continuity equation and momentum equation, describes the flow in porous media as follows:   

\begin{equation} \nabla (\rho u) = 0 \end{equation} (2)

  

\begin{align} &\frac{1}{\epsilon(\gamma)^{2}}\rho(u\nabla )u \notag\\ &\quad= -\nabla p + \frac{1}{\epsilon(\gamma)^{2}}u\nabla^{2}u - \frac{12\mu}{(H_{\text{gdl}} + H_{\text{ch}})^{2}}u - \frac{\mu}{\kappa(\gamma)}u \end{align} (3)

Where, u is the velocity vector, p is the relative pressure, ρ is the gas density and μ is the dynamic viscosity. Porosity $\epsilon (\gamma )$ is the first design-dependent feature. The third term of the right momentum equation considers the velocity gradient perpendicular to the two-dimensional plane.²⁸ The last item about the right side of the equation is Darcy velocity, which represents the viscous resistance of fluid in porous media, where κ is permeability. Therefore, the following Brinkman constant α(γ) will be the second design-related feature:   

\begin{equation} \alpha(\gamma) = \frac{\mu}{\kappa(\gamma)} \end{equation} (4)

2.2.2 Conservation of species

The equation of species conservation is:   

\begin{equation} (H_{\text{gdl}} + H_{\text{ch}})\rho(u\nabla )\omega_{\text{k}} = (H_{\text{gdl}} + H_{\text{ch}})\nabla j_{\text{k}}(\gamma) + R_{\text{k}} \end{equation} (5)

The convection term on the left consists of the mass fraction ω_k as a species, k = {O₂, N₂, H₂O}. The diffusion term on the right side of the equation term explains the species diffusion mass flux j_k(γ), which is approximately defined by mixed average diffusion:   

\begin{equation} j_{\text{k}}(\gamma) = D(\gamma)\left(\rho D_{\text{k}}^{\text{m}}\nabla \omega_{\text{k}} + \rho\omega_{\text{k}}D_{\text{k}}^{\text{m}}\frac{\nabla M}{M}\right) \end{equation} (6)

Where, D(γ) is the third design-related characteristic factor representing the diffusion coefficient. Under the assumption of isobaric and isothermal conditions, $D_{\text{k}}^{\text{m}}$ in the Eq. 6 can be used to derive by Maxwell Stefan equation:²⁹   

\begin{equation} D_{\text{k}}^{\text{m}} = \cfrac{1 - \omega_{\text{k}}}{\displaystyle\sum\nolimits_{\text{j${\neq}$k}}^{\text{n}_{\text{s}}}\cfrac{\omega_{\text{j}}}{D_{\text{kj}}}} \end{equation} (7)

Where n_s is the number of species, D_kj can be calculated using the reference value of the binary diffusion coefficient:   

\begin{equation} D_{\text{kj}} = D_{\text{kj}}^{\text{ref}}\left(\frac{P_{\text{ref}}}{P}\right)\left(\frac{T}{T_{\text{ref}}}\right)^{1.75} \end{equation} (8)

$D_{\text{kj}}^{\text{ref}}$, P_ref and T_ref are provided in Table 1.

Table 1. PEMFC performance and topological optimization parameters.

Symbol	Value	Description
Dimensional parameters
Rin	5 × 10−4	Inlet radius/m
Rout	1.5 × 10−2	Outlet radius/m
Hch	1 × 10−3	Channel height/m
Hgdl	5 × 10−4	GDL height/m
βgeo	60	Section angle/deg
PEMFC properties and parameters
Vcell	0.6	Cell voltage/V
Voc	1.23	Open-circuit voltage/V
T	343.15	Cell temperature/K
Tref	298.15	Reference temperature/K
Pref	101325	Reference pressure/atm
Rlump	7 × 10−5	Lumped resistance/Ω m2
Uin	2	Inlet velocity/m s−1
$\omega_{\text{O}_{2}}^{\text{in}}$	0.228	Inlet mass fraction of O2
$\omega_{\text{H}_{2}\text{O}}^{\text{in}}$	0.023	Inlet mass fraction of H2O
$\omega_{\text{N}_{2}}^{\text{in}}$	0.749	Inlet mass fraction of N2
$C_{\text{O}_{2}}^{\text{ref}}$	30	Reference concentration of O2/mol m−3
αc	0.5	Cathode transfer coefficient
i0	0.17	Exchange current density/A m−2
ns	3	Number of species k = {O2, N2, H2O}
ne	4	Number of electrons
$M_{\text{H}_{2}\text{O}}$	0.018	Molar mass of H2O/kg mol−1
$M_{\text{N}_{2}}$	0.028	Molar mass of N2/kg mol−1
$M_{\text{O}_{2}}$	0.032	Molar mass of O2/kg mol−1
$\nu_{\text{H}_{2}\text{O}}$	2	Chemometric coefficient of H2O
$\nu_{\text{N}_{2}}$	0	Chemometric coefficient of N2
$\nu_{\text{O}_{2}}$	−1	Chemometric coefficient of O2
$D_{\text{O}_{2},\text{H}_{2}\text{O}}^{\text{ref}}$	2.82 × 10−5	O2 diffusion coefficient in H2O/m2 s−1
$D_{\text{N}_{2},\text{H}_{2}\text{O}}^{\text{ref}}$	2.56 × 10−5	N2 diffusion coefficient in H2O/m2 s−1
$D_{\text{O}_{2},\text{N}_{2}}^{\text{ref}}$	2.2 × 10−5	O2 diffusion coefficient in N2/m2 s−1
κgdl	1.18 × 10−11	Permeability of GDL/m2
$\epsilon_{\text{gdl}}$	0.4	Porosity of GDL
μ	2.07 × 10−5	Dynamic viscosity/Pa s
ρ	1	Gas density/kg/m3
Rg	8.3144	Universal gas constant/J mol−1 K−1
F	96485.3	Faraday constant (constant)/C mol−1
Average-depth model parameters
κrib	3.8 × 10−12	Equivalent permeability of porous ribs/m2
κch	1000	Equivalent permeability of porous channels/m2
$\epsilon_{\text{rib}}$	0.22029	Equivalent porosity of porous ribs
$\epsilon_{\text{ch}}$	0.94493	Equivalent porosity of porous channels
Drib	0.19703	Equivalent diffusivity of porous ribs
Dch	0.92167	Equivalent diffusion coefficient of porous channels
Topology optimization parameters
q	1	Convex factor
p	4	Penalty factor
τmin	3.5 × 10−4	Filter radius/m
β	2, 4, 6, 8, 10	Projection slope
θβ	0.5	Projection point
γ0	0.5	Initial value of control variable
topt	1 × 10−7	Optimum tolerance
evmax	5000	Maximum number of model evaluations
tin	0.001	Internal tolerance factor
Ci	1000	Constrained penalty factor
Nin	10	Maximum internal iteration
Vmax	0.7	Maximum volume fraction

The average molar mass M can be written as:   

\begin{equation} \frac{1}{M} = \sum\nolimits_{\text{k}=1}^{\text{n}_{\text{s}}}\frac{\omega_{\text{k}}}{M_{\text{k}}} \end{equation} (9)

R_k in Eq. 5 represents the consumption or production of species given by¹⁴   

\begin{equation} R_{\text{k}} = \frac{M_{\text{k}}\nu_{\text{k}}i_{\text{loc}}}{\text{n}_{\text{e}}F} \end{equation} (10)

Where, M_k is the molar mass of the species, n_e is the number of electrons, ν_k is the stoichiometric coefficient of the reaction, and F is the Faraday constant.

2.2.3 Reaction kinetics

The reduction reaction of oxygen in the porous cathode is as follows:   

\begin{equation} \text{O}_{2} + \text{4H$^{+}$} + \text{4e$^{-}$} \leftrightarrow \text{2H$_{2}$O} \end{equation} (11)

Set the reaction that occurs in the electrolyte and electrode catalytic layer as the cathode boundary condition. This boundary condition means the mass flux due to the electrochemical reaction R_k. The electrochemical reaction is measured by i_loc means of local current density and Butler-Volmer equation according to the following dynamic expression:¹⁴   

\begin{equation} i_{\text{loc}_{\text{BV}}} = i_{0}\frac{\text{C}_{\text{O${_{2}}$}}}{\text{C}_{\text{O${_{\text{2ref}}}$}}}\left(\exp\left(\frac{\alpha_{c}}{R_{\text{g}}T}F\eta_{c}\right) - \exp\left(-\frac{1 - \alpha_{c}}{R_{\text{g}}T}F\eta_{c}\right)\right) \end{equation} (12)

The cathode overpotential η_c is calculated by:   

\begin{equation} \eta_{c} = V_{\text{OC}} - V_{\text{cell}} - R_{\text{lump}}i_{\text{loc}} \end{equation} (13)

Where lumped resistance R_lump is defined to describe voltage losses at the anode, MEA, and GDL due to charge transfer. It is noteworthy that $i_{\text{loc}} = i_{\text{loc}_{\text{BV}}}$. Eqs 12 and 13 are highly non-linear. To solve the current density, an additional differential equation can be defined:   

\begin{equation} \frac{\partial i_{\text{loc}}}{\partial i} = i_{\text{loc}} - i_{\text{loc}_{\text{BV}}} \end{equation} (14)

Among them, i_loc is the variable that needs to be solved and realized by using the ODE interface in COMSOL.

2.3 Topological optimization

Topology optimization is used to design the optimal flow field of the circular radial flow field, and a gradient-based optimization method is used. The sparse nonlinear optimization SNOPT is suitable for topology optimization problems with nonlinear dependent design variables. When the relative change value of the design variable is less than the optimality tolerance value set by the SNOPT solver, the SNOPT optimization solver ends with a value of 1 × 10⁻⁷. The parameters related to this algorithm are listed in Table 1. Because the discrete design variables have no well-defined derivative information, the optimization is based on continuous design variables, which makes the topology optimization problem ill-posed. Therefore, the material model allocates intermediate values between porous ribs (0) and porous channels (1).

2.3.1 Material modeling and regularization

The material model can describe continuous design variables from 0 to 1, and its formula defines the degree of sparsity. The material model of this study applies a specific method to each design parameter based on density. Indicator γ represents a continuous design variable field used to interpolate three design-related attributes $\epsilon (\gamma )$, α(γ), and D(γ). γ $ \in $ [0, 1] indicates porous rib (0), porous channel (1), or any intermediate value between 0 and 1 (0, 1). The latter means that the average properties can continuously change from porous ribs to porous channels, in an intermediate state of fluid-solid coupling.

First, porosity $\epsilon (\gamma )$ of the linear interpolation is as follows:¹⁴   

\begin{equation} \epsilon(\gamma) = \epsilon_{\text{rib}} + (\epsilon_{\text{ch}} - \epsilon_{\text{rib}})\gamma \end{equation} (15)

Among them, $\epsilon_{\text{ch}}$ and $\epsilon_{\text{rib}}$ are the average porosity of porous channels and ribs calculated by Eq. 1.

Second, the constant Brinkman α(γ) is discussed in the equation. Using the convex interpolation function³⁰ to interpolate:   

\begin{equation} \alpha(\gamma) = \alpha_{\text{rib}} + (\alpha_{\text{ch}} - \alpha_{\text{rib}})\frac{1 - \gamma}{1 + \text{q}\gamma} \end{equation} (16)

Where, the convexity factor q is an adjustable parameter related to the reverse permeability of the intermediate design variable. α_ch and α_rib are scalar values derived from Eqs. 1 and 4 represent the Brinkman constants for porous channels and porous ribs, respectively.

Third, the diffusion coefficient factor D(γ), expressed in Eq. 6, is interpolated by the solid isotropic material penalty (SIMP) method given below:   

\begin{equation} D(\gamma) = D_{\text{rib}} + (D_{\text{ch}} - D_{\text{rib}})\gamma^{\text{p}} \end{equation} (17)

Where, p is the penalty coefficient of the SIMP method. D_rib is the diffusion coefficient factor of porous ribs and D_ch is the diffusion coefficient factor of porous channels. D_rib and D_ch are calculated by the depth averaging program Eq. 1.

Regularization can solve the problem of grid dependency. In this study, Hermann von Helmholtz filter³¹ is used as follows:   

\begin{equation} \gamma_{\text{f}} = \text{r}_{\text{min}}^{2}\nabla^{2}\gamma_{\text{f}} + \gamma \end{equation} (18)

Where, γ_f is the filtered design variable. r_min is the minimum radius of the filter and its minimum value must not be smaller than the size of the grid element. Because the radius of the filter is measured in meters, it overcomes the constraint of grid size. Although filtering simplifies the solution of optimization problems, it leaves a large gray area in the topology. These areas contain unphysical intermediate values and therefore wish to use projection techniques to reduce grayscale areas. This study used hyperbolic tangent projection:³²   

\begin{equation} \gamma_{\text{h}} = \frac{\tanh(\beta\theta_{\beta}) + \tanh(\beta(\gamma_{\text{f}} - \theta_{\beta}))}{\tanh(\beta\theta_{\beta}) + \tanh(\beta(1 - \theta_{\beta}))} \end{equation} (19)

Where, γ_h is the projection design variable, θ_β is the projection point, and β is the projection slope. Increasing the projection slope can reduce the area of intermediate states and obtain clearer topological structures, but the discretization of the optimization results will also be stronger. In the design, it is important to perform parameterized scans to approach discrete solutions. The parameterized scan values for β are defined as follows:   

\begin{equation} \beta = \{2,4,6,8,10\} \end{equation} (20)

Optimization starts with setting the initial design variable γ₀ = 0.5 for the β = 2. When the solver runs the following β will get the initial design variables from the last solution of the previous one. Each step of β runs 100 iterations, resulting in a total of 500 iterations to complete the topology optimization process.

2.3.2 Objective function and constraints

The formula for multi-objective topology optimization is as follows:   

\begin{equation} \min\limits_{\gamma}{:}\ \int\nolimits_{\varOmega}V_{\text{cell}} \cdot i_{\text{loc}}d\varOmega + \int\nolimits_{\varGamma}\left(p + \frac{1}{2}\rho v^{2}\right)vnd\varGamma \end{equation} (21)

  

\begin{equation*} \text{aim}{:}\ \frac{\displaystyle\int\nolimits_{\Omega}\gamma_{\text{i}}\text{d}\Omega}{\displaystyle\int\nolimits_{\Omega}1\text{d}\Omega} \leq \text{V}_{\text{max}} \end{equation*}

  

\begin{equation*} 0 \leq \gamma_{\text{i}} \leq 1 \end{equation*}

  

\begin{equation} \text{i} \in \{1,2,3,\ldots,\text{N}\} \end{equation} (22)

Multi-objective function Eq. 21 is to maximize the operating power of the battery while minimizing the energy loss of the reaction gas. The first function is the generated power, which is minimized by a negative sign, where V_cell is the battery voltage, and Ω represents the surface of the design domain. The second function is power dissipation, which calculates the net energy flux at the inlet and outlet of the battery. p and v are the pressure and velocity of the gas at the corresponding boundary. Γ is the domain boundary. ρ is defined by the ideal gas law $\rho = \frac{p_{\text{ref}}M}{R_{\text{g}}T}$, where R_g is the general gas constant, T is the battery temperature, and the average molar mass M is calculated by Eq. 9.

The volume constraint expression (Eq. 22) indicates that the ratio of the volume of the porous channel to the volume of the total domain cannot exceed a certain percentage. Therefore, V_max is the maximum volume fraction that limits the number of materials used in the optimization process. γ_i is an independent design variable between 0 and 1. N is the number of design variables.

2.3.3 Topology optimization process

Firstly, the preprocessing of the structural model is completed through COMSOL, including geometric modeling, design parameters, boundary conditions, and grid settings. The entire process is controlled through MATLAB, and the initial structure is analyzed using the finite element method. Then, the results obtained from the first step of analysis are extracted and used as initial design parameters for the second step of optimization. The optimization results are exported by the SNOPT solver.

3. Results and Discussion

In this work, the radial flow field of the cyclic PEM fuel cell was optimized to maximize power output and minimize dissipated energy. The rationality of this program was verified by comparing the finite element simulation results of the 2D model of the dual channel serpentine flow channel using the average depth model with its 3D single cathode PEMFC model. The obtained topology optimization model is used for 3D modeling, and different inclinations and obstacles are applied to optimize the design at the longitudinal depth level.

3.1 Model validation

We designed a dual channel serpentine flow field with the same design parameters by studying the average depth model published by Behroo¹⁴ et al. The chemical composition of the reactant gas at the inlet of the model is oxygen, nitrogen, and water, with the boundary condition being the average velocity inlet and the outlet set as a fully developed airflow to suppress reflux and normal flow. The other parameters validated are displayed in Table 1. Therefore, in the design of the dual channel serpentine flow field, the same control equations, conditions, and parameters as the reference paper were used, and the average depth model was compared with the various results of the reference 3D model. It is found that although the results of 2D PEMFC cathode simulation using average depth model are somewhat different from those of 3D model, their results can be used for the research of radial flow field, so they can be applied to the topology optimization of 2D model, as shown in Fig. 2.

Figure 2.

Verification of the average depth model for a dual channel serpentine flow field.

3.2 Boundary conditions and mesh design of radial flow topology optimization model

The pre-optimized radial flow field has a velocity of 2 m s⁻¹ at the inlet, a mass fraction of 0.228 for oxygen, 0.023 for water, and 0.749 for nitrogen. The edge boundary has a periodic cycle condition. The outlet is set to zero stress conditions (∇σn = 0) and zero flux ($\nabla \omega_{\text{O}_{2}}\text{n} = 0,\nabla \omega_{\text{H}_{2}\text{O}}\text{n} = 0$) for the transported species, as shown in Fig. 3a. The attributes and parameters used in modeling are shown in Table 1.

Figure 3.

Boundary condition setting and grid division of the cathode model for optimization.

COMSOL uses the Hermann von Helmholtz filter to extract only the grid information of the required element, and filters the material volume factor by the product of the Laplace operator of the minimum filter radius and the volume factor, which improves the calculation efficiency. Therefore, the denser the grid, the smaller the minimum unit of the model’s topological logic structure, and the clearer the topological logic structure model can be obtained. Therefore, the model adopts an extremely fine free quadrilateral fluid dynamics grid distribution, as shown in Fig. 3b.

3.2.1 Model optimization analysis

This section shows the results of topology optimization of the model, as shown in Fig. 4. Among them, the dark region is the porous rib plate part, and the light region is the porous channel part. The cathode channel and GDL region are stretched to form a 3D model, respectively. The parameters of the 3D model are shown in Table 1.

Figure 4.

Radial PEMFC topological structure.

The 2D model diagram of topology optimization is directly compared with its 3D model diagram through the finite element method, as shown in Fig. 5. The physical field distribution of the 2D model is basically consistent with that of the 3D model. According to its velocity field graph, the speed of the 3D model is 41 % higher than that of the 2D model, because the 2D model calculates its average velocity while the 3D model calculates its peak velocity. In the oxygen concentration distribution map and water concentration distribution map, the oxygen concentration distribution difference and water concentration distribution difference of the 2D model are far less than those of the 3D model, which is due to the fact that the channel and GDL of the 3D model are calculated separately, and the attributes of the channel and GDL in the 2D model are mixed and averaged, which reduces the gas diffusion, so the oxygen concentration distribution difference and water concentration distribution difference of the 2D model are significantly lower than those of the 3D model, Therefore, the 3D model formed by stretching the 2D model of its topology optimization can more accurately reflect the optimization effect.

Figure 5.

Comparison of topological structure PEMFC 2D and 3D model results.

Therefore, by designing PEMFC models with different inclinations and different shapes of obstacles, we can use topology optimization to design the transverse structure and optimize the longitudinal structure at the same time. Affected by the longitudinal depth of the model channel, the inclination of the channel is designed to be 1°, 2° and 3°. The shape of the obstacle structure in the longitudinal section is rectangular, triangular, and semicircular, with a depth of 5 × 10⁻⁴ m, as shown in Fig. 6.

Figure 6.

PEMFC model with different inclinations and structural obstacles.

3.2.2 Velocity analysis of PEMFC models with different inclinations and obstacles

It can be seen from the speed distribution diagram of each case in Fig. 7 that the speed of each case has been improved to varying degrees. The peak speed of Case 3 is 2.97 m s⁻¹, which is 3 % higher than Case 1 and Case 2. The peak speed of Case 5 and Case 6 is 2.99 m s⁻¹, which is 4.1 % higher than Case 1 and Case 2. It is found from the comparison of Cases 1, 2, and 3 that the peak speed increases with the increase of inclination, and the speed distribution of Case 3 is more uniform. In Cases 4, 5 and 6, the speed of reactants will be increased when passing through obstacles, which is more obvious in the speed distribution diagram of Case 5, and the peak speed of Cases 4, 5, and 6 is larger than that of Cases 1, 2, and 3, which indicates that setting obstacles is more effective for speed improvement than increasing the gradient.

Figure 7.

Velocity distribution diagram of radial PEMFC cathode models for each case.

3.2.3 Oxygen concentration analysis of PEMFC models with different inclinations and obstacles

From Fig. 8, it can be observed that Cases 2, 3, 4, and 6 have the largest decrease in oxygen concentration values, at 1.7 mol m⁻³, while Cases 1 and 5 have a smaller decrease in oxygen concentration of 1.5 mol m⁻³. According to the oxygen concentration distribution, it can be observed that the oxygen concentration distribution of Case 1, 2, and 3 is more uniform than that of Case 4, 5, and 6, which indicates that the setting of inclination is more effective than the setting of barriers for the uniform distribution of oxygen concentration. The oxygen concentration distribution maps of Cases 1, 2, and 3 show that the oxygen concentration distribution uniformity increases with the increase of inclination.

Figure 8.

Oxygen concentration distribution diagram of radial PEMFC cathode model for each case.

3.2.4 Pressure analysis of PEMFC models with different inclinations and obstacles

From Fig. 9, it can be observed that the maximum pressure difference for Case 3 is 13.1 Pa, while the minimum pressure difference for Cases 1 and 6 is 3.0 Pa. From the pressure distribution diagrams of Cases 1, 2, and 3, it can be seen that as the inclination increases, the pressure difference of the reaction gas gradually increases, and the pressure decreases more uniformly. From the pressure distribution diagrams of Cases 4, 5, and 6, it can be seen that the maximum pressure difference in Case 5 is 7.57 Pa, while the minimum pressure difference in Case 6 is 2.98 Pa. After setting obstacles, the pressure of the reaction gas decreases gradually after passing through each obstacle. Cases 3 and 5 will cause significant pressure loss and are not conducive to drainage, so the layout of Cases 1, 2, 4, and 6 is more suitable.

Figure 9.

Pressure distribution diagram of radial PEMFC cathode models for each case.

3.2.5 Polarization curve analysis of PEMFC models with different inclinations and obstacles

As shown in Fig. 10a, the average current density and voltage polarization curves of parallel flow field, serpentine flow field, topology optimization model and optimization case are given. It can be seen that the radial flow field optimized by topology has higher efficiency than traditional flow fields at medium to high current densities. As shown in Fig. 10b, under a working voltage of 0.6 V, the maximum average current density of Case 2 is 2434 A m⁻², which is 5 % higher than Case 3, 4 % higher than the serpentine flow field, 28 % higher than the parallel flow field, and 3.8 % higher than the TOM radial flow field; At a working voltage of 0.5 V, the maximum current density of Case 3 is 3181 A m⁻², which is 5 % higher than Case 2, 22 % higher than the serpentine flow field, and 6 % higher than the TOM radial flow field. At the same time, at a working voltage of 0.6 V, the average current density of Cases 1, 3, 4, and 5 is lower than the snake shaped flow field and TOM radial flow field. Case 4 has the worst effect, but at a working voltage of 0.5 V, all cases have better effects than the snake shaped flow field. However, the average current density of Cases 4 and 6 is lower than the TOM radial flow field. Therefore, increasing the inclination has a better effect on improving the average current density than setting obstacles, and at operating voltages of 0.5 V and 0.6 V, Case 2 has the best overall effect.

Figure 10.

Polarization curve of radial PEMFC cathode model in each case and average current density histogram of each case under operating voltages of 0.5 V and 0.6 V.

4. Conclusion

The cathode model of PEMFC radial flow field was designed through topology optimization. This model was constructed using a 2D average depth program, which synthesized the channel, GDL, and cathode electrode into a planar model. The rationality of this program was verified by comparing the application of 2D and 3D snake shaped flow channels with the average depth model. 3D modeling of the obtained TOM model and optimizing longitudinal depth by applying different tilt angles and obstacles. By comparing the velocity, oxygen concentration, and pressure distribution maps and polarization curves of each case model, the following conclusions can be drawn:

1. (1)    By comparing the speed distribution map of the optimization case, it is found that the application of obstacles has a significant optimization effect on the speed. The peak speed of Case 3 is 2.97 m s⁻¹, which is 3 % higher than that of Cases 1 and 2, and the peak speed of Cases 5 and 6 is 2.99 m s⁻¹, which is 4.1 % higher than that of Cases 1 and 2.
2. (2)    By comparing the oxygen concentration distribution map of the optimization case, it is found that the setting of the inclination has a better impact on the oxygen concentration distribution. With the increase of the inclination, the uniformity of the oxygen concentration distribution is also increased.
3. (3)    By comparing the pressure distribution maps of the optimized cases, as the inclination increases, the pressure difference of the reaction gas gradually increases, and the pressure decreases more evenly. After setting obstacles, the pressure of the reaction gas decreases gradually after passing through each obstacle. Cases 3 and 5 will cause significant pressure loss, which is not conducive to drainage. Therefore, the layout of Cases 1, 2, 4, and 6 is more suitable.
4. (4)    By comparing the polarization curve of the optimized case and the average current density of each case under working voltages of 0.5 V and 0.6 V, it was found that increasing the inclination has a better effect on improving the average current density than setting obstacles. Moreover, under working voltages of 0.5 V and 0.6 V, Case 2 has the best overall effect.

According to the set objective equation, algorithm optimization of channel structure can be achieved, and the optimal distribution case can be found in the design space of uniformly distributed materials. The constraint equation supports multiple design variables, making it easy for designers to optimize the battery performance of PEMFC. As an algorithm optimization, topology optimization can obtain the optimization results more scientifically and save the design time and calculation cost of researchers. Therefore, the application of topology optimization in PEMFC has great room for development.

CRediT Authorship Contribution Statement

Cheng Qu: Conceptualization (Lead), Data curation (Lead), Formal analysis (Lead), Methodology (Lead), Software (Lead), Validation (Lead), Writing – original draft (Lead), Writing – review & editing (Equal)

Minggang Zheng: Funding acquisition (Lead), Resources (Equal), Writing – review & editing (Equal)

Conflict of Interest

The authors have no relevant financial or non-financial interests to disclose.

Funding

Jinan ‘20 New Universities’ Introduction Innovation Team Project: 2021GXRC075

References

• 1)  Y. Wang, D. F. R. Diaz, K. S. Chen, Z. Wang, and X. C. Adroher, Mater. Today, 32, 178 (2020).
• 2)  F. Xiao, Y. C. Wang, Z. P. Wu, G. Y. Chen, F. Yang, S. Q. Zhu, K. Siddharth, Z. J. Kong, A. L. Lu, J. C. Li, C. J. Zhong, Z. Y. Zhou, and M. H. Shao, Adv. Mater., 33, 2006292 (2021).
• 3)  M. Marappan, K. Palaniswamy, T. Velumani, K. B. Chul, R. Velayutham, P. Shivakumar, and S. Sundaram, Chem. Rec., 21, 663 (2021).
• 4)  J. K. Liu and Y. S. Ma, Adv. Eng. Softw., 100, 161 (2016).
• 5)  A. Vogel and P. Junker, Struct. Multidiscipl. Optim., 63, 95 (2021).
• 6)  Z. Y. Zhang, Y. Zhao, B. X. Du, X. Q. Chen, and W. Yao, Struct. Multidiscipl. Optim., 62, 3071 (2020).
• 7)  J. P. Wang, J. Wu, and R. Westermann, Struct. Multidiscipl. Optim., 65, 92 (2022).
• 8)  P. Gomes and R. Palacios, Struct. Multidiscipl. Optim., 62, 2117 (2020).
• 9)  N. Pollini, Renew. Energy, 195, 1015 (2022).
• 10)  F. Mathieu-Potvin and L. Gosselin, Int. J. Hydrogen Energy, 39, 7382 (2014).
• 11)  Q. Xun and Y. J. Liu, Int. J. Energy Res., 43, 4807 (2019).
• 12)  C. Kim and H. Sun, Int. J. Automot. Technol., 13, 783 (2012).
• 13)  Z. M. Zhang, J. Zhang, and T. Zhang, Sustainability (Basel), 14, 4730 (2022).
• 14)  R. Behrou, A. Pizzolato, and A. Forner-Cuenca, Int. J. Heat Mass Transfer, 135, 72 (2019).
• 15)  T. Nakanishi, K. Koseki, S. Maruyama, and H. Shundo, Solid electrolyte fuel cell, US, US4910100 A, Mar 20 (1990).
• 16)  R. Diethelm, J. Brun, and B. Barp, Module for a fuel cell battery: US (1993).
• 17)  B. Emad, High temperature fuel cell, US, AU5455696A, Jan 2 (1997).
• 18)  S. Alexander, S. John, J. Markus, and H. Philipet, Fuel cell battery with integrated heat exchanger, EP, EP1347529B1, May 7 (2008).
• 19)  D. Bruno and K. Martin, Fuel cell battery with afterburning at the periphery of a cell stack, US, US6432567B1, Aug 13 (2002).
• 20)  D. Bruno and B. Emad, Fuel cell battery with a stack of planar cells, US, AU771489B2, Mar 25 (2004).
• 21)  S. Cano-Andrade, A. Hernandez-Guerrero, M. R. von Spakovsky, C. E. Damian-Ascencio, and J. C. Rubio-Arana, Energy, 35, 920 (2010).
• 22)  G. Vladimir, B. Frano, and N. K. Jay, Fuel cell collector plates with improved mass transfer channels, US, US6551736B1, Apr 22 (2003).
• 23)  I. Perez-Raya, A. Hernandez-Guerrero, F. Elizalde-Blancas, D. Juarez-Robles, and L. Almanza-Huerta, Proceedings of the ASME International Mechanical Engineering Congress and Exposition - 2010, Vol. 5, p. 899 (2012).
• 24)  B. R. Friess and M. Hoorfar, Int. J. Hydrogen Energy, 37, 7719 (2012).
• 25)  Z. H. Zhang and M. G. Zheng, Energy Source Part A, 44, 8266 (2022).
• 26)  Comsol multiphysics, v. 5.3a, comsol ab, Stockholm, Sweden (2017). URL: www.comsol.com
• 27)  C. Wang, Q. L. Zhang, S. Y. Shen, X. H. Yan, F. J. Zhu, X. J. Cheng, and J. L. Zhang, Sci. Rep., 7, 43447 (2017).
• 28)  L. H. Olesen, F. Okkels, and H. Bruus, Int. J. Numer. Methods Eng., 65, 975 (2006).
• 29)  R. J. Kee, M. E. Coltrin, and P. Glarborg, Chemically Reacting Flow: Theory and Practice, John Wiley & Sons (2005).
• 30)  J. Alexandersen, N. Aage, C. S. Andreasen, and O. Sigmund, Int. J. Numer. Methods Fluids, 76, 699 (2014).
• 31)  B. S. Lazarov and O. Sigmund, Int. J. Numer. Methods Eng., 86, 765 (2011).
• 32)  F. Wang, B. S. Lazarov, and O. Sigmund, Struct. Multidiscipl. Optim., 43, 767 (2011).