质子交换膜燃料电池（proton exchange membrane fuel cell，PEMFC）能够将化学能直接转化为电能，具有转换效率高、噪声低、无污染等优点^[1]. 气体扩散层（gas diffusion layer，GDL）是PEMFC的重要组成部分，位于双极板和催化层(catalytic layer，CL)之间，多由碳材料构成. GDL的孔隙率分布对反应气体传输和PEMFC输出性能均具有显著影响^[2-3]. 国内外学者为得到具有良好机械性能、优异导电性和合适孔隙率的GDL，进行了诸多有关结构设计和材料改性的研究. 刘舜等^[4]发现，当GDL孔隙率从0.2增加到0.8时，PEMFC的电流密度逐渐增大，最大可达847 mA/cm². Abraham等^[5]基于田口方法对GDL孔隙率进行优化，仿真结果表明GDL优化后，PEMFC的功率密度提高了12.5%. Xia等^[6]研究了GDL孔隙率大小和厚度对反应气体传输的影响，结果显示当阳极和阴极GDL的厚度分别为80~120 μm和140~170 μm时，氧气传输效果更佳；当GDL孔隙率为35%~45%时，电池性能更优. 程植源等^[7]将GDL沿厚度方向进行分层处理，认为具有梯度结构的GDL可以有效防止“水淹”并能有效提升燃料电池的性能. Jha等^[8]研究发现，当GDL孔隙率从0.6增加到0.7时，PEMFC的阴极氧气消耗减少了7%.

上述方法虽然对GDL的孔隙率大小及分布进行了优化，但孔隙率设计方案相对简单，不能充分展现孔隙率分布对燃料电池传质及输出性能的影响. 本研究分别建立平行流场、蛇形流场和叶脉流场的三维PEMFC模型，探究不同流场模型中常规扩散层孔隙率分布方案对氧气摩尔分数分布的影响规律；基于不同流场模型催化层中氧气摩尔分数的分布特征，提出对应3种流场模型分布的GDL孔隙率分布方案；分析不同PEMFC的膜电流密度分布、极化曲线与功率密度曲线.

碳布和碳纸是常见的2种GDL材料. 碳布质地柔软、支撑性差；碳纸多由碳纤维压制而成，具有良好的导电性和支撑力. 合适的孔隙能为PEMFC内反应物的传输提供通路，如图1所示为不同孔隙率的GDL材料SEM图像，其中ε_p为碳纸孔隙率，ε_c为碳布孔隙率.

如图2所示，为了分析平行、蛇形、叶脉流场中孔隙率大小和分布对PEMFC传质特性的影响，基于COMSOL Multiphysics建立包含阴阳极流场、气体扩散层、催化层及质子交换膜的三维PEMFC模型，整体模型长、宽均为20 mm. PEMFC模型的几何参数和操作条件如表1所示.

对PEMFC模型进行如下假设：1）PEMFC处于稳定工作状态，2）反应气体均为理想气体，3）反应气体流动为层流，4）反应气体不可压缩，5）忽略重力影响.

在PEMFC的电化学反应过程中，质量、动量、能量、电荷、组分均守恒.

1）质量守恒方程：

(1)$ \frac{{\partial \left( {\varepsilon \rho } \right)}}{{\partial t}}+\nabla \cdot \left( {\varepsilon \rho {\boldsymbol{u}}} \right) = {S _{\text{m}}} . $

式中：ε为多孔介质孔隙率，ρ为反应气体密度， t为时间，S_m为质量源项，u为反应气体速度矢量.

2）动量守恒方程：

(2)$ \frac{{\partial \left( {\varepsilon \rho {\boldsymbol{u}}} \right)}}{{\partial t}}+\nabla \cdot \left( {\varepsilon \rho {\boldsymbol{uu}}} \right) = - \varepsilon \nabla p+\nabla \cdot \left( {\varepsilon \mu \nabla {\boldsymbol{u}}} \right)+{S _{\text{u}}} . $

式中：$\;\mu $为反应气体动力黏性系数，p为压力，S_u为动量源项.

3）能量守恒方程：

(3)$ \frac{{\partial \left( {\varepsilon \rho {c _{{p}}}T} \right)}}{{\partial t}}+\nabla \cdot \left( {\varepsilon \rho {c _{{p}}}{\boldsymbol{u}}T} \right) = \nabla \cdot \left( {{k^{{\text{eff}}}}\nabla T} \right)+{S _{\text{Q}}} . $

式中：c_p为比定压热容，T为温度， k^eff为有效热导率，S_Q为能量源项.

4）电荷守恒方程：

(4)$ \nabla \cdot \left( {{\sigma _{{\text{sol}}}}\nabla {{\varphi }_{{\text{sol}}}}} \right)+{S_{{\rm{ sol}}}} = 0 , $

(5)$ \nabla \cdot \left( {{\sigma _{{\text{mem}}}}\nabla {{\varphi }_{{\text{mem}}}}} \right)+{S_{{\rm{ mem}}}} = 0 . $

式中：σ_sol为固相电导率，σ_mem为电解质相电导率，$\varphi_{\mathrm{sol}} $为固相电势，$\varphi_{\mathrm{mem}} $为电解质相电势， S_sol为固相电势源项，S_mem为电解质相电势源项.

5）组分守恒方程：

(6)$ \nabla \cdot (\varepsilon {\boldsymbol{u}}{{c}_{i}})=\nabla \cdot (D_{i}^{\text{eff}}\nabla {{c}_{i}})+{{S} _{i}}. $

式中：c_i为组分i的浓度； $D_i^{\mathrm{eff}} $为气体组分的有效扩散系数；S_i为组分源项. 组分源项在流道、扩散层和微孔层中为0.

一般使用Butler-Volmer方程描述催化层发生的反应，表达式为

(7)$ {i_{\text{a}}} = {i_{{\mathrm{ref,a}}}}{\left(\frac{{{c_{{{\text{H}}_{\text{2}}}}}}}{{c_{{{\text{H}}_{\text{2}}}}^{{\text{ref}}}}}\right)^{0.5}}\cdot\left(\exp \frac{{{\alpha _{\text{a}}}F{\eta _{\text{a}}}}}{{RT}} - \exp \frac{{ - {\alpha _{\text{c}}}F{\eta _{\text{a}}}}}{{RT}}\right) , $

(8)$ {i_{\text{c}}} = {i_{{\text{ref,c}}}}\left(\frac{{{c_{{{\text{O}}_{\text{2}}}}}}}{{c_{{{\text{O}}_{\text{2}}}}^{{\text{ref}}}}}\right)\cdot\left( - \exp \frac{{{\alpha _{\text{a}}}F{\eta _{\text{c}}}}}{{RT}}+\exp \frac{{ - {\alpha _{\text{c}}}F{\eta _{\text{c}}}}}{{RT}}\right) . $

式中：η_a为阳极活化过电压，η_c为阴极活化过电压，F为法拉第常数，R为摩尔气体常数，$ {i_{{\text{ref,a}}}} $为阳极体积参考交换电流密度，${{c}_{{{\text{H}}_{\text{2}}}}} $为氢气的浓度，$c_{{{\text{H}}_{2}}}^{\text{ref}} $为氢气的参考浓度，${{c}_{{{\text{O}}_{\text{2}}}}} $为氧气的浓度，$c_{{{\text{O}}_{2}}}^{\text{ref}} $为氧气的参考浓度，$ {i_{{\text{ref,c}}}} $为阴极体积参考交换电流密度，i_a为阳极交换电流密度，i_c为阴极交换电流密度.

为了验证网格无关性，对平行流场PEMFC模型采用4种密集程度不同的网格划分方案，在相同操作条件下对模型极化曲线进行仿真，并基于模型 4的电流密度计算相对误差$\delta $，结果如图3所示，其中N为网格数量. 仿真结果表明，当工作电压U_w=0.54 V时，4种网格划分方案的电流密度误差均满足计算要求，且随着网格数量增多，相对误差呈先增大后减小的趋势. 在4种模型中，模型3相对误差最小为3.51%，模型4相对计算时间较长，综合考虑计算精度和计算时间，采用模型 3对PEMFC模型进行数值模拟.

为了充分验证本研究建立模仿真方法的可行性，根据不同文献中的PEMFC几何参数及操作条件，进行3组模型验证. 第1组PEMFC模型是根据文献[9]中测试用的PEMFC几何尺寸建立的直流道模型；第2组PEMFC模型是根据文献[10]中测试用PEMFC几何尺寸建立的平行流场模型；第3组PEMFC模型建立与文献[11]尺寸相同的PEMFC流道几何模型，并与文献中的实验条件保持一致对PEMFC模型进行极化曲线仿真. 如图4所示，3组仿真极化曲线分别与文献中的实验极化曲线相吻合，验证了本研究几何模型建立与仿真计算方法的可行性. 此外，本研究采用的几何模型参数与操作参数与第3组保持一致.

为了分析扩散层孔隙率大小对平行流场PEMFC输出性能的影响，分别建立孔隙率为0.40、0.50和0.60的3种扩散层模型. 如图5所示为平行流场PEMFC在U_w=0.54 V、采用3种不同孔隙率(ε₁=0.40、ε₂=0.50、ε₃=0.60)分布方案时，催化层中的氧气摩尔分数$x_{{\mathrm{O}}_2} $分布. 可以看出，对应流道中部的催化层内氧气摩尔分数较低，氧气供应明显不足，而边缘处氧气摩尔分数相对充足. 3种不同孔隙率方案均因参与化学反应消耗使氧气摩尔分数不断减小，但相较于方案1、2，方案3的氧气摩尔分数分布更加均匀，且平均值相对较大.

孔隙率对扩散层中的气体传输起着至关重要作用. 孔隙率越大催化层界面上氧气摩尔分数越高，但当孔隙率过高时，GDL电导率将急剧下降，导致欧姆极化损失增大^[12]. 为了研究GDL孔隙率对PEMFC传质特性及输出性能的影响，基于图5催化层中氧气摩尔分数分布特征提出3种GDL孔隙率梯度分布方案，如图6所示. 方案4~6的孔隙率分布比方案1~3的孔隙率分布更规则，氧气摩尔分数高则孔隙率相对较小，氧气摩尔分数低则孔隙率相对较大，有利于氧气从流道到扩散层的均匀传质. 如图7所示，方案1~6中催化层氧气平均摩尔分数分别为0.017 5、0.021 3、0.022 1、0.020 5、0.021 2和0.021 6，其中方案3的氧气平均摩尔分数最大. 大孔隙率的GDL虽然提高了氧气向CL的扩散效率，但对比图1(b)和图1(c)可知，碳纤维体积减小会导致GDL导电性下降，不利于PEMFC电流密度的提高^[13].

如图8所示为平行流场PEMFC在U_w=0.54 V、采用6种不同孔隙率分布方案时，质子交换膜水平方向中间截面的膜电流密度J分布，平均膜电流密度J_eme分别为0.907 5、0.945 5、0.928 1、0.956 1、0.973 4和0.984 4 A/cm²，如图9所示. 由此可见，基于催化层中氧气摩尔分数分布特征提出的孔隙率梯度分布方案（方案5、6）能够提高平均膜电流密度^[13].

如图10所示为采用6种不同孔隙率分布方案时，平行流场PEMFC的极化曲线与功率密度曲线，峰值功率密度P分别为0.491 2、0.512 3、0.503 2、0.517 9、0.527 4和0.533 4 W/cm²，当PEMFC终止电压U=0.39 V时电流密度分别为1.058、1.119、1.101、1.139、1.148和1.161 A/cm². 方案6的峰值功率密度较前5种方案分别提高了8.59%、4.12%、6.00%、2.99%和1.14%，电流密度较前5种方案分别提高了9.74%、3.75%、5.45%、1.93%和1.13%.

增大GDL孔隙率能够促进氧气扩散，提高催化层中氧气平均摩尔分数，但是当孔隙率过高时，由于GDL电导率急剧下降，欧姆极化损失的增加值大于浓差极化损失的减小值，反而使电池性能下降（方案1~3）. 基于平行流场PEMFC催化层中氧气摩尔分数分布特征提出的GDL孔隙率梯度分布方案可以明显提高膜电流密度^[13]，改善电池输出性能（方案4~6）. 方案6在平均孔隙率不变的前提下，通过增大催化层中氧气供应不足处的孔隙率改善了局部氧气供应，催化层中氧气供应相对充足的位置维持较低的孔隙率使扩散层电导率相对较高，因此方案6的峰值功率密度和电流密度均最优^[14].

为了分析扩散层孔隙率大小对蛇形流场PEMFC输出性能的影响，分别建立孔隙率为0.40、0.50和0.60的3种扩散层模型. 如图11所示为催化层中氧气摩尔分数分布特征，相应的孔隙率（ε₁=0.40、ε₂=0.50、ε₃=0.60、ε₄=0.70、ε₅=0.30、ε₆=0.45）分布方案如图12所示. 由图11、12可以看出，蛇形流场催化层入口处氧气摩尔分数较出口处更高，沿气体流动方向氧气摩尔分数梯度下降. 方案7~14中的氧气平均摩尔分数分别为0.017 5、0.021 6、0.022 8、0.025 6、0.019 6、0.031 0、0.023 4和0.024 0，方案6催化层的平均氧气摩尔分数分布最佳.

如图13所示为蛇形流场PEMFC在U_w=0.54 V、采用8种不同孔隙率分布方案时，质子交换膜水平方向中间截面的膜电流密度分布. 蛇形流场膜电流密度沿气体流动方向均匀下降，流道出口处由于氧气供应明显不足导致电流密度过低. 8种分布方案的平均膜电流密度分别为0.982 8、1.029 9、1.008 7、1.116 7、1.041 3、1.160 1、1.098 7和1.105 1 A/cm²，如图14所示，其中方案12的平均膜电流密度相对较高.

如图15所示为采用8种孔隙率分布方案时蛇形流场PEMFC的极化曲线与功率密度曲线，峰值功率密度分别为0.531 9、0.558 1、0.547 0、0.605 4、0.563 8、0.629 0、0.595 6和0.599 1 W/cm²，终止电压为0.39 V时电流密度分别为1.122、1.192、1.185、1.270、1.238、1.270、1.270和1.270 A/cm². 方案12的峰值功率密度较其他7种分别提高了18.26%、12.70%、14.90%、3.90%、11.56%、3.90%、5.61%和4.99%. 对于蛇形流场PEMFC，增大GDL孔隙率可以提高催化层中氧气平均摩尔分数，但是当孔隙率过高时，由于GDL电导率急剧下降，反而使膜电流密度降低（方案7~9）. 沿水平方向增加GDL孔隙率层数可在一定程度上提高催化层中氧气平均摩尔分数和膜电流密度，但随着层数的增加效果变得不明显（方案8、10、13和14）. 方案12增大了蛇形流场出口附近的孔隙率，明显改善了氧气供应，而入口处采用较优孔隙率，因此方案12中的峰值功率密度和电流密度均最优.

为了分析扩散层孔隙率大小对叶脉流场PEMFC输出性能的影响，分别建立孔隙率为0.40、0.50和0.60的3种扩散层模型. 如图16所示为催化层中氧气摩尔分数分布特征，相应的孔隙率梯度分布方案如图17所示. 由图16可以看出，与主流道对应的催化层中氧气摩尔分数较高且沿着气体流动方向逐渐降低，整个催化层中的氧气摩尔分数沿主流道对称分布，氧气摩尔分数最低位置出现在气体出口附近的几条支流道中. 由图16、17可知，6种分布方案的CL氧气摩尔分数分布均因参与化学反应消耗而使其摩尔分数不断减小，但方案20相较于前5种分布方案氧气摩尔分数分布更加均匀且平均值相对较大，原因是方案20的孔隙率分布更加规则，有利于氧气从流道到催化层的单方向均匀传质，也说明孔隙率梯度分布的阴极扩散层气体传质优于单一分布的阴极扩散层气体传质. 如图18所示，6种方案的催化层氧气平均摩尔分数分别为0.023 5、0.029 1、0.028 2、0.028 3、0.028 0和0.043 7，其中方案20催化层中氧气平均摩尔分数最高.

如图19所示为叶脉流场PEMFC在U_w=0.54 V、采用6种孔隙率分布方案时，质子交换膜水平方向中间截面的膜电流密度分布，其中方案20的高膜电流密度区域面积最大. 6种分布方案的平均膜电流密度分别为1.059 0、1.104 4和1.064 3、1.116 4、1.120 2和1.215 1 A/cm²，如图20所示，方案20的平均膜电流密度远高于其他5种方案.

如图21所示为采用6种不同孔隙率分布方案时叶脉流场PEMFC的极化曲线与功率密度曲线，峰值功率密度分别为0.572 6、0.597 8、0.576 2、0.604 2、0.606 9和0.661 1 W/cm²，当终止电压为0.39 V时，电流密度分别为1.259、1.342、1.285、1.328、1.356和1.367 A/cm². 由图可知，当孔隙率采用方案20时，PEMFC输出性能最佳，峰值功率密度较前5种方案分别提高了15.46%、10.59%、14.73%、9.42%和8.93%，电流密度较前5种方案分别提高了8.58%、1.86%、6.39%、2.94%、0.81%. 对于叶脉流场PEMFC，增大GDL孔隙率不能改善催化层中氧气平均摩尔分数以及平均膜电流密度（方案15~17）. 基于叶脉流场PEMFC催化层中氧气摩尔分数分布特征提出的GDL孔隙率梯度分布方案可以明显提高膜电流密度，改善电池输出性能（方案18~20）. 方案20通过增大气体出口附近对应的GDL孔隙率改善了局部氧气供应，而主流道相对应的GDL孔隙率较低使扩散层电导率相对较高，因此方案20的峰值功率密度和电流密度均最优.

为了研究扩散层孔隙率对传质特性的影响，构建平行、蛇形和叶脉3种不同流场GDL模型，通过调控孔隙率大小和分布，考察催化层氧气摩尔分数、膜电流密度以及极化曲线和功率密度曲线变化规律，得到如下结论. 1）针对平行流场提出的最优GDL孔隙率分布方案，通过增大氧气供应不足处的孔隙率，改善了局部氧气供应不足现象，减小了GDL中低氧气摩尔分数区域的面积，能够有效改善电池输出性能. 2）沿水平方向增加蛇形流场GDL的孔隙率层数可提高催化层中氧气平均摩尔分数和膜电流密度；最佳孔隙率分布方案提高了蛇形流场出口区域的氧气供应，能有效提高电池的电流密度. 3）通过增大叶脉流场气体出口附近对应的GDL孔隙率，改善了催化层中氧气摩尔分数分布均匀性；与常规孔隙率分布相比，该方案的峰值功率密度和电流密度均最优.