质子交换膜电解池内氧气泡输运过程特性

2021-03-16罗马吉

江苏大学学报（自然科学版） 2021年2期

何旭，罗马吉,3，陈奔

(1. 武汉理工大学现代汽车零部件技术湖北省重点实验室，湖北武汉 430070； 2. 武汉理工大学汽车零部件技术湖北省协同创新中心，湖北武汉 430070； 3. 不列颠哥伦比亚大学化学与生物工程系清洁能源研究中心，加拿大温哥华 V6T 1Z3)

随着对能源环境问题的日益重视，传统能源结构已难以满足社会发展的需要,急需一种新型能源来满足不断增长的需求.其中氢能具有来源丰富、质量轻和能量密度高等优点，被各国专家认为是最有发展前途的可再生能源载体之一.目前我国是世界第1大氢气生产国，已连续7年居世界第1位.主要受价格因素影响，其中超过95%的氢气来源于化石能源[1].

质子交换膜电解池(proton exchange membrane electrolysis cell，PEMEC)制氢技术可以在高电流密度和高压下工作，具有效率高、气体纯度高、绿色环保、能耗低、体积小、安全可靠、可实现更高的产气压力等优点，被公认为在制氢领域极具发展前景的电解制氢技术[2].但是质子交换膜电解池在大规模商业化过程中还有一系列的问题尚待解决，例如PEMEC在工作过程中，膜电极(MEA)产生的氧气通过多孔输运层(porous transport layer，PTL)进入阳极流道，与流道内的液态水形成复杂的气液两相流，会对PEMEC内部性能产生影响.相比单相流，两相流动现象更为复杂，比如自由界面的流体动力学、界面现象及质量传输等[3].

国内外对于PEMEC流道内两相流已有诸多研究.HAN B.等[4]研究了PEMEC内的两相流传输对其性能的影响，发现增加PTL孔隙率或减小其表面接触角将改善PEMEC的性能.S. S. LAFMEJANI等[5]研究了PEMEC流场中气液两相流，结果表明：由于出口流道底部到顶部相表面速度的差异，壁面膜厚度会发生变化.J. O. MAJASAN等[6]对平行流场及蛇形流场阳极流道内气泡运动现象进行了可视化研究，结果表明：较高的液态水流速延迟了泡状流到柱塞流的过渡，形成较小的气泡和较短的柱塞.王昭毅[7]采用可视化方法研究了不同电流密度和不同工作温度下的阳极两相流情况，发现电流密度增大会导致气泡在流道上部聚集形成气柱、堵塞总管.目前，对PEMFC(proton exchange membrane fuel cell)的研究也比较多，陈黎[8]、CAI Y. H.等[9]使用VOF(volume of fluid)方法对PEMFC流道中的单个液滴运动现象进行了分析，对比不同工况下液滴的运动过程.然而，对PEMEC的研究更多关注于流场的两相流流型分布现象，并未对氧气单气泡运动特性进行深入分析.

笔者以PEMEC阳极流道内氧气泡为研究对象，采用VOF方法研究流道内氧气泡在液态水中的运动特性，通过分析气泡脱离体积及流道内容积含气率等参数，研究不同PTL表面接触角、液态水流速、气体生成速率及PTL孔径尺寸对氧气泡运动的影响.

1 模型与计算方法

1.1 计算域及网格划分

模拟计算域如图1所示，其为矩形直通道，由作为PTL表面的底面、上壁面和左右两侧壁面组成.参考课题组所做项目的实际流道尺寸，模型大小设置为1 mm×1 mm×4 mm，网格划分采用六面体网格，数量为30万个左右.对网格无关性进行验证表明，当网格数目增大或者减小15%时，计算结果差异较小，因此认为网格数可行.氧气由位于PTL表面中心的圆孔进入流道，模拟实际情况中MEA产生的氧气穿过PTL，从流道下壁面逸出，液态水由左侧进入流道，其中氧气入口与液态水入口距离为1 mm.为保证库朗数小于1，将时间步长设置为10-6s.

图1 计算域

1.2 VOF计算方法

由于VOF方法能够考虑表面张力和壁面黏附作用，可以精确捕捉两相界面，并且VOF方法追踪的是网格单元内流体体积的变化，而非追踪流体质点的运动，所以计算量小，因此采用VOF方法计算气液两相的体积分数.VOF方法的质量和动量守恒方程分别为

(1)

(2)

式中：ρ为平均密度；u为流体的速度矢量；p为静态压力；μ为平均动力黏度；g为重力加速度；Fvol为由表面张力产生的动量源项.

VOF方法能够实现的条件是2种或多种流体(或相)没有相互穿插，所有相的体积分数在1个控制体积中的和等于1.用α代表控制体积中各相的体积分数，以下标1和2分别代表气相和液相，则有

α1+α2=1,

(3)

(4)

物质属性由控制体积中的每一分相决定.对于两相系统，容积比率平均密度为

(5)

表面张力被认为是动量方程中的1个源项，采用CSF(continuum surface force)模型处理，即

(6)

式中：σ12为表面张力系数；ki为气液界面表面曲率.

ki=·n=·(nwcosθ+twsinθ),

(7)

式中：n为两相界面处的单位法向量；nw和tw分别为壁面处的单位法向量和切向量；θ为壁面处的接触角，该接触角不是壁面本身的边界条件，只是用于调整壁面附近单元表面的法向方向.

1.3 模型假设与边界条件

气液两相入口均为速度入口，出口为自由出流，初始状态流道内充满液态水.壁面为无滑移边界.指定气相为第1相，液相为第2相，接触角等均是针对第2相指定.两侧壁面及上壁面是中性壁面，接触角均为90°.工作温度为70 ℃，水和氧气物性参数设置如下：氧气密度为1.1596 8 kg·m-3，黏度为2.266 4×10-5Pa·s；液态水的密度为977.8 kg·m-3，黏度为4.061×10-4Pa·s，表面张力系数为0.064 35 N·m-1.基准条件下液态水入口流速为0.5 m·s-1，即流量为30 mL·min-1.氧气的摩尔产生率N根据法拉第定律计算：

(8)

式中：I为电流密度；A为膜电极活化面积；F为法拉第常数，F=96 485.3 C·mol-1.

电流密度为2 A·cm-2时，氧气产生率为1.658×10-3kg·s-1·m-2，由于氧气生成速率非常慢，为节约计算资源，提高可视化进程，将氧气产生速率扩大2至3个数量级.QUAN P.等[10]指出，生成速率过小时，增大2个数量级并不会对结果产生显著影响，模拟结果仍可以为PEMEC优化设计提供定性指导.

2 结果与讨论

气泡在流道中的运动过程包含生长、脱离和排出流道3个阶段.气泡在生长脱离过程中受力情况复杂，主要受到浮力、表面张力、绕流曳力、壁面黏滞力、气体动量力及剪切升力的综合作用.其中壁面黏滞力和表面张力阻碍气泡的脱离，浮力、绕流曳力及剪切升力促进气泡从PTL表面脱离.下面将研究气泡在不同工作条件下的状态，对比流道中容积含气率φg(流道中气体体积与流道体积之比)和PTL表面气体覆盖率Ag(PTL表面三相接触线包围的面积与PTL面积之比)的变化，对气泡运动过程进行分析.

2.1 PTL壁面接触角对气泡运动过程的影响

在液态水流速vL=0.5 m·s-1、气体生成速率vg=1.0 m·s-1、PTL孔径D=0.2 mm、PTL表面接触角θ分别为40°、65°、90°、115°、140°时，气泡的生长脱离过程如图2所示，θ=40°时，氧气泡呈球形生长并上浮，形成缩颈与孔口相连.在绕流曳力作用下，气泡沿来流方向向下游生长.气泡浮力及绕流曳力随体积增大而升高，当其合力大于表面张力及壁面黏滞力作用时，气泡脱离，并随液态水一同排出；θ增大到65°至140°时，三相接触线外扩，气泡与PTL表面接触面积增加，表面张力增大，向上合力不足以克服表面张力及黏滞力作用，气泡从孔口脱离后未发生θ=40°时的上浮现象，而是黏附在PTL表面向出口移动；θ=115°时，氧气泡在向下游方向运动过程中沿底面铺展，接触到两侧壁面，由于PTL表面与侧面湿润性的差别，气泡趋向于两侧壁移动，中间形成凹陷；气泡在θ=140°时，氧气以气膜形式向前铺展运动，在运动过程中接触右侧壁面，由于PTL表面与侧面湿润性的差别更大，因此气泡迅速在右侧壁面上展开，沿右侧壁面向出口移动.

图2 不同PTL表面接触角下气泡的生长脱离过程

PTL表面接触角对流道中容积含气率的影响如图3所示，在θ=140°时，φg较小，因为θ=140°时，氧气以气膜形式运动，相比于亲水性壁面，气膜高度较低，因此φg较小.

图3 PTL表面接触角对流道中容积含气率的影响

PTL表面接触角对PTL表面气体覆盖率的影响如图4所示，Ag对PEMEC的性能有着重要的意义，若Ag过大，会导致反应物液态水无法通过PTL到达MEA，降低PEMEC的性能.随PTL疏水性增强，PTL表面三相接触线外扩速度增大，Ag增大.

图4 PTL表面接触角对PTL表面气体覆盖率的影响

气泡脱离时间t1和脱离体积V随PTL表面接触角的变化如图5所示，V与t1随θ的增大而增大，这是由于在气体逸出速度较小时表面张力起主导作用，与W. ABBASSI等[11]的研究结果一致.由于θ对Ag有很大的影响，需要考虑PTL材料的影响，才能准确预测气泡脱离时的时间和体积.由分析可知，PEMEC运行过程中，PTL作亲水处理时对其性能有促进作用，因此，以下对θ=40°进行研究.

图5 气泡脱离时间和脱离体积随PTL表面接触角的变化

2.2 液态水流速对气泡运动过程的影响

在θ=40°、vg=1.0 m·s-1、D=0.2 mm、液态水流速vL分别为0.2、0.5、0.8、1.1 m·s-1时，液态水流速对流道内容积含气率的影响如图6所示.

图6 液态水流速对流道内容积含气率的影响

从图6可以看出：φg随vL增大而减小,在vL=0.2 m·s-1、t=16.7 ms时，φg开始下降，即16.7 ms时，气泡才从流道中排出，时间大幅延长.因为vL越大，气泡受到的绕流曳力也就越大，气泡形变加剧，气泡在体积较小时，由孔口脱离，并随液态水快速排出，排出时间减小，φg越小.

液态水流速对PTL表面气体覆盖率的影响如图7所示，在各液态水流速条件下，Ag在一定范围内上下波动，从整体来看，vL越小，气泡不易排出，Ag也就越大，不利于反应的正常进行.

图7 液态水流速对PTL表面气体覆盖率的影响

在PEMEC运行时，适当增加液态水流速有利于气泡排出，J. O. MAJASAN等[6]研究也发现，增大液态水流速可以有效促进气体的排出，避免气体阻塞流道，对其性能有着积极的意义.但是在实际运行过程中，大的液态水流速会导致大的能耗，并使流场中的压力升高，可能会给MEA性能及密封带来影响，因此应综合考虑其影响，选取合适的流速.

2.3 气体生成速率对气泡运动过程的影响

在θ=40°、vL=0.5 m·s-1、D=0.2 mm、气体生成速率vg分别为0.5、1.0、1.5、2.0、2.5 m·s-1时，气体生成速率对流道内容积含气率的影响如图8所示，vg越大，向流道内逸出的气体越多，φg也就越大，王昭毅[7]可视化研究也发现，气体生成速率增大会导致气泡在流道内聚集，φg增大.

图8 气体生成速率对流道内容积含气率的影响

气体生成速率对PTL表面气体覆盖率的影响如图9所示，vg=0.5 m·s-1时，Ag相对较小，在较大流速时，Ag增大，并且波动较大，这是因为vg增大，流道内的气泡数量增多，各气泡运动相互影响，使运动情况更复杂.

图9 气体生成速率对PTL表面气体覆盖率的影响

气泡脱离时间和脱离体积随气体生成速率的变化如图10所示，随着vg增大，气体动量力增强，气泡生长速度加快，体积快速增大，其所受浮力及绕流曳力也快速增大，因此t1缩短，V反而增大.

图10 气泡脱离时间和脱离体积随气体生成速率的变化

气体生成速率的变化对应于PEMEC实际运行过程中电流密度的变化，从经济性考虑，适当增大电流密度可提高制氢率，降低单位成本，但是会导致流道内的气液两相运动状况更复杂，φg和Ag增大，又会对PEMEC性能带来不利影响.

2.4 PTL孔径对气泡运动过程的影响

在θ=40°、vL=0.5 m·s-1、vg=1.0 m·s-1、D分别为0.1、0.2、0.3、0.4 mm时，PTL孔径对PTL表面气体覆盖率的影响如图11所示.

图11 PTL孔径对PTL表面气体覆盖率的影响

从图11可以看出：Ag随着D增大而增大，这是由于D增大时，气泡颈部较粗，气泡向下游变形生长时，三相接触线在PTL表面包围的面积也就越大.

气泡脱离时间和脱离体积随PTL孔径的变化如图12所示，t1先减小后增大，V随D增大而增大，这是由于气体入口速度一致，D增大导致V增大.体积增大会产生更大的绕流曳力，在D为0.1～0.3 mm时，t1缩短，这与耿俊[12]的研究结果一致；D=0.4 mm时气泡的脱离时间略大于D=0.3 mm，这是由于气泡从孔口逸出时，气泡初始颈部较粗，需要更大的绕流曳力使其脱离，因此缩颈由初始状态至脱离的过程较长，导致D=0.4 mm时气泡的脱离时间延迟.PTL是多孔介质，其内部结构复杂，气孔尺寸也有差异.