文章快速检索    
  同济大学学报(自然科学版)  2019, Vol. 47 Issue (1): 104-112.  DOI: 10.11908/j.issn.0253-374x.2019.01.014
0

引用本文  

许澎, 许思传, 郭鑫, 刘鹏程, 刘星宇. 一种质子交换膜燃料电池冷启动分层数值模型[J]. 同济大学学报(自然科学版), 2019, 47(1): 104-112. DOI: 10.11908/j.issn.0253-374x.2019.01.014.
XU Peng, XU Sichuan, GUO Xin, LIU Pengcheng, LIU Xingyu. A Layered Numerical Model for Cold Start in Proton Exchange Membrane Fuel Cell[J]. Journal of Tongji University (Natural Science), 2019, 47(1): 104-112. DOI: 10.11908/j.issn.0253-374x.2019.01.014

基金项目

国家重点研发计划(2017YFB0102802)

第一作者

许澎(1991—),男,博士生,主要研究方向为燃料电池水热管理、冷启动. E-mail: xpxp1991@126.com

通信作者

许思传(1963—), 男, 教授, 博士生导师, 工学博士, 主要研究方向为燃料电池发动机系统、整车热管理. E-mail: scxu@tongji.edu.cn

文章历史

收稿日期:2018-03-20
一种质子交换膜燃料电池冷启动分层数值模型
许澎 1,2, 许思传 1,2, 郭鑫 1,2, 刘鹏程 1,2, 刘星宇 1,2     
1. 同济大学 汽车学院, 上海 201804;
2. 同济大学 新能源汽车工程中心,上海 201804
摘要:针对质子交换膜燃料电池在低温环境启动时伴随的复杂物理和化学现象,搭建了一种新的燃料电池单电池冷启动数值模型,该模型全面包含了燃料电池低温启动的电化学反应机理、水热传递机理、水相变机理等.研究了催化层聚合物体积分数、孔隙率、质子交换膜厚度三个因素对低温冷启动的影响.最终提出了储水量(WSC)这一预测燃料电池低温冷启动能力的指标.
关键词质子交换膜燃料电池    冷启动    相变    
A Layered Numerical Model for Cold Start in Proton Exchange Membrane Fuel Cell
XU Peng 1,2, XU Sichuan 1,2, GUO Xin 1,2, LIU Pengcheng 1,2, LIU Xingyu 1,2     
1. College of Automotive Studies, Tongji University, Shanghai 201804, China;
2. Clean Energy Automotive Engineering Center, Tongji University, Shanghai 201804, China
Abstract: In terms of complex physical and chemical phenomenon of proton exchange membrane fuel cell (PEMFC) in cold start, a novel numerical model of PEMFC in cold start was developed. The model comprehensively considered electrochemical reaction, thermal and water and phase change mechanism in cold start. Three factors, ionomer volume fraction in catalyst layer, porosity and membrane thickness, were investigated. The indicator of water storage capacity(WSC) was finally introduced to evaluate the cold start capability of PEMFC.
Key words: proton exchange membrane fuel cell (PEMFC)    cold start    phase change    

作为车用动力的质子交换膜燃料电池实现大规模产业化必须要经受启停、高电位、冷启动、电压循环、大电流、空气杂质、反极等复杂环境和严苛工况的考验[1].其中冷启动工况指燃料电池在0 ℃以下的低温环境中启动至常温(约60~80 ℃)运行的过程.美国能源部针对净功率80 kW燃料电池动力系统在冷启动方面提出了2020年的技术目标:动力系统在低温环境静置8 h后可以满足-30 ℃自启动和-40 ℃辅助启动,从-20 ℃的低温环境中运行至50 %额定功率用时30 s以内,并且耗能小于5 MJ[2].燃料电池低温冷启动技术的难点在于低温启动时冰的形成导致电化学反应活性位点的覆盖和反应气传输通道的阻碍,短期造成电池启动缓慢甚至失败.长期来看,由于水结冰的循环相变造成电池内部多孔介质结构不可逆的破坏,比如:质子交换膜和催化层间隙变大,催化层微孔的致密化和崩塌,催化剂颗粒的粗化等,影响燃料电池的耐久性和可靠性[3-4].提高燃料电池低温冷启动能力,首先需要从材料、结构方面来强化低温环境下电池的功率输出能力、排水能力、传质能力,比如丰田公司旗下Mirai燃料电池汽车搭载了全新3D精细网眼流场板的燃料电池来优化排水能力[5],质子交换膜轻薄化来提高功率输出能力等;再者,强化利用电化学反应产生的“废热”[6]或者外部热源来快速提高电池温度.

由于燃料电池低温冷启动试验的高成本、低可控性等特点,数值模拟成为低温冷启动研究的一种重要手段.目前有关燃料电池低温冷启动的数值模型从维度上主要分为一维[7-13]和三维[14-17].三维冷启动模型尤以宾夕法尼亚州州立大学王朝阳[14]的均相模型和滑铁卢大学李献国、天津大学焦魁[15]的分相模型为代表.由于综合考虑了电池厚度方向和平面方向的复杂物理、化学现象,三维冷启动模型在数值高精确度、后处理简易、直观性等方面更胜一筹.但是,多维度计算的网格数量较大,对时间成本和计算设备要求较高.而一维冷启动数值模拟由于计算效率高和拓展性高等特点,且可兼顾燃料电池冷启动伴随的一系列物理化学现象,受到关注.

燃料电池一维冷启动数值模型目前主要有:Sundaresan[8]、Khandelwal[9]以及周怡博[10-13]模型.李义[7]等建立了燃料电池系统一维冷启动模型,得出大电流密度启动模式有利于系统低温冷启动.该集总参数模型将电堆视为整体,忽略燃料电池内部的电化学反应过程、水热传递过程、相变过程和传质过程.Sundaresan[5]建立了一个单体燃料电池冷启动一维热模型,该模型主要考虑冷启动过程中电池一维方向上的温度分布特征,讨论了双极板的热容量和加热方法对冷启动的影响.Khandelwal[5]在前人工作基础上,通过一维热模型,从冷启动时间和耗能两个指标着重分析了低温冷启动能力,并且指出低温冷启动存在一个最佳的启动电流,同时建议采用端板绝缘的方法来加速电堆的温升过程.以上Sundaresan和Khandelwal两位研究学者的模型均聚焦热分析,忽略电池内部一系列结冰、冷凝等复杂的相变现象,因而不能全面反映冷启动过程中电池内部的水热传递机理.周怡博[5]采用有限体积法的思想通过求解一组偏微分方程组对燃料电池电堆进行冷启动数值模拟,研究可变加热和载荷功率的方法对冷启动的影响.

基于分层集总参数的思想,本文建立了一种新的质子交换膜燃料电池冷启动分层数值模型,该模型将单电池分为质子交换膜、阴阳极催化层、阴阳极气体扩散层、阴阳极流道、阴阳极冷却流道、阴阳极端板11层,将每一层假设为一个控制体,在控制体的中心集中了每一层的参数(温度、水蒸气浓度、冰体积分数、液态水体积分数、膜结合水含量和膜冻结水含量),全面考虑低温环境下电池内部水的液态水、冰、水蒸气、膜结合水和膜冻结水不同形态及其相变过程,研究燃料电池低温环境下水热传递机理,最后,提出了预测低温冷启动能力的WSC指标.

1 PEMFC冷启动分层数值模型

图 1所示,燃料电池的基本结构包含阴阳极端板、阴阳极冷却流道、阴阳极流场、阴阳极气体扩散层、阴阳极催化层和质子交换膜.由于具备多孔介质孔隙结构,水以水蒸气、液态水、冰三种形式存在于气体扩散层和催化层的孔隙中.催化层具有传导质子的聚合物,该聚合物和质子交换膜中聚四氟乙烯基质性质相同,所以催化层的聚合物中也存在膜结合水.在燃料电池常温运行下,水在质子交换膜中主要以水合氢离子的形式存在,低温情况下,水分子和氢离子的结合力较小,一部分水分子析出形成膜冻结水.电池启动条件下,由于气体流道中气体的吹扫作用,流道中冰的含量忽略不计.表 1列出了燃料电池内部不同组件中水的存在形式.本文建立的一维冷启动分层数值模型对温度和不同形态的水分别建立相应方程,计算域全局进行温度场求解,水蒸气计算域分布在阴阳极催化层、阴阳极气体扩散层和流道,液态水计算域分布在阴阳极催化层和气体扩散层,膜结合水计算域分布在阴阳极催化层和质子交换膜,膜冻结水计算域分布在质子交换膜内.

图 1 燃料电池冷启动过程水热传递 Fig.1 Thermal and water transport mechanism inside PEMFC during cold start
下载CSV 表 1 燃料电池内部低温环境下水的存在形式 Tab.1 Status of water for PEMFC in cold weather

模型有以下假设:①所有气体为理想气体;②电化学反应产生的水初始时刻以膜结合水形式存在;③忽略电池厚度方向的传热传质,扩散是主要传递形式;④电池内部压力变化忽略;⑤不考虑重力影响;⑥阴阳极催化层聚合物中的膜结合水饱和析出的冰全部存在孔隙中;⑦每一层的物理参数集中在各层中心位置.以下将对质子交换膜、阴阳极催化层、阴阳极气体扩散层、阴阳极流场、阴阳极冷却流道和阴阳极端板分别做详细的介绍.

1.1 质子交换膜

目前应用最广泛的全氟磺酸型质子交换膜由碳氟主链和带有磺酸基团的醚支链构成.磺酸基团紧紧依附在膜材料上很难移动,氢离子和磺酸集团之间由于吸引力形成亲水性离子团,氢离子的传导有两种形式[18]:输运机制和跳跃机制.其中,输运机制主要表现为水合氢离子在浓度梯度下的扩散行为,而当质子交换膜中的含水量大到使得聚合物链之间距离变短甚至相连,氢离子的传导即为跳跃机制.质子交换膜必须具备良好的湿润性以确保高质子传导率.一般,聚合物的湿润性用水含量λn来表示,即每个磺酸集团携带的水分子数,而质子交换膜聚合物中水的浓度和水含量是相关的,所以水含量可以表示如下:

$ {\lambda _n} = \frac{E}{{{\rho _{\rm{m}}}}}{c_{{\rm{h2o}}}} $ (1)

式中:ρm为膜的密度;E为膜的等效质量;ch2o为膜内水的摩尔浓度.需要注意的是,低温冷启动过程中膜含水量指的是膜结合水含量.

首先建立质子交换膜的能量方程如下:

$ {\rho _{\rm{m}}}{C_{\rm{m}}}\frac{{{\rm{d}}{T_{\rm{m}}}}}{{{\rm{d}}t}} = {Q_{{\rm{o}} - {\rm{m}}}} + {Q_{{\rm{ac}} - {\rm{m}}}} + {Q_{{\rm{cc}} - {\rm{m}}}} + {S_{{\rm{t}} - {\rm{m}}}} $ (2)

式中:Cm为比热容;Tm为温度;Qo-m为膜阻抗产生的焦耳热通量;Qac-m为阳极催化层向膜扩散的热通量;Qcc-m为阴极催化层向膜扩散的热通量;St-m为源项.以上各项表达式分别表示如下:

$ {Q_{{\rm{o}} - {\rm{m}}}} = \frac{{{I^2}{A_{\rm{m}}}}}{{{\delta _{\rm{m}}}}} $ (3)

式中:I为电流密度;δm为膜厚;Am为面积电阻,与质子交换膜的厚度和电导率有关,表示如下:

$ {A_{\rm{m}}} = \frac{{{\delta _{\rm{m}}}}}{{{\kappa _{\rm{m}}}}} $ (4)
$ {\kappa _{\rm{m}}} = \left( {0.513\;9\lambda - 0.326} \right)\exp \left[ {1\;268\left( {\frac{1}{{303.15}} - \frac{1}{{{T_{\rm{m}}}}}} \right)} \right] $ (5)

式中:κm为质子电导率.

$ {Q_{{\rm{cc}} - {\rm{m}}}} = \frac{{{T_{{\rm{cc}}}} - {T_{\rm{m}}}}}{{{\delta _{{\rm{cc}}}}/2 + {\delta _{\rm{m}}}/2}}\frac{{{k_{{\rm{cc}} - {\rm{m}}}}}}{{{\delta _{\rm{m}}}}} $ (6)
$ {Q_{{\rm{ac}} - {\rm{m}}}} = \frac{{{T_{{\rm{ac}}}} - {T_{\rm{m}}}}}{{{\delta _{{\rm{ac}}}}/2 + {\delta _{\rm{m}}}/2}}\frac{{{k_{{\rm{ac}} - {\rm{m}}}}}}{{{\delta _{\rm{m}}}}} $ (7)

式中:TccTacTm分别为阴阳极催化层和膜的温度;kcc-mkac-m分别为阴阳极催化层到质子交换膜的有效导热系数,表示如下:

$ {k_{{\rm{cc}} - {\rm{m}}}} = \frac{{{\delta _{\rm{m}}}/2 + {\delta _{{\rm{cc}}}}/2}}{{{\delta _{\rm{m}}}/2/{k_{\rm{m}}} + {\delta _{{\rm{cc}}}}/2/{k_{{\rm{cc}}}}}} $ (8)
$ {k_{{\rm{ac}} - {\rm{m}}}} = \frac{{{\delta _{\rm{m}}}/2 + {\delta _{{\rm{ac}}}}/2}}{{{\delta _{\rm{m}}}/2/{k_{\rm{m}}} + {\delta _{{\rm{ac}}}}/2/{k_{{\rm{ac}}}}0}} $ (9)

式中:δccδac分别表示阴阳极催化层的厚度;kmkcckcc分别表示膜、阴阳极催化层的导热系数;

$ {S_{{\rm{t}} - {\rm{m}}}} = {h_{{\rm{fusn}}}}\left( {{M_{{\rm{h2o}}}}{S_{{\rm{nf}} - {\rm{m}}}}} \right) $ (10)

式中:hfusn为相变潜热值;Mh2o为水的摩尔质量;Snf-m为膜结合水的源项.

其次建立膜结合水方程为

$ \frac{{{\rho _{\rm{m}}}}}{{EW}}w\frac{{{\rm{d}}{\lambda _{\rm{n}}}}}{{{\rm{d}}t}} = {\lambda _{{\rm{m}} - {\rm{ac}}}} + {\lambda _{{\rm{m}} - {\rm{cc}}}} - {S_{{\rm{nf}} - {\rm{m}}}} $ (11)

式中:w为催化层中聚合物的体积分数;λn为膜内结合水含量; λm-ccλm-ac分别为膜结合水向阴阳极催化层扩散的扩散通量.

$ {\lambda _{{\rm{m}} - {\rm{ac}}}} = \frac{{{\rho _{\rm{m}}}}}{{{\delta _{\rm{m}}}EW}}{D_{{\rm{n}} - {\rm{m}} - {\rm{ac}}}}\frac{{{\lambda _{{\rm{n}} - {\rm{ac}}}} - {\lambda _{\rm{n}}}}}{{{\delta _{\rm{m}}}/2 + {\delta _{{\rm{ac}}}}/2}} $ (12)
$ {\lambda _{{\rm{m}} - {\rm{cc}}}} = \frac{{{\rho _{\rm{m}}}}}{{{\delta _{\rm{m}}}EW}}{D_{{\rm{n}} - {\rm{m}} - {\rm{cc}}}}\frac{{{\lambda _{{\rm{n}} - {\rm{cc}}}} - {\lambda _{\rm{n}}}}}{{{\delta _{\rm{m}}}/2 + {\delta _{{\rm{cc}}}}/2}} $ (13)

式中:λn-ccλn-ac分别为阴阳极催化层的膜结合水含量;Dn-m-ccDn-m-ac分别为膜结合水向阴阳极催化层扩散的有效扩散系数,表示如下:

$ {D_{n - {\rm{m}} - {\rm{cc}}}} = \frac{{{\delta _{\rm{m}}}/2 + {\delta _{{\rm{cc}}}}/2}}{{{\delta _{\rm{m}}}/2/{D_{{\rm{n}} - {\rm{m}}}} + {\delta _{{\rm{cc}}}}/2/{D_{n - {\rm{cc}}}}}} $ (14)
$ {D_{n - {\rm{m}} - {\rm{ac}}}} = \frac{{{\delta _{\rm{m}}}/2 + {\delta _{{\rm{ac}}}}/2}}{{{\delta _{\rm{m}}}/2/{D_{{\rm{n}} - {\rm{m}}}} + {\delta _{{\rm{ac}}}}/2/{D_{n - {\rm{ac}}}}}} $ (15)

式中:Dn-mDn-ccDn-ac为膜的结合水扩散通量,表示如下:

$ {D_n} = \left\{ \begin{array}{l} 2.692\;66 \times {10^{ - 10}}\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;{\lambda _n} \le 2\\ {10^{ - 10}}\exp \left[ {2\;416\left( {\frac{1}{{303}} - \frac{1}{T}} \right)} \right]\left[ {0.87\left( {3 - {\lambda _n}} \right) + 2.95\left( {{\lambda _n} - 2} \right)} \right]\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;2 < {\lambda _n} \le 3\\ {10^{ - 10}}\exp \left[ {2\;416\left( {\frac{1}{{303}} - \frac{1}{T}} \right)} \right]\left[ {2.95\left( {4 - {\lambda _n}} \right) + 1.642\;454\left( {{\lambda _n} - 3} \right)} \right]\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;3 < {\lambda _n} \le 4\\ {10^{ - 10}}\exp \left[ {2\;416\left( {\frac{1}{{303}} - \frac{1}{T}} \right)} \right]\left[ {2.563 - 0.33{\lambda _n} + 0.026\;4{\lambda _n} - 0.000\;671\lambda _n^3} \right.\;\;\;\;{\lambda _n} > 4 \end{array} \right. $ (16)

最后膜的冻结水方程为

$ {\rho _{\rm{m}}}w\frac{{{\rm{d}}{\lambda _f}}}{{{\rm{d}}t}} = {S_{nf - {\rm{m}}}} $ (17)

式中:λf为膜的冻结水含量.

1.2 催化层 1.2.1 阴极催化层

燃料电池催化层兼具质子交换膜和多孔介质的材料属性,是构成电化学反应三相界面的区域.如图 2所示,催化层主要由传导质子的聚合物、传导电子的碳颗粒或碳纤维、催化剂颗粒组成.和质子交换膜中聚合物属性相似,在低温启动初始阶段,催化层聚合物中水以膜结合水形式存在,到达饱和状态以后从聚合物中析出,并且在催化层的多孔介质中生成冰.

图 2 催化层结构 Fig.2 Representation of structure in catalyst layer

首先建立阴极催化层能量方程如下:

$ {\overline {\rho C} _{{\rm{cc}}}}\frac{{{\rm{d}}{T_{{\rm{cc}}}}}}{{{\rm{d}}t}} = {Q_{{\rm{cg}} - {\rm{cc}}}} + {Q_{{\rm{m}} - {\rm{cc}}}} - {Q_{{\rm{ent}}}} + {Q_{{\rm{o}} - {\rm{c}}}} + {Q_{{\rm{ac}} - {\rm{c}}}} + {S_{{\rm{pc}} - {\rm{c}}}} $ (18)

式中:$ {\overline {\rho C} _{{\rm{cc}}}} $为催化层的热容;Tcc为阴极催化层的温度;Qcg-cc为从阴极扩散层到阴极催化层的热扩散通量;Qm-cc为从质子交换膜向阴极催化层的热扩散通量;Qent为可逆热通量[10]Qo-c为欧姆热通量;Qac-c为活化热通量;Spc-c为源项,以上各项表达式如下:

$ \begin{array}{l} {\overline {\rho C} _{{\rm{cc}}}} = \left( {\left( {{s_{{\rm{ice}} - {\rm{cc}}}}{\rho _{{\rm{ice}}}}{C_{{\rm{ice}}}} + {s_{{\rm{lq}} - {\rm{cc}}}}{\rho _{{\rm{lq}}}}{C_{{\rm{lq}}}}} \right){\varepsilon _{{\rm{cc}}}} + } \right.\\ \;\;\;\;\;\;\;\;\;\;\;\left. {{\rho _{{\rm{cc}}}}{C_{{\rm{cc}}}}\left( {1 - {\varepsilon _{{\rm{cc}}}} - {w_{\rm{c}}}} \right) + {w_{\rm{c}}}{\rho _{\rm{m}}}{C_{\rm{m}}}} \right) \end{array} $ (19)

式中:slq-cc为催化层多孔介质中液态水饱和度;sice-cc为催化层多孔介质中冰的体积分数; εcc为催化层孔隙率;CiceClq分别表示冰和液态水的比热容;ρiceρlq分别表示冰和水的密度;wc为阴极催化层聚合物体积分数.

$ {Q_{{\rm{m}} - {\rm{cc}}}} = \frac{{{T_{\rm{m}}} - {T_{{\rm{cc}}}}}}{{{\delta _{\rm{m}}}/2 + {\delta _{{\rm{cc}}}}/2}}\frac{{{k_{{\rm{m}} - {\rm{cc}}}}}}{{{\delta _{{\rm{cc}}}}}} $ (20)
$ {Q_{{\rm{cg}} - {\rm{cc}}}} = \frac{{{T_{{\rm{cg}}}} - {T_{{\rm{cc}}}}}}{{{\delta _{{\rm{cg}}}}/2 + {\delta _{{\rm{cc}}}}/2}}\frac{{{k_{{\rm{cg}} - {\rm{cc}}}}}}{{{\delta _{{\rm{cc}}}}}} $ (21)
$ {Q_{{\rm{ent}}}} = \frac{{{T_{{\rm{cc}}}}I\Delta S}}{{2F{\delta _{{\rm{cc}}}}}} $ (22)

式中:δccδcg为阴极催化层和扩散层的厚度;ΔS为可逆热;I为工作电流;F为法拉第常数;km-cckcg-cc分别为质子交换膜、阴极扩散层到阴极催化层的有效导热系数.

$ {Q_{{\rm{o}} - {\rm{cc}}}} = \frac{{{I^2}{A_{{\rm{cc}}}}}}{{3{\delta _{{\rm{cc}}}}}} $ (23)
$ {Q_{{\rm{ac}} - {\rm{c}}}} = \frac{{I\eta }}{{{\delta _{{\rm{cc}}}}}} $ (24)

式中:η为活化损失; Acc为阴极催化层面积电阻.

$ \begin{array}{l} {S_{{\rm{pc}} - {\rm{cc}}}} = \left( {{h_{{\rm{fusn}}}}\left( {{M_{{\rm{h2o}}}}{S_{{\rm{ni}} - {\rm{cc}}}} + {S_{{\rm{li}} - {\rm{cc}}}}} \right) + \left( {{h_{{\rm{fusn}}}} + {h_{{\rm{cond}}}}} \right){S_{{\rm{vi}} - {\rm{cc}}}} + } \right.\\ \;\;\;\;\;\;\;\;\;\left. {{h_{{\rm{cond}}}}\left( {{S_{{\rm{vi}} - {\rm{cc}}}} - {S_{{\rm{nv}} - {\rm{cc}}}}{M_{{\rm{h2o}}}}} \right)} \right) \end{array} $ (25)

式中:hcond为冷凝潜热值;Sni-cc为膜结合水向冰转化的源项;Sli-cc为液态水向冰转化的源项;Svi-cc为水蒸气向冰转化的源项;Snv-cc为膜结合水向水蒸气转化的源项.

然后建立膜结合水方程,如下:

$ \frac{{{\rho _{\rm{m}}}}}{{EW}}w\frac{{{\rm{d}}{\lambda _{{\rm{n}} - {\rm{cc}}}}}}{{{\rm{d}}t}} = \frac{i}{{2F{\delta _{{\rm{cc}}}}}} + {\lambda _{{\rm{n}} - {\rm{m}} - {\rm{cc}}}} - {S_{{\rm{eod}}}} - {S_{{\rm{nv}} - {\rm{cc}}}} - {S_{{\rm{ni}} - {\rm{cc}}}} $ (26)

式中:λn-cc为阴极催化层膜结合水含量;右端第一项为阴极催化层电化学反应生成的水含量;λn-m-cc为质子交换膜向阴极催化层的膜结合水扩散通量;Seod为电渗通量;Snv-ccSni-cc为源项,以上各项分别表示如下:

$ {\lambda _{{\rm{n}} - {\rm{m}} - {\rm{cc}}}} = \frac{{{D_{{\rm{n}} - {\rm{m}} - {\rm{cc}}}}}}{{{\delta _{{\rm{cc}}}}}}\frac{{{\lambda _{\rm{n}}} - {\lambda _{{\rm{n}} - {\rm{cc}}}}}}{{{\delta _{\rm{m}}}/2 + {\delta _{{\rm{cc}}}}/2}}\frac{{{\rho _{\rm{m}}}}}{{EW}} $ (27)
$ {S_{{\rm{eod}}}} = {n_{\rm{d}}}\frac{I}{{{\delta _{{\rm{cc}}}}F}}{\lambda _{{\rm{n}} - {\rm{ac}}}} $ (28)

式中:Dn-m-cc为质子交换膜到阴极催化层膜结合水的扩散通量;nd为电渗拖拽系数;λn-ac为阳极催化层的膜结合水含量.

其次建立催化层多孔介质冰方程,如下:

$ {\rho _{{\rm{ice}}}}{\varepsilon _{{\rm{cc}}}}\frac{{{\rm{d}}{s_{{\rm{ice}} - {\rm{cc}}}}}}{{{\rm{d}}t}} = {S_{{\rm{ice}} - {\rm{cc}}}} $ (29)
$ {S_{{\rm{ice}} - {\rm{cc}}}} = {S_{{\rm{vi}} - {\rm{cc}}}} + {S_{{\rm{li}} - {\rm{cc}}}} + {S_{{\rm{ni}} - {\rm{cc}}}}{M_{{\rm{h2o}}}} $ (30)

式中:sice-cc为阴极催化层冰体积分数.

最后建立催化层多孔介质水蒸气方程,如下:

$ {E_{{\rm{cc}}}}\frac{{{\rm{d}}{v_{{\rm{cc}}}}}}{{{\rm{d}}t}} = {S_{{\rm{nv}} - {\rm{cc}}}} - {S_{{\rm{vl}} - {\rm{cc}}}} - {S_{{\rm{vi}} - {\rm{cc}}}} - {v_{{\rm{d}} - {\rm{cc}}}} $ (31)

式中:vcc为阴极催化层水蒸气摩尔浓度;Ecc为多孔介质有效孔隙率;Svl-cc为水蒸气向液态水转化的源项;vd-cc-cg为阴极催化层到阴极扩散层水蒸气的扩散通量.表示如下:

$ {E_{{\rm{cc}}}} = \left( {1 - {s_{{\rm{lq}} - {\rm{cc}}}} - {s_{{\rm{ice}} - {\rm{cc}}}}} \right){\varepsilon _{{\rm{cc}}}} $ (32)
$ {v_{{\rm{d}} - {\rm{cc}} - {\rm{cg}}}} = \frac{{{D_{{\rm{v}} - {\rm{cc}} - {\rm{cg}}}}}}{{{\delta _{{\rm{cc}}}}}}\frac{{{v_{{\rm{cc}}}} - {v_{{\rm{cg}}}}}}{{{\delta _{{\rm{cc}}}}/2 + {\delta _{{\rm{cg}}}}/2}} $ (33)
1.2.2 阳极催化层

首先建立阳极催化层能量方程,如下:

$ {\overline {\rho C} _{{\rm{ac}}}}\frac{{{\rm{d}}{T_{{\rm{ac}}}}}}{{{\rm{d}}t}} = {Q_{{\rm{ag}} - {\rm{ac}}}} + {Q_{{\rm{m}} - {\rm{ac}}}} + {Q_{{\rm{o}} - {\rm{ac}}}} + {S_{{\rm{pc}} - {\rm{ac}}}} $ (34)

式中:$ {\overline {\rho C} _{{\rm{ac}}}} $为阳极催化层的热容,表达式可参考式(19);Qag-ac为阳极扩散层到催化层的扩散热通量,参考式(21);Qm-ac为质子交换膜到催化层的扩散热通量,参考式(20);Qo-ac为阳极催化层的欧姆热通量,参考式(23);Spc-ac为源项,其表达式可参考式(25).以上表达式分别表示如下:

然后建立膜结合水方程,如下:

$ \frac{{{\rho _{\rm{m}}}}}{{EW}}w\frac{{{\rm{d}}{\lambda _{{\rm{n}} - {\rm{ac}}}}}}{{{\rm{d}}t}} = {\lambda _{{\rm{n}} - {\rm{m}} - {\rm{ac}}}} + {S_{{\rm{eod}}}} - {S_{{\rm{nv}} - {\rm{ac}}}} - {S_{{\rm{ni}} - {\rm{ac}}}} $ (35)

式中:λn-m-ac为质子交换膜向阳极催化层的膜结合水扩散通量,参考式(27);Snv-acSni-ac为源项.

其次建立催化层多孔介质冰方程,如下:

$ {\rho _{{\rm{ice}}}}{\delta _{{\rm{ac}}}}\frac{{{\rm{d}}{s_{{\rm{ice}} - {\rm{ac}}}}}}{{{\rm{d}}t}} = {S_{{\rm{ice}} - {\rm{ac}}}} $ (36)

式中:sice-ac为阳极催化层的冰体积分数;Sice-ac为源项.

$ {S_{{\rm{ice}} - {\rm{cc}}}} = {S_{{\rm{vi}} - {\rm{ac}}}} + {S_{{\rm{li}} - {\rm{ac}}}} + {S_{{\rm{ni}} - {\rm{ac}}}}{M_{{\rm{h2o}}}} $ (37)

最后建立催化层多孔介质水蒸气方程,如下:

$ {E_{{\rm{ac}}}}\frac{{{\rm{d}}{v_{{\rm{ac}}}}}}{{{\rm{d}}t}} = {S_{{\rm{nv}} - {\rm{ac}}}} - {S_{{\rm{vl}} - {\rm{ac}}}} - {S_{{\rm{vi}} - {\rm{ac}}}} - {v_{{\rm{d}} - {\rm{ac}} - {\rm{ag}}}} $ (38)

式中:vac为阳极催化层水蒸气摩尔浓度;Eac为多孔介质有效孔隙率,参考式(32);vd-ac-cg为阳极催化层到阳极极扩散层水蒸气的扩散通量,表达式可参考式(33).

1.3 气体扩散层模型

多孔介质中水以气态、液态、冰三种形式存在.根据水的三相图可知,温度在高于零度时,水的相变发生在液态和水蒸气之间,而当温度在零度以下,由于水结成冰,本文假设多孔介质中不存在过冷水.饱和水蒸气压力和当地水蒸气压力的差值决定着相变发生在液态和气态之间还是气态和冰之间.当温度大于冰点时,如果水蒸气压力大于饱和压力,水蒸气液化,反之蒸发;当温度低于冰点,如果水蒸气压力高于饱和压力,水蒸气凝华,反之不会升华,因为燃料电池的工作压力在1~2个大气压强左右.由于阴阳极多孔介质均承担传热、传质功能且传递机理相同,以下公式适用于阴极和阳极气体扩散层,没有对阴阳极分开阐述.首先建立扩散层能量方程,如下:

$ {\overline {\rho C} _{\rm{g}}}\frac{{{\rm{d}}{T_{\rm{g}}}}}{{{\rm{d}}t}} = {Q_{{\rm{c}} - {\rm{g}}}} + {Q_{{\rm{ch}} - {\rm{g}}}} + {Q_{{\rm{o}} - {\rm{g}}}} + {S_{{\rm{pc}} - {\rm{g}}}} $ (39)

式中:Tg为扩散层的温度;Qc-g为从催化层到扩散层的热扩散通量,表达式可参考式(21);Qch-g为从流场到扩散层的热扩散通量;Qo-g为扩散层的欧姆热通量, 表达式可参考式(23);Spc-g为源项.

$ {\overline {\rho C} _{\rm{g}}} = \left( {\left( {{s_{{\rm{ice}} - {\rm{g}}}}{\rho _{{\rm{ice}}}}{C_{{\rm{ice}}}} + {s_{{\rm{lq}} - {\rm{g}}}}{\rho _{{\rm{lq}}}}{C_{{\rm{lq}}}}} \right){\varepsilon _{\rm{c}}} + {\rho _{\rm{g}}}{C_{\rm{g}}}\left( {1 - {\varepsilon _{\rm{g}}}} \right)} \right) $ (40)
$ {S_{{\rm{pc}} - {\rm{g}}}} = \left( {{h_{{\rm{fusn}}}}{S_{{\rm{li}} - {\rm{g}}}} + \left( {{h_{{\rm{fusn}}}} + {h_{{\rm{cond}}}}} \right){S_{{\rm{vi}} - {\rm{g}}}} + {h_{{\rm{cond}}}}{S_{{\rm{vl}} - {\rm{g}}}}} \right) $ (41)
$ {Q_{{\rm{ch}} - {\rm{g}}}} = \frac{{{T_{{\rm{ch}}}} - {T_{\rm{g}}}}}{{{\delta _{{\rm{ch}}}}/2 + {\delta _{\rm{g}}}/2}}\frac{{{k_{{\rm{ch}} - {\rm{g}}}}}}{{{\delta _{\rm{g}}}}} $ (42)

然后建立冰方程,如下:

$ {\rho _{{\rm{ice}}}}{\varepsilon _{\rm{g}}}\frac{{{\rm{d}}{s_{{\rm{ice}} - {\rm{g}}}}}}{{{\rm{d}}t}} = {S_{{\rm{vi}} - {\rm{g}}}} + {S_{{\rm{li}} - {\rm{g}}}} $ (43)

式中:εg为扩散层的孔隙率;sice-g为扩散层冰体积分数;Svi-gSli-g分别为扩散层中水蒸气到冰的源项和液态水到冰的源项.

最后建立水蒸气方程,如下:

$ {E_{\rm{g}}}\frac{{{\rm{d}}{v_{\rm{g}}}}}{{{\rm{d}}t}} = {v_{{\rm{c}} - {\rm{g}}}} - {v_{{\rm{g}} - {\rm{ch}}}} - {S_{{\rm{vi}} - {\rm{g}}}} - {S_{{\rm{vl}} - {\rm{g}}}} $ (44)

式中:Eg为扩散层有效孔隙率;vg为水蒸气摩尔浓度;vc-g为催化层到扩散层水蒸气扩散通量,表达式可参考式(33);vg-ch为扩散层到流道水蒸气扩散通量.

$ {v_{{\rm{c}} - {\rm{g}}}} = \frac{{{D_{{\rm{v}} - {\rm{c}} - {\rm{g}}}}}}{{{\delta _{\rm{g}}}}}\frac{{{v_{\rm{c}}} - {v_{\rm{g}}}}}{{{\delta _{\rm{g}}}/2 + {\delta _{\rm{c}}}/2}} $ (45)
$ {v_{{\rm{g}} - {\rm{ch}}}} = \frac{{{D_{{\rm{v}} - {\rm{g}} - {\rm{ch}}}}}}{{{\delta _{\rm{g}}}}}\frac{{{v_{\rm{g}}} - {v_{{\rm{ch}}}}}}{{{\delta _{\rm{g}}}/2 + {\delta _{{\rm{ch}}}}/2}} $ (46)

式中:Dv-c-gDv-g-ch分别为水蒸气从催化层到扩散层、扩散层到流场的扩散系数,表达式可参考式(14).

1.4 流场模型

由于阴阳极流场公式相同,以下公式适用于阴极和阳极流场,没有对阴阳极分开阐述.首先建立流场能量方程,如下:

$ {\rho _{{\rm{ch}}}}{C_{{\rm{ch}}}}\frac{{{\rm{d}}{T_{{\rm{ch}}}}}}{{{\rm{d}}t}} = {Q_{{\rm{g}} - {\rm{ch}}}} + {Q_{{\rm{cp}} - {\rm{ch}}}} + {Q_{{\rm{o}} - {\rm{ch}}}} $ (47)

式中:ρch为流场的密度;Cch为流场的比热容;Qg-ch为气体扩散层到流场的热扩散通量,表达式可参考式(20);Qcp-ch为双极板到流场的热扩散通量, 表达式可参考式(21);Qo-ch为流场的欧姆热通量,表达式可参考式(23).

然后建立流场水蒸气方程,如下:

$ \frac{{{\rm{d}}{v_{{\rm{ch}}}}}}{{{\rm{d}}t}} = {v_{{\rm{g}} - {\rm{ch}}}} - {v_{{\rm{ch}} - {\rm{out}}}} $ (48)

式中:vch为流场中水蒸气摩尔浓度;vg-ch为扩散层向流场水蒸气扩散通量,表达式可参考式(46).

$ {v_{{\rm{ch}} - {\rm{out}}}} = {v_{{\rm{ch}}}}{q_{{\rm{ch}}}} $ (49)

式中:qch为流道中气体流量[3].

1.5 冷却流道模型

该部分适用于阴阳极流道.建立冷却流道能量方程,如下:

$ {\rho _{{\rm{cp}}}}{C_{{\rm{cp}}}}\frac{{{\rm{d}}{T_{{\rm{cp}}}}}}{{{\rm{d}}t}} = {Q_{{\rm{ep}} - {\rm{cp}}}} - {Q_{{\rm{cp}} - {\rm{cc}}}} $ (50)

式中:ρcp为冷却流道的密度;Ccp为冷却流道的比热容;Tcp为冷却流道的温度;Qep-cp为端板向冷却流道的热扩散通量,表达式可参考式(21);Qcp-cc为冷却流道向冷却液的热扩散通量.

$ {Q_{{\rm{cp}} - {\rm{cc}}}} = \left( {{T_{{\rm{cp}}}} - {T_{{\rm{cc}}}}} \right)\frac{{{h_{\rm{c}}}}}{{{\delta _{{\rm{cp}}}}}} $ (51)

式中:Tcc为冷却液的温度;hc为冷却液与冷却流道的换热系数.

1.6 端板模型

该部分适用于阴阳极端板.建立能量方程,如下:

$ {\rho _{{\rm{ep}}}}{C_{{\rm{ep}}}}\frac{{{\rm{d}}{T_{{\rm{ep}}}}}}{{{\rm{d}}t}} = - {Q_{{\rm{ep}} - {\rm{cp}}}} + {Q_{{\rm{bd}} - {\rm{ep}}}} $ (52)

式中:ρep为端板的密度;Cep为端板的比热容;Tep为端板的温度;Qbd-ep为外界和端板的换热量.

$ {Q_{{\rm{bd}} - {\rm{ep}}}} = \left( {{T_{{\rm{bd}}}} - {T_{{\rm{ep}}}}} \right)\frac{{{h_{{\rm{surr}}}}}}{{{\delta _{{\rm{ep}}}}}} $ (53)

式中:Tbd为环境温度;hsurr为换热系数;δep为端板厚度.

1.7 相变源项

对于相变源项一一说明.当多孔介质中局部温度大于冰点,且水蒸气浓度大于饱和水蒸气浓度时,水蒸气向液态水转变,反之液态水向水蒸气转变,当局部温度小于冰点时,假设孔隙中不存在固体冰,不发生水蒸气和液态水到冰的相变,所以水蒸气和液态水之间转化的源项表达式为

$ {S_{{\rm{v}} - {\rm{l}}}} = \left\{ \begin{array}{l} \left\{ \begin{array}{l} {\gamma _{{\rm{cond}}}}\varepsilon \left( {1 - {s_{{\rm{lq}}}} - {s_{{\rm{ice}}}}} \right)\left( {{c_{{\rm{vp}}}} - {c_{{\rm{sat}}}}} \right)\;\;\;\;\;\left( {若\;{c_{{\rm{vp}}}} \ge {c_{{\rm{sat}}}}} \right)\\ {\gamma _{{\rm{evap}}}}\varepsilon {s_{{\rm{lq}}}}\left( {{c_{{\rm{vp}}}} - {c_{{\rm{sat}}}}} \right)\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\left( {若\;{c_{{\rm{vp}}}} < {c_{{\rm{sat}}}}} \right) \end{array} \right.\;\;\;\;若\;\;T \ge {T_{{\rm{freeze}}}}\\ 0\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;若\;\;T < {T_{{\rm{freeze}}}} \end{array} \right. $ (54)

当孔隙中的局部温度高于冰点时,冰融化成水,源项数值为负;当温度低于冰点时,液态水转化成冰,数值为正.液态水和冰之间转化的源项表达式为

$ {S_{{\rm{l}} - {\rm{i}}}} = \left\{ \begin{array}{l} {\gamma _{{\rm{fusn}}}}\varepsilon {s_{{\rm{lq}}}}{\rho _{{\rm{lq}}}}\;\;\;\;\;\;\;\;若\;T \le {T_{{\rm{freeze}}}}\\ - {\gamma _{{\rm{melt}}}}\varepsilon {s_{{\rm{ice}}}}{\rho _{{\rm{ice}}}}\;\;\;\;若\;T \ge {T_{{\rm{freeze}}}} \end{array} \right. $ (55)

当孔隙中的局部温度低于冰点时,由于不考虑液态水的存在,水蒸气在达到饱和状态后直接凝华成冰,数值为正;对于温度大于冰点时,水蒸气永远不会发生凝华现象.水蒸气和冰之间转化的源项表达式为

$ {S_{{\rm{v}} - {\rm{i}}}} = \left\{ \begin{array}{l} \left\{ \begin{array}{l} {\gamma _{{\rm{desb}}}}\varepsilon \left( {1 - {s_{{\rm{lq}}}} - {s_{{\rm{ice}}}}} \right)\left( {{c_{{\rm{vp}}}} - {c_{{\rm{sat}}}}} \right)\;\;\;\;\;\left( {若\;{c_{{\rm{vp}}}} \ge {c_{{\rm{sat}}}}} \right)\\ 0\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\left( {若\;{c_{{\rm{vp}}}} < {c_{{\rm{sat}}}}} \right) \end{array} \right.\;\;\;\;若\;\;T < {T_{{\rm{freeze}}}}\\ 0\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;若\;\;T \ge {T_{{\rm{freeze}}}} \end{array} \right. $ (56)

当聚合物的水含量高于平衡水含量时,数值为正,表示聚合物上的水会析出变成水蒸气;反之,数值为负,表示水蒸气被聚合物吸收.膜结合水和水蒸气之间转化的源项表达式为

$ {S_{n - {\rm{v}}}} = {\zeta _{n - {\rm{v}}}}\frac{{{\rho _{\rm{m}}}}}{E}\left( {{\lambda _{{\rm{nf}}}} - {\lambda _{{\rm{equil}}}}} \right) $ (57)

当环境温度高于冰点,不会出现膜结合水和孔隙中冰的相变;反之,聚合物中的膜结合水在达到饱和状态以后,会在孔隙中生成冰,如果水含量还没有达到饱和状态则不会析出冰.膜结合水和冰之间转化的源项表达式为

$ {S_{{\rm{n}} - {\rm{i}}}} = \left\{ \begin{array}{l} {\zeta _{{\rm{n}} - {\rm{i}}}}\frac{{{\rho _{\rm{m}}}}}{E}\left( {{\lambda _{{\rm{nf}}}} - {\lambda _{{\rm{sat}}}}} \right)\;\;\;\;若\;{\lambda _{{\rm{nf}}}} \ge {\lambda _{{\rm{sat}}}}\\ 0\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;若\;{\lambda _{{\rm{nf}}}} \ge {\lambda _{{\rm{sat}}}} \end{array} \right. $ (58)

当环境温度低于冰点时,聚合物中的膜结合水在达到饱和状态以后,转化成聚合物上的膜冻结水,如果膜结合水含量还没有达到饱和状态且聚合物上含有冻结水,冻结水会转变成膜结合水;当环境温度高于冰点时,没有膜冻结水的存在,不会出现膜结合水和冻结水的相变.膜结合水和冻结水之间转化的源项表达式为

$ {S_{{\rm{n}} - {\rm{f}}}} = \left\{ \begin{array}{l} {\zeta _{{\rm{n}} - {\rm{f}}}}\frac{{{\rho _{\rm{m}}}}}{E}\left( {{\lambda _{{\rm{nf}}}} - {\lambda _{{\rm{sat}}}}} \right)\;\;\;\;若\;{\lambda _{{\rm{nf}}}} \ge {\lambda _{{\rm{sat}}}}\\ {\zeta _{{\rm{n}} - {\rm{f}}}}\frac{{{\rho _{\rm{m}}}}}{E}{\lambda _{\rm{f}}}\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;若\;{\lambda _{{\rm{nf}}}} \ge {\lambda _{{\rm{sat}}}} \end{array} \right. $ (59)

注意, 在本小节表达式中, ζn-vζn-iζn-f表示相变速率; λnf表示阴阳极催化层和质子交换膜内的膜结合水含量; λf为质子交换膜中的膜冻结水含量。

1.8 边界和初始条件

由于膜结合水的求解区域为催化层和膜,几何边界位于催化层和气体扩散层的交界面,采用第二类边界条件,表达式如下[19]

$ \frac{{\partial {\lambda _{{\rm{nf}}}}}}{{\partial x}}\left| {_{{\rm{CL}}\left\| {{\rm{GDL}}} \right.}} \right. = 0 $ (60)

考虑到膜冻结水仅存在膜聚合物中,其控制方程的边界条件仍然采用第二类边界条件,表示如下:

$ \frac{{\partial {\lambda _{\rm{f}}}}}{{\partial x}}\left| {_{{\rm{CL}}\left\| {{\rm{mem}}} \right.}} \right. = 0 $ (61)

多孔介质孔隙中液态水和冰的求解域为催化层和气体扩散层,故液态水和冰的几何边界在气体扩散层和双极板的交界面至催化层和膜的交界面,第二类边界条件表达式如下:

$ \frac{{\partial {s_{{\rm{lq}},{\rm{ice}}}}}}{{\partial x}}\left| {_{{\rm{CL}}\left\| {{\rm{mem}}} \right.}} \right. = \frac{{\partial {s_{{\rm{lq}},{\rm{ice}}}}}}{{\partial x}}\left| {_{{\rm{GDL}}\left\| {{\rm{BP}}} \right.}} \right. = 0 $ (62)

考虑在低温环境下冷启动过程中,通入电池的的气体相对湿度为0,因此水蒸气在边界的第二类边界条件,表示如下:

$ \frac{{\partial {c_{{\rm{vp}}}}}}{{\partial x}}\left| {_{{\rm{CL}}\left\| {{\rm{mem}}} \right.}} \right. = \frac{{\partial {c_{{\rm{vp}}}}}}{{\partial x}}\left| {_{{\rm{GDL}}\left\| {{\rm{BP}}} \right.}} \right. = 0 $ (63)

由于阴极的水含量较多,电池阴极出口的水蒸气默认处于饱和状态,阳极侧水含量较低,进入阳极流道的流体流动阻力忽略不计,出口的流速和进口的流速一致.因此,燃料电池阴阳极出口处的水蒸气的通量为

$ {{\dot n}_{{\rm{cout}},{\rm{vp}}}} = {\xi _{\rm{c}}}\frac{{IA}}{{4F}}\frac{1}{{0.21}}\frac{{{p_{{\rm{sat}}}}}}{{{p_{\rm{c}}}}} $ (64)
$ {{\dot n}_{{\rm{aout}},{\rm{vp}}}} = \frac{{{\xi _{\rm{a}}}IA}}{{2F{c_{{{\rm{H}}_{\rm{2}}}}}}}{c_{{\rm{vp}}}} $ (65)

其中:ξcξa分别表示阴阳极的化学计量比;I为电池的电流密度;pc为阴极进口压力,压强的影响忽略不计,默认为定值;psat为饱和水蒸气压.

对于能量守恒方程,由于电池会与外界环境发生换热,采用第三类边界条件,换热量表达式为

$ \dot Q = hA\left( {{T_{{\rm{amb}}}} - {T_{\rm{w}}}} \right) $ (66)

式中:Tamb为电池外部环境温度;Tw为电池表面温度,通常指端板周围的环境温度.

根据求解方程的数量,需要设定6组初始条件,分别为:电堆中初始水蒸气含量、孔隙中液态水的体积分数、冰的体积分数、聚合物上的初始膜结合水含量,冻结水含量和初始温度.与三维模型中求解电子和质子传输方程来计算电压损失的方法不同,本模型直接基于巴特勒-沃尔默方程和反应气体低温下的传输特性求解电压损失.

2 结果分析与讨论

为了验证模型的有效性,模型的计算结果与Tajiri等人[18]实验结果进行了对比.如图 3所示,选取金属双极板单电池环境温度分别为-20 ℃和-30 ℃的两个工况,膜结合水的初始值为6.2,膜厚30 um,阴阳极反应气无加湿,启动模式为0~80 s间电流密度线性增长到40 mA·cm-2后保持稳定.从图中可以看出,两种工况下电压的变化趋势模拟结果和实验数据吻合较好,表明建立的模型可以用来预测燃料电池低温环境下启动性能.

图 3 低温冷启动试验与数值模拟对比 Fig.3 Comparison of voltage between simulation results and experiment

考虑到在电池启动过程中,如果燃料电池阴极催化层电化学反应生成的水可以及时排出或者被其他材料吸收,催化层多孔介质孔隙中将腾出更多的空间允许反应气到达催化剂表面,提高冷启动能力.由于本文假设初始阶段电化学反应生成的水以膜结合水形式出现,膜结合水存在于催化层的聚合物中,反应气传质在多孔介质孔隙中进行.显然,如果电化学反应生成的水被催化层和质子交换膜中的聚合物吸收,多孔介质孔隙体积较大,反应气传质通道不会被阻碍,电化学反应便可持续进行.因此,本文提出WSC这一指标,其和聚合物体积分数、多孔介质孔隙率和质子交换膜厚度三者呈正相关性,以上三个参数值越大,WSC值越高,电池内部容纳冰的体积越大,反应气传质得到改善,低温启动能力更佳.因次,本文从质子交换膜膜厚度、催化层聚合物体积分数、催化层孔隙率三个参数出发,研究WSC指标对低温冷启动的影响.由于目前国内燃料电池汽车冷启动温度目标已经下探到-20 ℃左右,本文聚焦-20 ℃低温环境下质子交换膜燃料电池单电池冷启动数值模拟研究.

多孔介质孔隙是传质通道,冷启动失败的直接原因是多孔介质的孔隙被堵,反应气无法源源不断的进入催化剂表面.图 4图 5表示的是在启动温度-20 ℃、初始含水量3.4和启动电流密度从0~80 s线性增长到800 mA·cm2的工况.显而易见,75 s之前,孔隙率对启动电压是没有影响的,因为此时催化层多孔介质还没有冰,电化学反应生成的膜结合水未至饱和.此后,催化层孔隙率越大,容纳冰的能力越强,启动电压越高.图 5中,当催化层孔隙率为0.45时,冷启动时间比0.25孔隙率结冰速率慢,冷启动时间延长近80 s.

图 4 电压与催化层孔隙率的关系 Fig.4 Relationship between voltage and catalyst layer porosity
图 5 冰体积分数和催化层孔隙率的关系 Fig.5 Relationship between ice volume fraction and catalyst layer porosity

质子交换膜的水传递除了电迁移以外,还有反扩散现象.在阴阳极水浓度梯度的驱使下,电化学反应生成的水会从阴极催化层经过质子交换膜反扩散至阳极.本文对膜厚度为0.01、0.02和0.03 mm分别进行了数值模拟研究.图 6~图 8表示的是启动温度-20 ℃、初始含水量3.4和启动电流密度为0.1 A·cm2的工况.从图 6可以看出,质子交换膜厚度越薄,初始阶段启动电压越高,之后最低.原因是膜越薄,阴阳极水浓度梯度越大,反扩散越充分,膜吸水速率更快,源源不断的给阳极的质子电迁移提供水分子来源.对于膜厚0.05 mm的情况来说,冷启动时间比膜厚0.01 mm延长10 s左右.原因是膜的吸水速率和吸水量是不同,膜越厚吸水更多,WSC值更高.

图 6 电压与质子交换膜厚度关系 Fig.6 Relationship between voltage and membrane thickness
图 7 冰体积分数和质子交换膜厚度关系 Fig.7 Relationship between ice volume fraction and membrane thickness
图 8 电压与催化层聚合物体积分数的关系 Fig.8 Relationship between voltage and ionomer volume fraction in catalyst layer

燃料电池电化学反应发生在反应物、催化剂和电解质三相界面,因此,阴极催化层电化学反应生成的水起初位于聚合物上,即膜结合水.当膜结合水达到饱和状态以后,膜结合水才析出到多孔介质里,从而阻碍扩散层中的反应气到达催化层.增加催化层中聚合物的体积分数,可以有效延长膜结合水达到其平衡的时间,从而为低温冷启动争取更多的时间和空间.从图 8可以看出,阴极催化层聚合物体积分数0.4对应的启动电压越高,启动时间更久,冷启动能力更高.原因是在冷启动阶段,反应气进气湿度为零,聚合物体积分数越高,吸收的水分多,对于较为干燥的膜有湿润作用,有利于质子传导,电压输出高.如图 9所示,当聚合物体积分数为0.4时,聚合物吸水饱和时间为25 s, 明显大于其他两种工况,这一点区别于催化层孔隙率对冷启动的影响,同时对应启动时间比聚合物体积分数为0.2时延长近30 s左右.说明聚合物体积分数越大,WSC变高,冷启动能力越强.

图 9 冰体积分数和催化层聚合物体积分数的关系 Fig.9 Relationship between ice volume fraction and ionomer volume fraction in catalyst layer
3 结论

建立了一种新的燃料电池单电池冷启动多相数值模型,通过与实验结果的对比,验证了该模型的有效性.研究了质子交换膜厚度、催化层聚合物体积分数、催化层孔隙率三个因素对燃料电池冷启动影响,最终提出了WSC这一指标来表征低温冷启动能力,得到以下结论:

(1) 由于初始阶段电化学反应生成的水位于阴极催化层聚合物上,阴极催化层聚合物体积分数越大,吸收的水更多,低温冷启动能力提高

(2) 质子交换膜的吸水速率和吸水量对低温冷启动能力的影响不同.膜越薄,吸水快,初始启动电压高,膜越后,吸水量大,低温冷启动能力提高;

(3) 催化层多孔介质孔隙体积越大,容纳冰的能力增强,冷启动能力升高;

(4) WSC储水量这一指标可以用来判断燃料电池低温冷启动能力,WSC值越高有利于低温冷启动.

参考文献
[1]
张剑波, 黄福森, 黄俊, 等. 质子交换膜燃料电池零下冷启动研究进展[J]. 化学通报, 2017, 80(6): 507
ZHANG Jianbo, HUANG Fusen, HUANG Jun, et al. A review on subzero startup of proton exchange membrane fuel cell[J]. Chemistry, 2017, 80(6): 507
[2]
DOE technical targets for fuel cell systems and stacks for transportation applications[EB/OL][2018-03-01] https://www.energy.gov/eere/fuelcells/doe-technical-targets-fuel-cell-systems-and-stacks-transportation-applications.
[3]
许澎, 高源, 许思传. 质子交换膜燃料电池停机后吹扫仿真[J]. 同济大学学报(自然科学版), 2017, 45(12): 1873
XU Peng, GAO Yuan, XU Sichuan. Numerical simulation on gas purge after shutdown in proton exchange membrane fuel cell[J]. Journal of Tongji University (Natural Science), 2017, 45(12): 1831
[4]
YANG X G, TABUCHI Y, KAGAMI, et al. Durability of membrane electrode assemblies under polymer electrolyte fuel cell cold-start cycling[J]. Journal of the Electrochemical Society, 2008, 155(7): B752 DOI:10.1149/1.2926505
[5]
KONNO N, MIZUNO S, NAKAJI H, et al. Development of compact and high-performance fuel cell stack[J]. SAE International Journal of Alternative Powertrains, 2015, 4(1): 123
[6]
MANABE K, IMANISHI H, OGAWA T, et al. Development of fuel cell hybrid vehicle rapid start-up from sub-freezing temperatures[J]. SAE World Congress & Exhibition, 2010, 41: 1379
[7]
李义, 许思传, 许澎. 质子交换膜燃料电池系统低温起动仿真[J]. 电源技术, 2016, 40(6): 1202
LI Yi, XU Sichuan, XU Peng. Simulation of cold start of PEM fuel cell system[J]. Chinese Journal of Power Sources, 2016, 40(6): 1202 DOI:10.3969/j.issn.1002-087X.2016.06.014
[8]
SUNDARESAN M, MOORE R M. Polymer electrolyte fuel cell stack thermal model to evaluate sub-freezing startup[J]. Journal of Power Sources, 2005, 145(2): 534 DOI:10.1016/j.jpowsour.2004.12.070
[9]
KHANDELWAL M, LEE S, MENCH M M. One-dimensional thermal model of cold-start in a polymer electrolyte fuel cell stack[J]. Journal of Power Sources, 2007, 172(2): 816 DOI:10.1016/j.jpowsour.2007.05.028
[10]
ZHOU Y B, LUO Y, YU S, et al. Modeling of cold start processes and performance optimization for proton exchange membrane fuel cell stacks[J]. Journal of Power Sources, 2014, 247(2): 738
[11]
LUO Y, JIAO K, JIA B. Elucidating the constant power, current and voltage cold start models of proton exchange membrane fuel cell[J]. International Journal of Heat and Mass Transfer, 2014, 77(4): 489
[12]
DU Q, JIA B, LUO Y, et al. Maximum power cold start mode of proton exchange membrane fuel cell[J]. International Journal of Hydrogen Energy, 2014, 39(16): 8390 DOI:10.1016/j.ijhydene.2014.03.056
[13]
LUO Y, JIA B, JIAO K, et al. Catalytic hydrogen-oxygen reaction in anode and cathode for cold start of proton exchange memgrane fuel cell[J]. International Journal of Hydrogen Energy, 2015, 40(32): 10293 DOI:10.1016/j.ijhydene.2015.06.094
[14]
MAO L, WANG C Y. Analysis of cold start in polymer electrolyte fuel cells[J]. Journal of the Electrochemical Society, 2006, 154(2): B139
[15]
JIAO K, LI X. Three-dimensional multiphase modeling of cold start processes in polymer electrolyte membrane fuel cells[J]. Electrochimica Acta, 2009, 54(27): 6876 DOI:10.1016/j.electacta.2009.06.072
[16]
MENG H. A PEM fuel cell model for cold-start simulations[J]. Journal of Power Sources, 2008, 178(1): 141 DOI:10.1016/j.jpowsour.2007.12.035
[17]
JIAO K. Experimental and modeling studies of cold start processes in proton exchange membrane fuel cells[D]. Waterloo: University of Waterloo, 2011.
[18]
JIAO K, LI X. Water transport in polymer electrolyte membrane fuel cells[J]. Progress in Energy & Combustion Science, 2011, 37(3): 221
[19]
TAJIRI K, TABUCHI Y, KAGAMI F. Effects of operating and design parameters on PEFC cold start[J]. Journal of Power Sources, 2007, 165(1): 279 DOI:10.1016/j.jpowsour.2006.12.017