0 引 言
质子交换膜燃料电池（proton exchange membrane fuel cell, PEMFC）具有能量转化效率高、启动快、零排放等优点,应用前景广阔[1]。然而,PEMFC的规模化使用仍然面临诸多挑战。例如,电池自零下温度环境启动时,可能会因结冰导致反应位点被覆盖或者反应（产）物传输通道被堵塞而无法启动。
数值模型是研究PEMFC内多场耦合过程的重要手段和工具。MAO等[2]首先建立了PEMFC的三维冷启动模型,研究了冷启动时电流密度、含水量和冰的三维分布特征;JIANG等[3]为了解PEMFC冷启动时非等温特性,研究了因温度升高导致电-热-气（水）传输物性的改变及其对电池冷启动性能的影响。尽管宏观尺度冷启动模型[3,4,5]已经实现PEMFC冷启动过程的仿真,但是尚不能够对结冰形貌等进行动态可视化,宏观模型中催化层（catalyst layer, CL）内含冰量对物质输运以及电化学有效反应面积的影响关系式都是基于纯粹的经验关系式,缺乏必要的微观基础。相对于宏观均质模型,介孔模型基于电极的介孔微结构建模[6,7],模型结果可揭示电极的构效关系,可用于设计或优化电极微孔结构,还可获得更准确的有效传输物性等宏观参数（或关系式）[8]。WANG等[9]首次利用理想规则结构的催化层建立了PEMFC过程介孔尺度数值模型,并研究了催化层形貌（孔径、厚度和孔隙率）对电池性能的影响,随后模型进一步发展为基于随机催化层微结构[10]。
本文基于随机网格法重建催化层介孔结构,针对PEMFC冷启动时催化层内发生的电-热-气（水）传输、电化学反应和结冰的多场耦合过程,建立了介孔尺度模型。研究了冷启动过程中,冰的形成和分布以及冰的形成对传质及电化学反应的影响,详细分析了PEMFC冷启动行为特性。
1 阴极催化层微结构重建
Table 1 Composition of the catalyst layer表 1 催化层的各相组成 |
| Component | Phase | Volume fraction / % |
|---|---|---|
| Carbon black | MS | 34.72 |
| Pt | MS | 1.24 |
| Ionomer | MS | 24.04 |
| Pore | P | 40.00 |
Fig. 1 SEM micrograph of cathode catalyst layer图1 阴极催化层的SEM图 |
Fig. 2 Result from the SGM reconstruction (red for mixed solid and light blue for the pore)图 2 重建结果（红色为固体混合相,淡蓝色为孔相） |
Fig. 3 Evolution of volume fraction of each phase along the direction of electrode thickness图3 各相体积分数沿电极厚度方向的变化 |
Fig. 4 Pore diameter distribution for the reconstructed CL and compared with the experimental result图4 重建催化层的孔径分布与实验结果的对比 |
2 数值模型
2.1 模型假设
阴极催化层位于质子交换膜（proton exchange membrane, PEM）和气体扩散层（gas diffusion layer, GDL）之间（如图5示）,x = 0处为催化层与气体扩散层交界面,x = xL（xL = 15 μm）处为催化层与PEM交界面;GDL为氧气、水蒸气和电子传输提供通道,PEM为质子和膜态水传输提供输运路径。从GDL侧进入的氧气与阳极反应产生的电子（经外电路后由GDL侧进入）和质子（通过PEM侧进入）在催化层内相遇,在催化剂作用下在活性位点（相界面）发生电化学反应。产物水一部分以水蒸气形式向GDL侧扩散;另一部分以膜态水的形式向PEM侧扩散（回扩散）。此外,电渗作用会使阳极及PEM内部分膜态水迁移到阴极催化层内。当催化层内的水（水蒸气）出现局部过饱和时,将发生相变而形成冰。
Fig. 5 Schematic of charge and species transfer in cathode catalyst layer图5 阴极催化层模型的电荷及物质输运示意图 |
2.2 控制方程
模型求解5个控制方程,分别是氧气浓度、水蒸气浓度、膜含水量方程及质子、电子负荷传输方程。
孔相中的氧气和水蒸气传输方程为:
$\frac{\partial \left( \varepsilon {{C}^{i}} \right)}{\partial t}\text{+}\nabla \cdot \left( \vec{u}{{C}^{i}} \right)=\nabla \cdot \left( D_{\text{eff}}^{i}\nabla {{C}^{i}} \right)+S_{\text{c}}^{i}$ (1)
式中:上标i代表组分（氧气或水蒸气）;C为组分浓度;t为时间;$\vec{u}$为速度矢量;S为对应源项。组分的有效扩散系数$D_{\text{eff}}^{i}$由Bruggeman关系给出:
$D_{\text{eff}}^{i}={{D}^{i}}{{\left( 1-{{{{\varepsilon }'}}_{\text{ice}}} \right)}^{b}}$ (2)
其中:${{{\varepsilon }'}_{\text{ice}}}$为孔相数值单元内的含冰量;b为Bruggeman因子,b = 1.5。物质扩散系数Di可由该组分的分子扩散系数$D_{\text{b}}^{i}$和努森扩散系数$D_{\text{k}}^{i}$计算:
$\frac{1}{{{D}^{i}}}\text{=}\frac{1}{D_{\text{b}}^{i}}\text{+}\frac{1}{D_{\text{k}}^{i}}$ (3)
其中:$D_{\text{k}}^{i}(\text{c}{{\text{m}}^{\text{2}}}\text{/s})=4850\overline{d}\sqrt{{T}/{{{M}^{i}}}\;}$;$\overline{d}$为平均孔径;T为温度;${{M}^{i}}$为物质的摩尔质量。
对于氧气,$S_{\text{c}}^{{{\text{O}}_{\text{2}}}}$代表了电化学反应的氧气消耗量,其大小与阴极传输电流密度jc的关系为:
$S_{\text{c}}^{{{\text{O}}_{\text{2}}}}=-\frac{{{s}^{{{\text{O}}_{\text{2}}}}}{{j}_{\text{c}}}}{nF}$ (4)
其中,${{s}^{{{\text{O}}_{\text{2}}}}}$和n分别为该组分的化学计量数和参与电极反应时的转移电子数。
针对离子聚合物内的水传输,膜含水量方程为:
${{\varepsilon }_{\text{m}}}\frac{\partial C_{\text{m}}^{{{\text{H}}_{\text{2}}}\text{O}}}{\partial t}=\nabla \cdot \left( D_{\text{m,eff}}^{{{\text{H}}_{\text{2}}}\text{O}}\nabla C_{\text{m}}^{{{\text{H}}_{\text{2}}}\text{O}} \right)\text{+}S_{\text{m,c}}^{{{\text{H}}_{\text{2}}}\text{O}}$ (5)
其中:${{\varepsilon }_{\text{m}}}$为离子聚合物在固体混合相内的体积分数;而$C_{\text{m}}^{{{\text{H}}_{\text{2}}}\text{O}}={{\rho }_{\text{m}}}\lambda /\text{EW}$,且
$D_{\text{m,eff}}^{{{\text{H}}_{\text{2}}}\text{O}}\text{=}\varepsilon _{\text{m}}^{1.5}\left\{ \begin{align} & 5.93\times {{10}^{-5}}\lambda \left( {{\text{e}}^{0.28\lambda }}-1 \right){{\text{e}}^{\left( -4269/T \right)}}\text{ } \\ & \text{ if 0}\lambda \le \text{3} \\ & 7.97\times {{10}^{-6}}\lambda \left( 1+161{{e}^{-\lambda }} \right){{\text{e}}^{\left( -4269/T \right)}} \\ & \text{ otherwise} \\ \end{align} \right.$ (6)
离子聚合物内与孔隙中的水浓度可用Springer关系式[18]描述:
λ = 0.043 + 17.18a - 39.85a2 + 36.0a3 0﹤a≤ 1 (7)
上式中,水活度$a=C_{{}}^{{{\text{H}}_{\text{2}}}\text{O}}/C_{\text{sat}}^{{{\text{H}}_{\text{2}}}\text{O}}$。
对于水传输方程,其源项$S_{\text{c}}^{{{\text{H}}_{\text{2}}}\text{O}}$不仅包括水的产生,还包括电渗和相变过程引起的水的增加或减少,对于孔相和固体混合相中的水传输方程,其源项分别为:
$S_{\text{g,c}}^{{{\text{H}}_{\text{2}}}\text{O}}=-\frac{{{s}^{{{\text{H}}_{\text{2}}}\text{O}}}{{j}_{\text{c}}}}{nF}-{{S}_{\text{ice}}}\text{+}{{S}_{\text{balance}}}$ (8)
$S_{\text{m,c}}^{{{\text{H}}_{\text{2}}}\text{O}}=-\nabla \cdot \left( \frac{{{n}_{\text{d}}}}{F}{{\mathsf{i}}_{\text{e}}} \right)-{{S}_{\text{balance}}}$ (9)
其中:nd为电渗系数;ie为质子电流;Sbalance为水蒸气和含水量方程在相界面处因热力学平衡过程所导致的源项。
${{S}_{\text{balance}}}=\frac{1}{V}\int\limits_{{{A}_{\text{inter}}}}{D_{\text{m},\text{g}}^{{{\text{H}}_{\text{2}}}\text{O}}\nabla {{C}^{{{\text{H}}_{\text{2}}}\text{O}}}\text{d}A}$ (10)
其中:Ainter为界面面积;V为控制体体积;$D_{\text{m},\text{g}}^{{{\text{H}}_{\text{2}}}\text{O}}$为固体混合相与孔相水传输系数的调和平均数。取决于局部热力学状态,孔相中部分的水蒸气可能会瞬间凝华成冰,Sice可写成:
${{S}_{\text{ice}}}=\left\{ \begin{align} & 0\text{ if }{{C}^{{{\text{H}}_{\text{2}}}\text{O}}}\le C_{\text{sat}}^{{{\text{H}}_{\text{2}}}\text{O}} \\ & {{r}_{\text{freeze}}}\left( C_{\text{sat}}^{{{\text{H}}_{\text{2}}}\text{O}}-{{C}^{{{\text{H}}_{\text{2}}}\text{O}}} \right)\text{ } \\ & \text{ if }{{C}^{{{\text{H}}_{\text{2}}}\text{O}}}>C_{\text{sat}}^{{{\text{H}}_{\text{2}}}\text{O}}\text{ and }T<T_{\text{freeze}}^{{{\text{H}}_{\text{2}}}\text{O}} \\ \end{align} \right.$ (11)
其中:rfreeze为结冰速率;$C_{\text{sat}}^{{{\text{H}}_{\text{2}}}\text{O}}$为水蒸气的饱和浓度;$T_{\text{freeze}}^{{{\text{H}}_{\text{2}}}\text{O}}$为冰点温度。
当冰形成时,Sbalance的热力学平衡面积部分被冰覆盖,假设被覆盖后面积${{A}_{\text{inter}}}$与${{{\varepsilon }'}_{\text{ice}}}$的关系为:
${{A}_{\text{inter}}}={{A}_{0,\operatorname{int}\text{er}}}\left( 1-{{{{\varepsilon }'}}_{\text{ice}}} \right)$ (12)
其中,${{A}_{0,\operatorname{int}\text{er}}}$为无冰时界面的面积。
电荷方程分为质子势ϕe方程（13）与电子势ϕs方程（14）:
$0=\nabla \cdot \left( \kappa _{\text{eff}}^{\text{e}}\nabla {{\phi }_{\text{e}}} \right)+{{j}_{\text{c}}}$ (13)
$0=\nabla \cdot \left( \sigma \nabla {{\phi }_{\text{s}}} \right)-{{j}_{\text{c}}}$ (14)
式中:σ为电子电导率;$\kappa _{\text{eff}}^{\text{e}}$代表质子传导率。MAO等[2]给出了${{\kappa }^{\text{e}}}$（S/cm）与λ的关系:
${{\kappa }^{\text{e}}}=\exp \left[ 2222\left( \frac{1}{303}-\frac{1}{{{T}_{\text{a}}}} \right) \right]\left( 0.005139\lambda -0.00326 \right)$(15)
模型的主控方程与源项归纳总结于表2中。
Table 2 Governing equations and source terms of the model表2 模型主控方程及源项 |
| Categories | Governing equations | Source terms | |
|---|---|---|---|
| P | MS | ||
| Species | $\frac{\partial \left( \varepsilon {{C}^{{{\text{O}}_{\text{2}}}}} \right)}{\partial t}=\nabla \cdot \left( D_{\text{eff}}^{{{\text{O}}_{\text{2}}}}\nabla {{C}^{{{\text{O}}_{\text{2}}}}} \right)+S_{\text{c}}^{{{\text{O}}_{\text{2}}}}$ | $S_{\text{c}}^{{{\text{O}}_{\text{2}}}}=-\frac{{{s}^{{{\text{O}}_{\text{2}}}}}{{j}_{\text{c}}}}{nF}$ | / |
| $\frac{\partial \left( \varepsilon {{C}^{{{\text{H}}_{\text{2}}}\text{O}}} \right)}{\partial t}=\nabla \cdot \left( D_{\text{eff}}^{{{\text{H}}_{\text{2}}}\text{O}}\nabla {{C}^{{{\text{H}}_{\text{2}}}\text{O}}} \right)+S_{\text{c}}^{{{\text{H}}_{\text{2}}}\text{O}}$ | $S_{\text{c}}^{{{\text{H}}_{\text{2}}}\text{O}}=-\frac{{{s}^{{{\text{H}}_{\text{2}}}\text{O}}}{{j}_{\text{c}}}}{nF}-{{S}_{\text{ice}}}\text{+}{{S}_{\text{balance}}}$ | $S_{\text{c}}^{{{\text{H}}_{\text{2}}}\text{O}}=-\nabla \cdot \left( \frac{{{n}_{\text{d}}}}{F}{{\mathrm{i}}_{\text{e}}} \right)-{{S}_{\text{balance}}}$ | |
| Charge | $0=\nabla \cdot \left( \kappa _{\text{eff}}^{\text{e}}\nabla {{\phi }_{\text{e}}} \right)+{{S}_{\phi ,\text{e}}}$ | / | ${{S}_{\phi ,\text{e}}}={{j}_{\text{c}}}$ |
| $0=\nabla \cdot \left( \sigma \nabla {{\phi }_{\text{s}}} \right)+{{S}_{\phi ,\text{s}}}$ | / | ${{S}_{\phi ,\text{s}}}=-{{j}_{\text{c}}}$ | |
在阴极的电化学反应动力学可用Tafel动力学表示:
${{j}_{\text{c}}}=-\frac{1}{V}\int\limits_{A_{\text{inter}}^{\text{r}}}{g\left( {{{{\varepsilon }'}}_{\text{ice}}} \right){{i}_{0,\text{ref}}}\frac{{{C}^{{{\text{O}}_{\text{2}}}}}}{C_{\text{ref}}^{{{\text{O}}_{\text{2}}}}}\exp \left( -\frac{{{\alpha }_{\text{c}}}}{RT}F\eta \right)\text{d}A}\text{ }$ (16)
其中:V为控制体体积;$A_{\text{inter}}^{\text{r}}$为电化学活性反应面积;αc为阴极传输系数;i0,ref为阴极参考交换电流密度;$C_{\text{ref}}^{{{\text{O}}_{\text{2}}}}$为参考氧气浓度;$g\left( {{{{\varepsilon }'}}_{\text{ice}}} \right)$为孔隙中局部含冰量对有效反应面积的影响关系式;η表示过电势,有:
$\eta ={{\phi }_{\text{s}}}-{{\phi }_{\text{e}}}-{{U}_{0}}$ (17)
Table 3 Physical properties and other parameters of the material表3 材料的物性及部分计算参数 |
| Parameters | Value | Units |
|---|---|---|
| Cathode exchange current density i0,ref[2] | 104 | A/m3 |
| Cathode equilibrium potential U0[2] | U0=1.23-9.0×10-4(T-298.15) | V |
| Electro-osmotic drag coefficient nd[21] | ${{n}_{\text{d}}}={{a}_{T}}\left\{ \begin{align} & 0\text{ }\lambda \le \text{3} \\ & {{\left( \frac{\lambda -3}{3} \right)}^{4}}\text{ 3}\lambda \le \text{6} \\ & 1\text{ others} \\ \end{align} \right.$${{a}_{T}}=\left( -\frac{1484.21}{{{T}_{\text{a}}}}+6.062974 \right)\text{ ; if }{{a}_{T}}\text{0}\text{.2, }{{a}_{T}}\text{=0}\text{.2 }$ | H2O/H+ |
| Cathode transfer coefficient αc[2] | 1 | |
| Reference oxygen concentration $C_{\text{ref}}^{{{\text{O}}_{\text{2}}}}$[2] | 40 | mol/m3 |
| Ambient temperature Ta | 243.15(set as default)/248.15/253.15/258.15 | K |
| Current loading | 2 000/3 000(set as default)/4 000 | A/m2 |
| Diffusion coefficient of oxygen/vapor $D_{\text{b}}^{i}$[20] | 1.327 × 10-7/1.28 × 10-5 | m2/s |
| Bruggeman coefficient b[2] | 1.5(set as default)/2.5 | |
| Faraday constant F | 964 87 | C/mol |
| Universal gas constant R | 8.314 | J/(mol·K) |
| Conductivity of carbon black σ[20] | 104 | S/m |
| *Equivalent weight of ionomer EW | 1.0 | kg/mol |
| *Volume fraction of ionomer at mix solid phase ${{\varepsilon }_{\text{m}}}$ | 0.400 7 | |
| Net-water transport coefficient α[22] | -0.1/-0.2(set as default)/-0.4/-0.6 | |
| Density of O2/N2/ice/vapor/carbon/ionomer ρ[2] | 1.4/1.2/0.6/1 800/1 820 | kg/m3 |
| **Specific surface area of phase interface $A_{\text{inter}}^{{}}/\text{V}$ | 2.67 × 1011 | 1/m |
| *Specific electro-chemical active area $A_{\text{inter}}^{\text{r}}/\text{V}$ | 1.064 × 1011 | 1/m |
| *Average pore size $\overline{d}$ | 40 | nm |
| *Porosity of the GDL ${{\varepsilon }_{\text{GDL}}}$ | 0.513 | |
| *Thickness of the GDL ΔX | 200 | μm |
Notes: * data from collaborators; ** results from the reconstructed structure. |
2.3 初边界条件和求解方法
当t = 0 s时,电极内氧气浓度为10 mol/m3,孔相及固相中水浓度依据环境条件（0.1 MPa,启动温度）下相对湿度RH = 38%给定。
氧气由气流通道向CL扩散,水蒸气由CL向气流通道扩散;在PEM/CL边界考虑零流条件,则有:
$\begin{align} & {{\left. D_{\text{eff}}^{{{\text{O}}_{\text{2}}}}\frac{\partial {{C}^{{{\text{O}}_{2}}}}}{\partial x} \right|}_{x=0}}=D_{\text{GDL}}^{{{\text{O}}_{\text{2}}}}\frac{\left( {{\left. {{C}^{{{\text{O}}_{\text{2}}}}} \right|}_{x=0}}-C_{\text{a}}^{{{\text{O}}_{\text{2}}}} \right)}{\Delta X} \\ & {{\left. D_{\text{eff}}^{{{\text{O}}_{\text{2}}}}\frac{\partial {{C}^{{{\text{O}}_{\text{2}}}}}}{\partial x} \right|}_{x={{x}_{\text{L}}}}}=0 \\ \end{align}$ (18)
$\begin{align} & {{\left. D_{\text{g,eff}}^{{{\text{H}}_{\text{2}}}\text{O}}\frac{\partial {{C}^{{{\text{H}}_{\text{2}}}\text{O}}}}{\partial x} \right|}_{x=0}}=D_{\text{GDL}}^{{{\text{H}}_{\text{2}}}\text{O}}\frac{\left( {{\left. {{C}^{{{\text{H}}_{\text{2}}}\text{O}}} \right|}_{x=0}}-C_{\text{a}}^{{{\text{H}}_{\text{2}}}\text{O}} \right)}{\Delta X} \\ & {{\left. D_{\text{g,eff}}^{{{\text{H}}_{\text{2}}}\text{O}}\frac{\partial {{C}^{{{\text{H}}_{\text{2}}}\text{O}}}}{\partial x} \right|}_{x={{x}_{\text{L}}}}}=0 \\ \end{align}$ (19)
其中:$D_{\text{GDL}}^{i}$为扩散层内氧气或水蒸气的扩散系数;ΔX为扩散层厚度。
在PEM/CL的边界上,引入净水传输系数α[20]来衡量催化层与质子交换膜之间水的传输量;在CL/GDL界面,膜中水含量传输流为0,则有:
${{\left. D_{\text{m,eff}}^{{{\text{H}}_{\text{2}}}\text{O}}\frac{\partial C_{\text{m}}^{{{\text{H}}_{\text{2}}}\text{O}}}{\partial x} \right|}_{x=0}}\text{=}0$ ${{\left. D_{\text{m,eff}}^{{{\text{H}}_{\text{2}}}\text{O}}\frac{\partial C_{\text{m}}^{{{\text{H}}_{\text{2}}}\text{O}}}{\partial x} \right|}_{x={{x}_{\text{L}}}}}\text{=}-\alpha \frac{I}{2F}$ (20)
其中,I/2F代表电化学反应产生水的速率。
在x = 0和x = xL处加载电流I,x = xL处的电子电势固定为0 V,故电荷方程的边界条件为:
${{\left. \frac{\partial {{\phi }_{\text{e}}}}{\partial x} \right|}_{x=0}}\text{=}0$ ${{\left. \kappa \frac{\partial {{\phi }_{\text{e}}}}{\partial x} \right|}_{x={{x}_{\text{L}}}}}\text{=}\frac{I}{{{A}_{\text{s}}}}$ (21)
${{\left. \sigma \frac{\partial {{\phi }_{\text{s}}}}{\partial x} \right|}_{x=0}}\text{=}-\frac{I}{{{A}_{\text{s}}}}$ ${{\left. {{\phi }_{\text{s}}} \right|}_{x={{x}_{\text{L}}}}}\text{=}0$ (22)
利用商业计算流体软件Fluent求解,通过用户自定义函数（user-defined functions, UDFs）定义方程的各种源项和物性参数等;基于有限体积法离散求解,瞬态项离散采用完全隐式格式,空间导数项采用中心差分格式差分,采用代数多重网格（algebraic multi-grid, AMG）迭代法加速收敛。模型中数值单元为重建结构（图2）的微元。
3 结果与分析
所求解的算例及其条件如表4所示,电池在恒电流下启动,首先以基础算例分析和讨论模型结果,然后分析讨论变参数（算例1 ~ 算例4）结果。
Table 4 Cases summary表4 算例汇总 |
| Series | α | $\frac{I}{{{A}_{\text{s}}}}$/ (A/m2) | Ta / K | Others |
|---|---|---|---|---|
| Base case | -0.2 | 3 000 | 243.15 | $g({{{\varepsilon }'}_{\text{ice}}})=1-{{{\varepsilon }'}_{\text{ice}}}$; b = 1.5 |
| Case 1 | -0.1/-0.4/-0.6 | 3 000 | 243.15 | $g({{{\varepsilon }'}_{\text{ice}}})=1-{{{\varepsilon }'}_{\text{ice}}}$; b = 1.5 |
| Case 2 | -0.2 | 2 000/4 000 | 243.15 | $g({{{\varepsilon }'}_{\text{ice}}})=1-{{{\varepsilon }'}_{\text{ice}}}$.; b = 1.5 |
| Case 3 | -0.2 | 3 000 | 258.15/253.15/248.15 | $g({{{\varepsilon }'}_{\text{ice}}})=1-{{{\varepsilon }'}_{\text{ice}}}$; b = 1.5 |
| Case 4 | -0.2 | 3 000 | 243.15 | $g({{{\varepsilon }'}_{\text{ice}}})={{\left( 1-{{{{\varepsilon }'}}_{\text{ice}}} \right)}^{2}}$; b = 1.5 $g ({{{\varepsilon }'}_{\text{ice}}})=1-{{{\varepsilon }'}_{\text{ice}}}$; b = 2.5 |
3.1 结冰形貌（基础算例）
模型能预测PEMFC冷启动过程中冰在催化层内的形貌和分布状态。在5 s时开始在孔隙中有冰生成,图6选择绘制了不同时刻孔隙中含冰量的三维分布。9.5 s时冰在催化层内已经初步形成,随着时间的推移,冰的份额增加、逐渐覆盖催化层内孔隙,而且冰占据的区域向PEM侧扩展。
Fig. 6 Ice formation in CL and its time-change图6 催化层内含冰量分布及其随时间演化 |
图7中x截面上的含冰量分布显示出了冰在孔隙内的局部非均匀生成状态。较多冰在孔隙中靠近MS相处生成,孔隙中部位置冰的份额明显较少,这与结晶过程中晶体会优先在固体表面形成的成核理论描述一致。
Fig. 7 Ice distribution on the central x-sectional-plane at 18.5 s图7 t = 18.5 s时催化层中部x截面的含冰量分布 |
Fig. 8 Ice fraction profile along the x-axis and its time-change图8 沿x轴方向的含冰量分布及其随时间演化 |
3.2 结冰对传质和电化学反应的影响（基础算例）
Fig. 9 Distribution of mean oxygen concentration (a) and water content (b) along the x-axis and their evolution over time图9 氧气平均浓度（a）和含水量（b）在沿x轴方向的分布及其随时间演化 |
为比较不同时刻催化层内电化学反应的强度以及含冰量对反应的影响,图10a显示了6.5 s与18.5 s时的输运电流密度（jc）分布。6.5 s与18.5 s时电化学反应均集中在靠近PEM侧（X = 1）的区域,而在阴极CL和阴极GDL的交界面处（X = 0）,jc接近于0;与6.5 s相比,18.5 s时含冰量更高,伴随着冰覆盖活性反应面积,电化学反应区域有向X = 0端转移的趋势。
Fig. 10 Distribution of the transfer current density (a) and protonic current density (b) in CL图10 催化层内输运电流密度（a）与质子电流密度（b）分布 |
当阴极催化层内结冰时,冰不仅会覆盖活性反应面积导致电化学极化损失增加,还会堵塞气体通道增加了浓差极化损失。定义催化层内平均过电势$\left| \eta \right|$为催化层内反应过电势的绝对值以局部输运电流密度（${{{j}'}_{\text{c}}}$）为权数的加权平均值。图11所示的平均过电势随含冰量的演化曲线证实了这一点。
Fig. 11 Evolution of the average overpotential with ice fraction图11 平均过电势随含冰量的演化 |
3.3 电化学有效反应面积与结冰量的关系式
当冰在催化层内形成时,部分活性反应面积将会被覆盖,影响电化学反应的进行。为了研究电化学有效反应面积与含冰量的关系,定义电化学有效反应面积比Seff为催化层内没有被冰覆盖的活性反应面积与总活性面积的比值:
${{S}_{\text{eff}}}=\frac{\int\limits_{A_{\operatorname{int}\text{er}}^{\text{r}}}{g\left( {{{{\varepsilon }'}}_{\text{ice}}} \right)\text{d}A}}{A_{\operatorname{int}\text{er}}^{\text{r}}}$ (23)
Fig. 12 Plots for different conditions of total ice fraction to ECA图12 不同条件下电化学有效反应面积与总含冰量的关系 |
从图12中可知,仿真结果中的Seff均随着含冰量的增加而减少,各Seff曲线的重合度很高（不同算例之间的偏差小于10%）,但与常用的线性关系式有明显偏离;仿真结果的有效反应面积的减少速率更快,覆盖所有的活性反应面积不需要完全堵塞孔隙。根据仿真结果拟合得到:
${{S}_{\text{eff}}}=1-1.249{{s}^{3}}+0.261{{s}^{2}}-0.096s$ (24)
式（24）具有一定的介观尺度基础,预期可改进宏观冷启动模型的计算精度。
4 结 论
基于PEMFC催化层重建结构建立了电池冷启动时多场耦合过程的介孔尺度数值模型。模型考虑了冷启动时电池内的反应物/产物和电荷传输、电化学反应过程,能够用以深入研究电池冷启动时阴极催化层内的电-气-汽（水）传输、电化学反应和结冰等。
数值仿真结果揭示了冷启动时冰的生成对催化层内电化学反应、水的生成和传输、氧气的传输等的影响,再现了冷启动过程中电极内生成冰的形貌和分布,分析并推演得到电化学活性面积与冰含量的宏观关系式。该关系式与常用的线性关系式有较大差异,如果用于宏观模型,预期可以提高模型计算精度。
符号表:
a 水活度
A 面积,m2
B Bruggeman指数
C 摩尔浓度,mol/m3
$\overline{d}$ 平均直径,nm
D 扩散系数,m2/s
EW 离子聚合物的当量质量,kg/mol
F 法拉第常数,96 485 C/mol
I 电流负载,A/m2
i0 交换电流密度,A/m2
ie 质子电流密度,A/m2
j 输运电流密度,A/m3
M 摩尔质量,kg/mol
n 反应电子迁移数
nd 电渗系数
r 速率,s-1
R 通用气体常数,8.314 J/(mol·K)
S 源项,反应面积比
s 整体含冰量
T 温度,K
U0 阴极电化学反应热力学平衡电势,V
α 传输系数,净水传输系数
ε 体积分数
η 过电势,V
κ 离子传导率,S/m
λ 离子聚合物内含水量
ρ 密度,kg/m3
ϕ 电势,V
σ 电子传导率,S/m
角标:
a 环境
c 阴极
e 电解质
eff 有效的
freeze 结冰
g 气体
i 组分
inter 界面
ice 冰
L 总长度
m 模型,离子聚合物
r 反应
ref 参考值
s 固体
sat 饱和值
