质子交换膜燃料电池(proton exchange membrane fuel cell, PEMFC)具有能量转化率高、零排放等优点，被认为是未来具有潜力的高效绿色发电装置^[1-2]。PEMFC运行过程中，燃料的输送和水管理是影响PEMFC性能的关键因素。质子交换膜含水量过低会导致膜电阻增大，但是在较高的电流密度下，PEMFC由于电化学反应加快，会产生较多的水导致流道和多孔结构的膜电极(membrane electrode assembly，MEA)出现水淹^[3]，致使PEMFC中燃料的质量传输受限，电池性能下降^[4-5]。

很多研究者建立了数值模型研究PEMFC内部的水传输和燃料传输情况。Springer等开发了一个等温一维模型，结果表明膜电阻随电流密度的增加而增大，较薄的质子交换膜在缓解膜电阻较高问题方面具有巨大优势^[6]；文[7]提出了一个包含铝板、集电器、气体通道、催化层(catalyst layer，CL)和质子交换膜的稳态一维模型，预测了电流密度和相对湿度对水输送系数的影响，表明输送系数较低甚至为负值的工况有助于避免膜的脱水。文[8]通过建立的一维模型，发现非均匀温度和压力对预测阳极和阴极扩散层(gas diffusion layer, GDL)中的蒸气和液态水通量有很大影响。J X Liu等建立了二维全电池分析模型，考虑了阴阳极催化层上的电化学反应、反应物的扩散效应和流道状态，研究了极化对PEMFC性能的影响，结果表明极化主要发生在阴极侧，但是阳极过电位也不可忽略^[9]。Chevalier等提出了一种无量纲准二维稳态电流密度分布模型，并发现沿流道分布的3个无量纲参数控制电流密度分布：流道和气体扩散层界面的Peclet数、气体扩散层和催化层界面的Damköhler数以及催化层和膜界面的Wagner数^[10]。Lei等用电渗阻力、反扩散和水力渗透的综合机理描述了水的运移，结果表明在较低的电压下，电流密度的迅速下降是由于通过质子交换膜的氧扩散阻力增加所致^[11]。

由于二维模型的局限，更多的研究者建立了三维模型研究PEMFC内部情况。单通道的三维多相PEMFC模型能预测工作压力、温度、空气化学计量和反应物相对湿度对PEMFC性能的影响^[12]。计算多孔扩散层的相对含水饱和度的扩散系数模型，能模拟PEMFC内的水淹情况^[13]。

流道的设计直接影响反应气体的分布和输送、水管理以及PEMFC的燃料利用率，最终影响电池的输出性能。强制对流具有优越的传质优势，可以促进燃料的输送。H C Liu等开发了一个模型来研究挡板以阵列形式排列在氧气传输通道中对电池性能的影响，研究表明隔板的合理布置可以提高PEMFC的局部电流密度^[14]。Su等分别设计了适用于平行流场和蛇形流场的阶梯流道，他们发现，阶梯流道对平行流场有明显影响，但对蛇形流场和扩散层传质性能影响不大^[15]。然而，目前障碍物流场数值模拟的模型为了计算简便，简略了一些影响PEMFC性能的因素，难以准确描述电池内部的真实情况。

本文采取模拟与实验相结合的研究手段，建立了与实验条件一致的直流道多相三维模型，利用实验的极化曲线对模拟结果进行拟合，确定模型有效性参数，在此基础上建立了波形流道进行模拟，分析了两种流道内氧气、液态水、速度以及电流密度分布，对比研究了直流道和波形流道对PEMFC性能提升的机理。此多相三维模型增加了以下因素：1) 考虑了阴阳极微孔层(micro porous layer，MPL)对PEMFC的影响；2) 增加了阴极粒子模型，考虑了催化剂微孔结构中的质量传递阻力，包括由聚合物薄膜引起的电阻以及催化剂中颗粒物质周围液态水膜导致的阻力；3) 考虑了MEA中水的产生和运输情况；4) 通过黏性动量阻力模拟压降，求解气体通道中的液体饱和度。

1 数学模型 1.1 守恒方程 1.1.1 质量守恒方程^[16-18]

$ \frac{\partial(\varepsilon \rho)}{\partial t}+\nabla \cdot(\varepsilon \rho \boldsymbol{u})=0. $

(1)

式中：ε为孔隙率；ρ为密度，kg/m³；u为速度，m/s；t为时间，s。

1.1.2 动量守恒方程^[16-18]

$ \frac{\partial(\varepsilon \rho \boldsymbol{u})}{\partial t}+\nabla \cdot(\varepsilon \rho \boldsymbol{u} \boldsymbol{u})=-\varepsilon \nabla p+\nabla \cdot(\varepsilon \mu \nabla \boldsymbol{u})+S_{\mathrm{u}}. $

(2)

式中：p为压强，Pa；μ为黏度，Pa·s；S_u为源项，在气体扩散层内可以简化为Darcy定律，

$ \varepsilon \boldsymbol{u}=\frac{K_{\mathrm{gdl}}}{\mu} \nabla p. $

(3)

其中K_gdl为GDL渗透率，m²。对于气体和液态两相流，

$ \mu=s_{1} \mu_{1}+\left(1-s_{1}\right) \mu_{\mathrm{g}} . $

(4)

其中：下标g代表气体，l代表液体；s是饱和度，为多孔材料中液态水的体积(V_l)与孔隙体积(V_p)的比值，

$ s_{1}=\frac{V_{1}}{V_{\mathrm{p}}}. $

(5)

式中V为体积，m³。

1.1.3 能量守恒方程^[16-18]

$ \begin{gathered} \frac{\partial\left(\varepsilon \cdot \rho c_{p} T\right)}{\partial t}+\nabla \cdot\left(\varepsilon \cdot \rho \cdot c_{p} \boldsymbol{u} T\right)= \\ \nabla \cdot\left(k^{\text {eff }} \nabla T\right)+S_{\mathrm{h}}. \end{gathered} $

(6)

其中：c_p为定压比热容，J/(kg·K)；T为温度，K；k^eff为有效导热系数；S_h为源项，在GDL、MPL、流道、CL和膜内的计算公式分别为

$ S_{\mathrm{h}}=\left\{\begin{array}{l} i_{\mathrm{m}}^{2} / \sigma_{\mathrm{sol}}(\mathrm{GDL}+\mathrm{MPL}), \\ 0(\text { 流道 }), \\ i_{\mathrm{s}}^{2} / \sigma_{\mathrm{sol}}-S_{\mathrm{gl}} \cdot L_{\mathrm{w}}(\text { 集电器 }), \\ R_{\mathrm{an}}\left(\eta_{\mathrm{an}}-T \Delta S_{\mathrm{an}} /(2 \mathrm{~F})\right)+i_{\mathrm{s}}^{2} / \sigma_{\mathrm{sol}}+i_{\mathrm{m}}^{2} / \sigma_{\mathrm{mem}}- \\ \ \ \ \ \left(S_{\mathrm{ld}}+S_{\mathrm{gl}}\right) L_{\mathrm{w}}(\text { 阳极 } \mathrm{CL}), \\ R_{\mathrm{cat}}\left(\eta_{\mathrm{cat}}-T \Delta S_{\mathrm{cat}} /(2 F)\right)+i_{\mathrm{s}}^{2} / \sigma_{\mathrm{sol}}+i_{\mathrm{m}}^{2} / \sigma_{\mathrm{mem}}- \\ \ \ \ \ \left(S_{\mathrm{ld}}+S_{\mathrm{gl}}\right) L_{\mathrm{w}}(\text { 阴极 } \mathrm{CL}), \\ i_{\mathrm{m}}^{2} \sigma_{\mathrm{mem}}(\text { 膜 }) . \end{array}\right. $

(7)

其中：i_s、i_m分别为固相和膜相的电流密度，A/cm²；L_w为由于水凝结引起的潜热，J/mol；σ_sol和σ_mem分别为固体电导率和膜的电导率，1/(Ω·m)，η_an、η_cat分别为阳极、阴极的表面过电位，V；S_gl为气相和溶解相的相互转化速率，S_ld为液相和溶解相的相互转化速率, kg/(m³·s)。

1.1.4 物料守恒方程

$ \frac{\partial\left(\varepsilon c_{i}\right)}{\partial t}+\nabla \cdot\left(\varepsilon \boldsymbol{u} c_{i}\right)=\nabla \cdot\left(D_{i}^{\text {eff }} \nabla_{c_{i}}\right)+S_{i}. $

(8)

其中c_i是物料i的质量分数。S_i为物料i的源项：

$ S_{\mathrm{H}_{2}}=-\frac{M_{\mathrm{w}}}{2 F} R_{\mathrm{an}}<0 , $

(9)

$ S_{\mathrm{O}_{2}}=-\frac{M_{\mathrm{w}}}{2 \mathrm{~F}} R_{\mathrm{cat}}<0 , $

(10)

$ S_{\mathrm{H}_{2} \mathrm{O}}=\frac{M_{\mathrm{w}}}{4 \mathrm{~F}} R_{\mathrm{cat}}>0. $

(11)

R_an、R_cat分别为阳极、阴极的交换电流密度，A/cm²，具体表达式见式(24)和(25)。D_i^eff为气体组分的扩散率，通过稀释近似方法计算，

$ D_{i}^{\mathrm{eff}}=\varepsilon^{1.5}(1-s)^{\gamma_{\mathrm{s}}} D_{i}^{0}\left(\frac{p_{0}}{p}\right)^{\gamma_{\mathrm{p}}}\left(\frac{T}{T_{0}}\right)^{\gamma_{\mathrm{t}}} \text {. } $

(12)

其中：D_i⁰是组分i在参考温度和压力下的质量扩散率^[18]；p₀为参考压强，Pa；T₀为参考温度，K。γ_s、γ_p、γ_t为特定参数，在模型中定义为：

$ \gamma_{\mathrm{s}}=1.0, \quad \gamma_{\mathrm{p}}=1.5, \quad \gamma_{\mathrm{t}}=2.0. $

(13)

膜内液态水传输的扩散系数为

$ D_{\mathrm{w}}=g(\lambda) \frac{\rho}{\mathrm{EW}} f(\lambda). $

(14)

式中EW为聚合物膜的当量，kg/kmol。f(λ) 和g(λ) 是含水量λ的函数：

$ f(\lambda)=4.1 \mathrm{e}^{-10}\left(\frac{\lambda}{25}\right)^{0.15}\left[1+\tanh \left(\frac{\lambda-2.5}{1.4}\right)\right] . $

(15)

$ g(\lambda)=2.5 \frac{\lambda}{22} . $

(16)

1.2 电化学反应

电化学反应背后的驱动力是表面过电位，即实际电动势和可逆电动势的差值，在此模型中体现为固相的相电位和膜的相电位的差异。其中：固相中主要体现的是电子在导电材料中的传输，包括催化剂层、多孔介质、双极板和集电器；膜内主要体现的是质子即氢离子在膜内的传输。

$ \nabla \cdot\left(\sigma_{\mathrm{sol}} \nabla \phi_{\mathrm{sol}}\right)=S_{\mathrm{sol}}, $

(17)

$ \nabla \cdot\left(\sigma_{\mathrm{mem}} \nabla \phi_{\mathrm{mem}}\right)=S_{\mathrm{mem}}. $

(18)

其中：ϕ_sol和ϕ_mem分别为固相电势和膜电势，V，

$ \sigma_{\text {mem }}=(0.514 \lambda-0.326)^{\omega_{i}} \mathrm{e}^{1\ 268\left(\frac{1}{303}-\frac{1}{T}\right)}. $

(19)

其中：ω_i为质子传导指数；λ的计算公式为

$ \lambda=0.043+17.18 a-39.85 a^{2}+36 a^{3} \quad(a<1), $

(20)

$ \lambda=14+1.4(a-1) \quad(a>1) . $

(21)

其中a是水活度，由水蒸气压强和饱和蒸气压强决定：

$ a=\frac{p_{\mathrm{wv}}}{p_{\mathrm{sat}}}, $

(22)

$ \begin{gathered} \lg p_{\mathrm{sat}}=-2.179\ 4+0.029\ 53(T-273.17)- \\ 9.183\ 7 \times 10^{-5}(T-273.17)^{2}+ \\ 1.445\ 4 \times 10^{-7}(T-273.17)^{3} . \end{gathered} $

(23)

其中：p_wv为水蒸气压强，p_sat为饱和蒸气压强，Pa。

电流守恒方程中的源项S_sol、S_mem具有以下定义：对于固相中的电势方程，阳极侧S_sol=-R_an(S_sol < 0)，阴极侧S_sol=+R_cat(S_sol>0)；对于膜相中的电势方程，阳极侧S_mem=+R_an(S_mem>0)，阴极侧S_mem=-R_cat(S_mem < 0)；R_an、R_cat具有以下定义^[19-20]：

$ R_{\mathrm{an}}=\zeta_{\mathrm{an}} j_{\mathrm{an}}^{\mathrm{ref}}\left(\frac{\left[\mathrm{H}_{2}\right]}{\left[\mathrm{H}_{2}\right]_{\mathrm{ref}}}\right)^{\gamma_{\mathrm{an}}}\left(\mathrm{e}^{\frac{a_{\mathrm{an}}^{\mathrm{an}} F \eta_{\mathrm{an}}}{R T}}-\mathrm{e}^{\frac{a_{\mathrm{cat}}^{\mathrm{an}} F_{\mathrm{an}}}{R T}}\right), $

(24)

$ R_{\mathrm{cat}}=\zeta_{\text {cat }} j_{\text {cat }}^{\mathrm{ref}}\left(\frac{\left[\mathrm{H}_{2}\right]}{\left[\mathrm{H}_{2}\right]_{\mathrm{ref}}}\right)^{\gamma_{\mathrm{cat}}}\left(\mathrm{e}^{\frac{\alpha_{\mathrm{a n}}^{\mathrm{cat}} F \eta_{\mathrm{cat}}}{R T}}-\mathrm{e}^{\frac{\alpha_{\mathrm{cat}}^{\mathrm{cat}} F \eta_{\mathrm{cat}}}{R T}}\right) . $

(25)

其中：ζ为比表面积，1/m；[ ]、[ ]_ref分别为物质的量浓度和参考条件下物质的量浓度。γ_an为阳极浓度系数，γ_cat为阴极浓度系数；α_an^an、α_cat^an分别为阳极电极上的阳极、阴极转移系数；α_an^cat、α_cat^cat分别为阴极电极上的阳极、阴极转移系数。j_an^ref和j_ref^cat分别为在指定温度下阳极、阴极的参考电流密度，A/m²；F为Faraday常数，C/mol；R为理想气体常数，J/(mol·K)。

表面过电位η的单位为V，计算公式如下：

$ \eta_{\mathrm{an}}=\phi_{\mathrm{sol}}-\phi_{\mathrm{mem}}-U_{\mathrm{an}}^{0}, $

(26)

$ \eta_{\mathrm{cat}}=\phi_{\mathrm{sol}}-\phi_{\mathrm{mem}}-U_{\mathrm{cat}}^{0}, $

(27)

U_an⁰、U_cat⁰分别为阴极和阳极处的半电池电位，V。

$ U_{\mathrm{an}}^{0}=E_{\mathrm{an}}^{0}-\frac{\Delta S_{\mathrm{an}}}{2 F}\left(T-T_{0}\right)-\frac{R T}{2 F} \ln \left(\frac{p_{\mathrm{H}_{2}}}{p_{0}}\right), $

(28)

$ U_{\mathrm{cat}}^{0}=E_{\mathrm{cat}}^{0}+\frac{\Delta S_{\mathrm{cat}}}{2 F}\left(T-T_{0}\right)-\frac{R T}{2 F} \ln \left(\frac{p_{\mathrm{H}_{2} \mathrm{O}}}{p_{\mathrm{sat}} \sqrt{p_{\mathrm{O}_{2}} / p_{0}}}\right). $

(29)

其中：E_an⁰、E_cat⁰分别为阳极、阴极的可逆电位，V；ΔS_an、ΔS_cat分别为阳极、阴极的反应熵，J/(kg·mol·K)。

1.3 阴极粒子模型

阴极层内采用阴极粒子模型，将离聚物薄膜引起的阻力和催化层内颗粒物周围液态水导致的阻力考虑在内，

$ R_{\mathrm{cat}}=4 F \frac{c_{\mathrm{O}_{2}}}{c_{\mathrm{O}_{2}} / j_{\mathrm{O}_{2}}^{\mathrm{ideal}}+R_{\mathrm{ion}}+R_{\mathrm{liq}}}. $

(30)

其中：c_O2是壁面的氧气摩尔浓度；R_ion是阴极离聚物电阻，R_liq为颗粒周围的液态水膜导致的阻力，s·m^-1，

$ R_{\text {liq }}=\frac{\zeta_{\text {cat }} r_{\mathrm{p}}^{2}}{K_{\mathrm{w}} D_{\mathrm{w}}} \cdot \frac{\sqrt[3]{1+s \varepsilon /(1-\varepsilon)}-1}{3(1-\varepsilon)}. $

(31)

其中：ζ_cat为阴极催化剂层比活性表面积，1/m；r_p催化剂颗粒半径，m；K_wD_w在液态水中的O₂溶解度和扩散率的乘积，m²·s^-1。

式(30)中j^ideal_O2计算公式为

$ j_{\mathrm{O}_{2}}^{\text {ideal }}=\frac{R_{\mathrm{cat}}}{4 F}. $

(32)

其中R_cat由方程(25)计算。

1.4 液态水输运方程

液态水在PEMFC中呈现气态、溶解相等，溶解相存在于催化剂层和膜中。采用式(33)描述溶解相在阴阳极的产生和运输，

$ \begin{array}{c} \frac{\partial}{\partial t}\left(\varepsilon_{\mathrm{i}} M_{\mathrm{w}} \frac{\rho}{\mathrm{EW}} \lambda\right)+\nabla \cdot\left(i_{\mathrm{m}} \frac{n_{\mathrm{d}}}{F} M_{\mathrm{w}}\right)= \\ \nabla \cdot\left(M_{\mathrm{w}} D_{\mathrm{w}} \nabla \lambda\right)+S_{\lambda}+S_{\mathrm{gd}}+S_{\mathrm{ld}}. \end{array} $

(33)

式中：i_m为离子电流密度，A/m²；n_d为渗透阻力系数；M_w为水的摩尔质量，g/mol；S_λ为催化剂层阴极侧反应水生成速率，S_gd为气体和溶解相之间的转化速率，kg/(m³·s)。S_gd、S_ld表达式如下：

$ S_{\mathrm{gd}}=\left(1-s^{\theta}\right) \gamma_{\mathrm{gd}} M_{\mathrm{w}, \mathrm{H}_{2} \mathrm{O}} \frac{\rho_{i}}{\mathrm{EW}}\left(\lambda_{\mathrm{eq}}-\lambda\right), $

(34)

$ S_{\mathrm{ld}}=s^{\theta} \gamma_{\mathrm{ld}} M_{\mathrm{w}, \mathrm{H}_{2} \mathrm{O}} \frac{\rho_{i}}{\mathrm{EW}}\left(\lambda_{\mathrm{eq}}-\lambda\right). $

(35)

其中：λ_eq为平衡含水量；θ为液体覆盖指数；γ_gd、γ_ld分别为气体和液体质量交换率常数。平衡含水量λ_eq计算公式为

$ \begin{gathered} \lambda_{\mathrm{eq}}=0.3+0.6 a(1-\tanh (a-0.5))+ \\ 0.69\left(\lambda_{a=1}-3.52\right) a^{0.5}\left(1+\tanh \left(\frac{a-0.89}{0.23}\right)\right)+ \\ s\left(\lambda_{s=1}-\lambda_{a=1}\right) . \end{gathered} $

(36)

λ_a=1为水活度为a=1时的平衡含水量，λ_s=1为水饱和度为s=1时的平衡含水量，均为无量纲量。

多孔电极和膜中液态水运输的驱动力是液体的压力梯度，具体体现为毛细力引发的压强p_c和气体压强p的作用：

$ \frac{\partial}{\partial t}\left(\varepsilon \rho_{1} s\right)=\nabla \cdot\left(\frac{\rho_{1} K K_{\mathrm{r}}}{\mu_{1}}\right) \nabla\left(p+p_{\mathrm{c}}\right)+S_{\mathrm{gl}}-S_{\mathrm{ld}}. $

(37)

式中：K为绝对渗透率，K_r为相对渗透率，m²；p为气体压强，p_c为毛细管压强，Pa。

在多孔的气体扩散层GDL和微孔层MPL中，相对渗透率为

$ K_{\mathrm{r}}=s^{d}. $

(38)

d为相对渗透率指数。

膜中的相对渗透率为

$ K_{\mathrm{r}}=\left(\frac{\frac{M_{\mathrm{w}}}{\rho_{1}} \lambda_{s=1}+\frac{\mathrm{EW}}{\rho_{i}}}{\frac{M_{\mathrm{w}}}{\rho_{1}} \lambda+\frac{\mathrm{EW}}{\rho_{i}}} \cdot \frac{\lambda}{\lambda_{s=1}}\right)^{2}. $

(39)

其中ρ_i为物种i的密度。

基于单向扩散理论计算气相和液相之间的转化速率S_gl，

$ S_{\mathrm{gl}}=\left\{\begin{array}{l} \gamma_{\mathrm{e}} \varepsilon s D_{\mathrm{gl}} \frac{M_{\mathrm{w}}}{R T} p \ln \left(\frac{p-p_{\mathrm{sat}}}{p-p_{\mathrm{wv}}}\right), p_{\mathrm{wv}} \leqslant p_{\mathrm{sat}} ;\\ \gamma_{\mathrm{c}} \varepsilon s D_{\mathrm{gl}} \frac{M_{\mathrm{w}}}{R T} p \ln \left(\frac{p-p_{\mathrm{sat}}}{p-p_{\mathrm{wv}}}\right), p_{\mathrm{wv}}>p_{\mathrm{sat}}. \end{array}\right. $

(40)

其中：γ_c和γ_e分别是冷凝率系数和蒸发率系数。D_gl有如下形式：

$ D_{\mathrm{gl}}=\left\{\begin{array}{l} 0.365 \times 10^{-4}\left(\frac{T}{343}\right)^{2.334} \times\left(\frac{10^{5}}{p}\right), \text { 阴极；} \\ 1.79 \times 10^{-4}\left(\frac{T}{343}\right)^{2.334} \times\left(\frac{10^{5}}{p}\right), \text { 阳极. } \end{array}\right. $

(41)

液态水在气体扩散层的毛细压力下进入气体通道，在内部不断聚集便会导致压降的存在，通过式(42)计算跟踪通道中的液态水，

$ \frac{\partial}{\partial t}\left(\rho_{1} s\right)+\nabla \cdot\left(\rho_{1} \boldsymbol{u}_{1} s\right)=\nabla \cdot\left(D_{\mathrm{liq}} \nabla s\right). $

(42)

其中：D_liq是气体通道中水扩散系数；u_l是液体速度，m/s，假定其为气体速度的一部分，

$ \boldsymbol{u}_{1}=\chi \boldsymbol{v}_{\mathrm{g}}. $

(43)

其中χ为液体与气体的速度比。

液态水在流道内的通量由毛细压力和气体通道中的动态压力驱动，

$ f_{\mathrm{liq}}=\varepsilon s \max \left[\left(p_{\mathrm{c}}+\frac{1}{2} \rho V^{2}\right), 0\right] . $

(44)

式(44)中p_c的计算公式为

$ p_{c}= \begin{cases}\psi \cos \theta_{c} \sqrt{\frac{\varepsilon}{K}} J(1-s), & \theta_{c}>90^{\circ} ; \\ \psi \cos \theta_{c} \sqrt{\frac{\varepsilon}{K}} J(s), & \theta_{c}<90^{\circ}.\end{cases} $

(45)

其中：ψ为表面张力，N/m；θ_c为接触角，(°), 当θ_c=90°时，饱合度s=0.5。J(x)为

$ J(x)=a x-b x^{2}+c x^{3}. $

(46)

其中：a，b，c分别为特定参数，在模型中定义为：

$ a=1.417, \quad b=2.12, \quad c=1.263 . $

(47)

2 模型验证

PEMFC涉及多种物理现象，包括电化学反应、传热、传质和相变等，内部反应的复杂性导致昂贵的实验费用，而数值研究成本较低，还能够得到实验无法得到的结果，因此被广泛应用到PEMFC的研究中。但是，数值研究必须在实验的标定和验证下才能证明数值仿真结果的正确性。

2.1 实验与仿真模型

实验采用平行直流道的PEMFC进行性能测试，反应面积为11.56 cm²(34 mm×34 mm)。平行流道上方的共用通道宽度设置为2 mm，流道宽度为0.8 mm，忽略共用通道后，采用L=30 mm的长度进行研究，如图 1a所示。采用15 μm的GORE膜；采用德国西格里SIGRACET 29BC的碳纸，厚度为0.235 μm；铂负载在阴极侧为0.4 mg/cm²，在阳极侧为0.2 mg/cm²；空气和氢气的化学计量比固定为2.5和2；阳极和阴极反应气体的相对湿度均为100%；电池温度为60 ℃；流道宽度、深度均为0.8 mm。建立与实验中所有尺寸一致的单根流道的计算区域模型，如图 1b所示。

实验测试装置由PEMFC单电池、电子负载、数据采集系统、加湿系统、氢气供气系统、空气供气系统等组成。采用的电子负载是由日本菊水公司生产的KFM2150和PLZ664WA。加湿采用鼓泡加湿。

2.2 模型有效性

为了验证模型的有效性，对所建立的单流道模型采取与实验测试条件一致的参数进行极化曲线的计算，包括化学计量比、电池温度、相对湿度、阴极侧通入空气、阳极侧通入氢气以及各自的化学计量比。结果表明，极化曲线和功率密度曲线的模拟结果与图 2中的实验结果吻合较好，证明了该模型的有效性。

3 结果与分析

本文设计了一种三维波形流道，如图 3所示，在节2的有效模型下进行相同模型参数(表 1)的计算。与常规直流道进行对比，波形流道内的波浪曲线为余弦曲线，流道内最高处高度为0.8 mm，最低处高度为0.4 mm，曲线方程为

$ y=0.3 \cos (2 x). $

(48)

计算区域包含的部件如图 3c所示。采用Ansys meshing将该区域划分为六面体单元，计算域截面示意图如图 3c右侧图所示。采用扫略划分网格，整个计算域共包含49万个网格单元。

3.1 电池性能曲线

图 4是直流道和三维波形流道的极化曲线以及功率密度曲线，可以看到三维波形流道比传统设计有更好的性能。在低电流密度下，极化曲线之间的差异较小，但在高电流密度下，极化曲线之间的差异更为明显。根据计算，在0.446 V的电压下，采用波形流道的PEMFC在峰值功率密度上比传统流道高10.16%。

3.2 氧气分布

图 5所示为不同流道(channel，CH)与GDL、MPL和CL 3个界面上的氧气平均质量分数随电流密度的变化。可以看出，氧气平均质量分数随着电流密度的增大而减小，这是由于反应速率增大导致需消耗更多的燃料。3个界面中CH-GDL界面上氧气平均质量分数最高，而GDL-MPL界面的氧气平均质量分数要小很多，这是因为气体扩散层虽然是多孔层，气体可以穿过它到达另一面，但是气体扩散层具有一定的厚度，导致气体在扩散过程中存在一定的阻力。MPL两侧的氧气浓度变化不大，这是由于MPL是传统GDL的附加疏水多孔层，能够改善PEMFC内的水管理和热量管理，而它的厚度较小，因此对氧气的传输影响较小。波形流道在较高的电流密度下，3个界面上的氧气浓度平均值都高于直流道，并且随着电流密度的增大，优越性增大，这是因为波形流道存在强制对流效应，就物质传输而言，对流效应比扩散效应有效得多。

图 6是电压为0.549 8 V时3个界面上的氧气质量分数。图 6a和6b是波形流道和直流道的CH-GDL界面处的氧气质量分数。随着气体进入流道，氧气质量分数越来越低，这是因为气体在狭小通道内会产生压降效应，这将导致到达催化层表面的气体浓度分布不均匀。在入口，两种流道的氧气质量分数差异不明显，但是在出口段，直流道中的平均氧气浓度明显低于波形流道，表明波形流道能够促进到达催化层表面气体浓度分布的均一性。在CH-GDL界面处，波形流道和直流道的氧气平均质量分数分别为0.179 59和0.169 24。

图 6c和6d是两种流道的GDL-MPL界面的氧气质量分数。与CH-GDL界面相比，氧气质量分数有了明显的降低，并且脊处的氧气质量分数低于流道下方。两种流道相比，波形流道狭窄区域内氧气质量分数出现突变，更加有利于气体从流道向气体扩散层的传输，这是由于凹凸的流道对气体产生强制对流效果，不仅使燃料更加容易透过气体扩散层，也有利于向脊处的扩散，从而进一步提高了PEMFC的性能。

图 6e和6f显示了MPL-CL界面的氧气质量分数分布，波形流道和直流道的氧气平均质量分数分别为0.087 03和0.099 51。波形流道将更多的燃料传输到催化层表面，增强了传质效率，提高了催化层的反应速率。

3.3 液态水去除性能分析

图 7是阴阳极流道内液态水的含量随电流密度的变化。结果表明，随着电流密度的增大，阴极流道内液态水的含量呈现增大的趋势，这是因为氧气在阴极催化层表面和氢离子以及电子反应生成水，而电流密度越大，生成的水就越多，当生成水量到一定程度时，水会突破气体扩散层进入流道。阳极处的液态水含量也呈现增加的趋势，这是由于阴极生成的水反扩散到阳极。流道内液态水的聚集会阻碍气体传输，影响PEMFC的性能。波形流道内聚集的液态水远少于直流道内的液态水，波形流道能够更加有效地将生成的水排出。

图 8是流道内液态水分布图。可以看出，液态水饱合度从入口处到出口处逐渐增加，这是因为入口处气体作用在液态水表面，液态水向前运动，会与流道中部以及后部的液态水聚集成大液滴，甚至是更大的贴壁大液滴，使得液态水运动更加困难。直流道内的液态水在流道中部和后部区域主要聚集在侧壁上，而波形流道对于侧壁上液态水的去除具有优越性，且在流道狭窄区域的液态水饱合度明显低于直流道。

3.4 气体速度与电流密度分布

图 9a展示了流道内气体速度分布。可以看到，波形流道内的气体速度在狭窄区域有较大的突变，这也揭示了波形流道的液态水去除优越性的原因。图 9b是Y轴方向的气体速度分布。在波形流道内的狭窄区域，气体具有垂直于气体扩散层的速度，并且在呈现周期性的波形流道内，速度也呈现周期性的突变，流道内Y轴方向的速度绝对值大于2 m/s。相比之下，直流道内垂直于GDL的方向几乎没有通量，流道内气体的浓度和催化层界面气体浓度的梯度产生的扩散效应是直流道GDL中传质的唯一途径。

氢在阳极催化层失去电子，失去电子的氢离子通过质子交换膜到达阴极催化层表面，电子则通过外电路到达阴极催化层，阴极催化层的氧气与氢离子和电子反应生水。随着电子的消耗，更多的电子通过外电路到达阴极侧，电子的定向移动形成电流。图 10是两种流道中的CL-MPL界面在电压为0.549 8 V下的电流分布。电化学反应主要发生在燃料流动的区域，在脊下方的催化层区域电流较小，这是因为在没有燃料气的地方，部分气体只能在浓度梯度的作用下扩散到流道脊下方的催化层，而波形流道中不规则的通道尺寸引发了强制对流效应，增强了燃料向催化层的传输，从而在阴极产生更多的电子，致使电子通过外电路流向阳极从而产生更多的电流。

4 结论

为了提高PEMFC的性能，本研究设计了一个横截面积在气体流动方向上周期性变化的波形流道, 并分析了电池性能提升的机理，主要的研究结果如下：

1) 建立三维、稳态、多相、等温、层流模型，通过实验标定模型，验证了模型的准确性后，设计了一种波形流道并在相同条件下进行计算，结果表明采用波形流道的PEMFC峰值功率密度提升10.16%。

2) 分析了各个组件界面的氧气分布，结果表明波形流道能够促进氧气向催化层的输送，增强了燃料气的强制对流效应。

3) 波形流道的狭窄区域提高了PEMFC流道内的气体流速，改善了PEMFC的水管理。