PEM电解槽流场的多物理场耦合建模与分析

童灵华,孙雷波,应芳义,姚剑琪,叶夏明,谢长君

(1.国网浙江慈溪市供电有限公司,浙江 慈溪 315300;2.国网浙江省电力有限公司宁波供电公司,浙江 宁波 315100;3.武汉理工大学 自动化学院,武汉 430070)

0 引言

氢储能具有能量密度高、容量大、寿命长等特点,是实现碳达峰、碳中和目标的重要大规模储能技术[1-4]。作为氢储能技术的重要设备,质子交换膜水电解槽(proton exchange membrane water electrolyzer,PEMWE)通过聚合物薄膜传导质子,对功率波动反应迅速,与风电、光伏具有良好的匹配性[5]。同时,PEMWE电流密度较高,可显著提高电解效率,发展前景良好。

PEMWE的内部结构和运行机理复杂,其性能受多种因素影响,许多学者对其数值建模进行了研究。Gómez等[6]总结了过去十年的PEMWE数值建模研究,认为多数模型都是基于半经验或经验方程构建,高度依赖假设和近似条件,考虑电解槽动态操作的模型较少。Nafchi等[7]提出一种一维数值模型的建立方法,比较了工作温度、阴极压力、膜厚度、通道宽度和高度、电流密度等不同参数对PEMWE性能的影响;Toghyani等[8]采用田口法降低PEMWE的输入电压;李建林等[9]采用“组件-单体-阵列-系统”的模型层级分析4种数值模型的优缺点和模型的适用性,并提出PEMWE未来数值建模的重点研究方向;刘承锡等[10]提出一种可并网的大容量电解槽动态模型,具备为电网提供快速频率响应的能量。

然而,PEMWE经验建模并不能完全模拟PEMWE内部运行过程,特别是流体流动、热量传递等难以简化表述。近年来,基于多物理场耦合的计算流体动力学(computational fluid dynamics,CFD)模型逐渐应用于PEMWE建模。Kaya等[11]利用CFD模型研究了不同负极催化剂对PEMWE性能的影响,发现Pt-Ir催化剂作为负极催化剂时产氢速率更高。何旭等[12]利用CFD模型分析了多孔运输层表面接触面积、进水流速、氧气生成速率及孔径大小对氧气和氢气泡的生成速率的影响。林楠等[13]建立了一种一维多物理场耦合PEMWE模型,验证了其电化学特性和传质传热特性。

在PEMWE三维CFD建模方面,Zhang等[14]建立了一种单通道PEMWE三维CFD模型研究水流向和流道宽度对电解槽各项指标的影响,认为减小通道深度和增加通道宽度可以改善水传热性能。Khatib等[15]建立了3种不同几何形状的电解槽三维CFD模型,探究在质子交换膜中使用开孔细胞泡沫材料对电解槽性能的影响。一些学者对异形流场的三维CFD建模开展相关研究,Ojong等[16]创建了一个具有三维多孔传输层的半经验全耦合模型,模拟电解槽单电池在多孔金属棉流场中的反应过程。Park等[17]对比通道流场和多孔运输层(PTL)流场的性能,发现通道流场性能优于PTL流场。Olesen等[18]使用ANSYS对圆平面高压PEMWE阳极流场中的两相流进行多物理场耦合分析,发现减少通道弯曲可以提高电解槽的性能。

在当前PEMWE三维建模研究中,以局部组件建模和单一流场建模为主,缺乏多种流场结构的对比分析,尚未明确不同流场结构与PEMWE电化学性能之间的关系,因此有必要通过建立PEMWE的三维CFD模型来研究电解槽运行时内部单电池的传热传质情况,为优化流道设计提供理论依据。综上,鉴于目前研究现状与不足,研究了传热、传质和电化学动力学等多物理场耦合下的电解槽工作特性,建立了PEMWE单电池三维CFD模型,仿真分析了平行流场、单通道蛇形流场、双通道蛇形流场、叶栅形流场等不同流场结构对PEMWE性能的影响。根据物质流速和分布情况,以及电解槽的极化曲线确定PEMWE单电池的最佳流场结构。最后,通过改变温度、压力等条件,验证温度与压力对电解槽性能的影响,发现温度的升高或压强的降低都能有效提升电解槽的性能。

1 模型建立

PEMWE单电池由双极板(bipolar plates,BP)、气体扩散层(gas diffusion layers,GDL)、催化剂层(catalyst layers,CL)和质子交换膜(polymer exchange membrane,PEM)组成。图1为PEMWE单电池结构示意图。在阳极,直流电压将水解离成氧气、质子和电子。质子通过质子交换膜迁移,与从阴极侧的外部电路传导的电子结合形成氢气。

图1中的阳极反应方程为

2H2O(l)→O2(g)+4H+(aq)+4e-

(1)

阴极反应方程为

4H+(aq)+4e-→2H2(g)

(2)

1.1 多物理场模型

在分析PEMWE的多相流情况时,主要考虑电化学动力学、流体动力学和热传递过程。

1.1.1电化学动力学

PEMWE单电池的工作电压由开路电压Voc、活化过电压Vact、浓度差过电压Vcon和欧姆过电压Vohm组成,由以下公式表示[19]:

Vcell=Voc+Vact+Vcon+Vohm

(3)

开路电压由Nernst方程得到:

(4)

式中:R为理想气体常数;F为法拉第常数;T为工作温度;pH2和pO 2分别为氢气和氧气的分压;aH2O为电极与膜之间的水活度;E0为电解槽的标准电动势,E0值受外界环境的温度与压力影响,本文中计算式为

E0=1.229-0.9×10-3(T-298)

(5)

活化极化是由电极反应动力学的缓慢性造成的,使用电极表面反应的Butler-Volmer方程,可以将活化过电压用电流密度表示为

(6)

式中:i0表示交换电流密度;αan和αca表示阳极和阴极的电荷转移系数。

活化过电压表达式为

(7)

由于电极表面的反应物被消耗,电极表面附近的反应会与溶液本身产生浓度差,造成电解槽电压的损失,即浓差极化,故浓度差过电压的本质是粒子扩散的外部表现形式。结合Nernst方程和Fick定律,浓度差过电压的表达式为:

(8)

式中:CO2和CH2分别表示膜和电极交界处的氧气和氢气浓度;CO2,O和CH2,O分别表示其对应参考值。

欧姆过电压由离子阻力和电极与电极板本身存在的等效电阻引起。根据欧姆定律,欧姆过电压与电解过程中的电解电流成线性比例关系:

Vohm=(2RBP+2RGDL+2RCL+RPEM+Rin)iA

(9)

式中:A表示膜面积。PB和GDL的欧姆电阻是根据材料的电阻率计算的,而CL和PEM的厚度很薄,其欧姆电阻一般忽略不计。Rin表示离子穿过膜时的电阻,本文中参考文献[20]相关研究计算:

(10)

式中:δmem、σmem和λ分别表示质子交换膜的厚度、导电率和含水量。

1.1.2流体动力学

质量守恒和动量守恒方程用于描述混合物在PEMWE单电池内的传质过程[21]:

(11)

式中:ε表示多孔介质的孔隙率;μ和ρf表示混合物的黏度和平均密度。Sv表示源项,V表示混合物的体积平均速度[22],其表达式为

(12)

式中:αk和Vk分别为物质k的体积分数和速度。

多孔介质中气液流动的各分量的对流和扩散由Maxwell-Stefan方程表示:

▽·(εVkCk)=▽·(ε1.5Dk▽Ck)+Sk

(13)

式中:Sk是物质k的质量源项;Ck和Dk分别是物质k的摩尔浓度和第k项有效扩散系数,而Dk是关于温度和压力的函数。

(14)

(15)

氧、氢和水的扩散系数定义如下:

(16)

1.1.3热传递

PEMWE的热量来源主要有入口进水带来的热量和电流做功时产生的热量,热量损失的主要方式为对流散热、辐射散热、出口产物以及未反应水带走的热量。电解槽的热平衡可以用其能量方程来描述[23]:

(17)

式中:ρeff、Cp,eff、keff表示有效密度、有效热容和有效热导率,表达式为

(18)

式中:ρs、Cp,s和ks分别表示固体区域的密度、热容和导热率;ρf、Cp,f和kf分别表示流体混合物的密度、热容和导热率。ST是可以通过计算的源项,主要包括不可逆热Qirrev、电化学反应产生的熵热源Qs、双极板多孔介质和质子交换膜产生的欧姆热源Qohm,表达式为:

(19)

多物理场耦合模型中的参数见表1。

表1 理论模型的参数

1.2 三维集合模型

以平行流场为例,图2为本文中应用的PEMWE的几何模型示意图,其结构由上到下依次为阳极流场、阳极电极(包括气体交换层和催化剂)、质子交换膜、阴极电极。

图2 PEMWE单电池几何模型示意图

如图3所示,对3种常见流场形式搭建物理模型,分别是平行流场、单通道蛇形流场和双通道蛇形流场。另外,设计叶栅形流场,增大流道与气体扩散层的接触面积,探讨其电化学性能。平行流场是市面上普遍使用的流场形状;蛇形流场相较于平行流场更利于反应物的接触,且相对于叶栅形流场具有较少的转角;叶栅形流场具有最大的接触面积。除流场形态外,其他几何参数相同,物理参数如表2所示。电极边长是2 cm,流道的宽度和高度是1 mm,流道间隔是1 mm。

表2 物理参数

图3 不同流场形式示意图

1.3 模型假设与参数

为使实验结果更直观,在不影响电解水反应过程的前提下,假设:

1) 反应物和产物都被认为是不可压缩的液体或理想气体,液态水不会蒸发;

2) 只考虑阳极侧的纯水和氧气混合流体流动,只在阳极侧供水,阴极侧的氢气产生后直接排出;

3) 电极的多孔介质存在各向同性和均匀性;

4) 反应物和产物在PEM中的扩散可以忽略不计;

5) 假设膜含有足够的水,且电导率仅是温度的函数。

2 结果与讨论

2.1 模型验证

为验证模型的准确性,将本文模型通过COMSOL Multiphysics仿真计算得到的数据与Weiβ等[24]的研究数据进行比较,结果如图4所示。验证实验选取平行流场PEMWE仿真曲线,温度为25 ℃,压力为0.1 MPa。可以看到,实验数据与仿真曲线基本拟合,误差小于3%,证明了模型的有效性。以该模型为基础,讨论不同流场对PEMWE性能的影响。

图4 PEMWE单电池极化曲线(T=25 ℃,P=0.1 MPa)

2.2 流场仿真验证

2.2.1流速分布

图5展示了在电流密度为2 A/cm2的情况下,不同流场中电解槽的流速分布情况。可以看出,流场的入口处流体流速最大,随后流速沿着通道减小。

图5 不同流场中的速度分布

对比各流场的速度分布,发现平行流场、单通道蛇形流场、多通道蛇形流场和叶栅形流场中的速度最大值分别为3.94、4.56、4.11和4.34 m/s。平行流场与叶栅形流场的速度峰值主要在进出口,在流道内的速度较小,有利于水在阳极充分反应,且进出口较快的流速易于水的流入和氧气的流出。单通道蛇形流场在经过多个转角后流速有明显的下降,且流速大于其他2种流场,这会导致反应物通过多孔介质的迁移较快,虽然较高的流速可以使生成物较快排出,但也会出现反应不充分的情况。而双通道蛇形流场流速分布较为均匀,但在转角处流速从1.5 m/s减小为0.5 m/s,减速情况更为明显。

2.2.2氧气浓度分布

在PEM水电解反应中,氢气和氧气的生成速率主要由电流密度决定,根据法拉第定律可得:

(20)

式中:n为物质的生成速率,由式(20)可知阳极氧气的生成速率与电流密度成正比。氧气浓度分布是析氧速率与混合物流动共同作用的结果。不同流场形态下,更高的析氧速率表明了更好的电化学性能。同时,流场形态会影响氧气在阳极的囤积情况,较高的氧气浓度会抑制反应的持续进行。

图6展示了电流密度为2 A/cm2时不同流场电解槽的阳极氧气浓度分布,各流场从进口端到出口端都有氧气浓度增大的趋势。其中,平行流场水电解池的氧气浓度最大达到了0.63 mol/m3,位于临近出口一侧的中部,说明水进入多孔介质进行了充分反应,整体上有最高的氧气生成率,但平行流场存在由于流速较慢、生成物囤积而使氧气浓度较高的情况;2种蛇形流场的氧气浓度最大达到0.42和0.41 mol/m3,且从进口到出口均匀增大。同时可以看出,叶栅形流场的氧气浓度最大达到0.43 mol/m3,主要位于出口一侧的末端和电极的边缘位置,但并未出现平行流场发生的氧气堆积在流道某一部分的情况。

图6 不同流场中的阳极氧气分布浓度

2.2.3压强分布

根据Navier-Stokes方程中对流体流速与压强关系的描述可知,流速越大则压强越大。结合前文中对流速的分析,对压强分布进行分析。

图7展示了在电流密度为2 A/cm2的情况下不同流场中电解槽的阳极压力,各流场从进口端到出口端都有压力减小的趋势。

图7 不同流场中的阳极压力分布

平行流场的压力减小趋势几乎垂直于进出口流道,且减小幅度较为均匀,相邻流道间的压力差值约为25 kPa。单通道蛇形流场的压力减小趋势平行于长流道,由于单通道蛇形的流速最大,故在此流场形状下的阳极压力整体大于其他3种流道形状,入口处压力接近300 kPa。双通道蛇形流场的压力近似于沿流道方向减小,且每经过一个转角压力降幅约为60 kPa。叶栅形流场的压力沿对角线减小,此流场由于反应较为充分,综合气相和液相的压力减小幅度最大。

2.2.4极化曲线

平行、单通道蛇形、多通道蛇形、叶栅形流场水电解池的极化曲线如图8所示。可以看出,平行流场中的反应物受到流动的限制,导致反应点集中,在电流密度为2 A/cm2时极化电压达1.78 V;2种蛇形流场的性能相差不大,这是因为在中低电流密度时扩散过电位在极化曲线中的占比较小,在电流密度为2 A/cm2时双通道蛇形流场极化电压仅比单通道蛇形流场低0.005 V,这是由于双通道蛇形流场的结构更有利于水在极板通道中流动,不会造成生成物的囤积而导致性能降低;而叶栅形流场中的反应物流动性相对最佳,反应最充分,在电流密度为2 A/cm2的条件下极化电压仅有1.72 V。在其他工作条件相同的情况下,制取相同量的氢气,叶栅形流场的PEMWE消耗的电能最少,即电解槽性能最佳。

图8 不同流场下的极化曲线

2.3 参数敏感性分析

对于同一PEMWE单电池结构,工作温度和堆内压力是影响其电化学性能的主要参数。选取20、40、60和80 ℃这4个不同温度分别计算压强为0.4 MPa条件下的极化曲线,如图9所示。可以看出,当PEMWE工作温度较高时,极化电压更低。这是由于高温能促进反应离子的扩散和向阴极的迁移,因此较高的温度可以增强质子交换膜的传输特性,降低电解槽所需电压。在电流密度为2 A/cm2时,其极化电压分别为1.727、1.721、1.714和1.707 V。

图9 不同温度下的极化曲线

选取0.1、0.2、0.4和0.8 MPa这4个不同压强在80 ℃工作温度下的极化曲线,如图10所示。可以看出,阴极室的压力增大,阻碍氢离子向阴极移动,从而阻碍电解反应的进行,因此电解电压会随着压力的升高而升高。在电流密度为2 A/cm2时,其极化电压分别为1.686、1.696、1.707和1.717 V。

图10 不同压强下的极化曲线

3 结论

1) PEMWE电氢转化率受到流场结构的影响,其中平行流场流速为3.94 m/s,相较于其他流场速度较慢,析氧反应在膜电极中分布不均匀,产生了较大的电化学与浓差过电位,在电流密度为2 A/cm2时极化电压达到1.78 V。在单通道蛇形流场和双通道蛇形流场中流速较快,但过快的流速导致反应不充分,且压力减小较为缓慢。叶栅形流场出口流速较快,为4.34 m/s,反应充分,产生的电动势低于平行流场3.5%、低于蛇形流场2.3%,压力减小幅度最大,性能相对最好。

2) 在叶栅形流场的条件下,温度每升高20 ℃,电解槽极化电压下降0.4%;槽内压力每增大1倍,电解槽极化电压上升0.59%,即升高电解槽温度或降低反应室压力都有助于电解反应的进行,提升PEMWE的性能。