数值模拟质子交换膜燃料电池的方法

韩丽辉,于春梅,冯根生

(北京科技大学 冶金与生态工程学院，北京 100083)

燃料电池是将燃料的化学能直接转换成电能的装置。燃料电池具有效率高、无污染、重量轻等优点，被称为21世纪的绿色环保能源，具有很大的发展潜力和应用前景[1]。在燃料电池的研究中，实验室的物理实验和数值计算必不可少，可靠的计算程序不仅为物理实验提供理论指导，还可以优化设计燃料电池、节约实验成本、提高科研工作效率。数值模拟软件Fluent在解决流体计算问题方面应用很广泛，并且该软件具有燃料电池模块，但其质子交换膜燃料电池模块是针对氢氧燃料电池开发的，当用户的燃料电池是甲醇或其他质子交换膜燃料电池时，该模块就不能完全适用，用户必须使用UDF编写自己的燃料电池程序。本文从Fluent数值求解机理、燃料电池问题基本设置、编写用户UDF、程序调试等方面，详细阐述了如何利用Fluent软件来数值求解质子交换膜燃料电池问题。

1 Fluent求解机理

1.1 基本数据结构

了解Fluent数据存取方式及程序求解机理是编写UDF的前提和基础。用户自定义函数需要从Fluent求解器中读取数据，数据的基本载体为网格。每个网格单元称为cell，每个cell都有节点(node)、面(face)、网格中心点(cell center)，见图1。每个网格都由面包围而成。Thread 是一种数据结构，可以用来存储一组网格和一组面的信息。与网格相关的数据类型为Node、face_t、cell_t 、Thread、Domain，通过这些数据类型定位到数据载体。本文中求解燃料电池问题的UDF程序会用到Fluent中已定义的一些基本数据宏[2-4]见表1，如网格中心点坐标，网格体积、流体密度、流体温度、流体压力、混合物组分、变量梯度、多孔介质孔隙度等。 表1中，ND ND为Fluent中全局变量，对于二维问题，该值为2，对于三维问题，该值为3；另外宏C_P(c,t)读取到的压力是表压力(Gauge Pressure)，其绝对压力等于该值加上操作压力，Fluent中全局操作压力变量为op_pres。

图1 Fluent 中基本数据结构

表1 Fluent中部分基础数据宏

1.2 方程求解过程

求解燃料电池问题需要求解自定义标量传输方程，Fluent中以压力为基本项的耦合求解器求解过程如图2所示。在迭代计算外要执行两步运算：一是初始化，即用Fluent初始化界面中的数据来初始化求解变量，之后调用用户自定义的初始化程序，后者数据会覆盖前者数据；二是调用用户自定义边界条件程序(PROFILE UDF)。

求解迭代循环过程：首先从执行用户定义的调整程序(ADJUST UDF)开始，之后耦合求解连续方程和动量方程；最后顺序求解能量方程、组分传输方程、湍流方程、用户自定义标量传输方程；所有方程求解结束后执行物料属性更新程序(用户自定义PROPERTY UDF)；最后检查是否满足收敛标准，不满足进行下一次迭代计算，满足则退出循环计算结束。

了解fluent控制方程求解过程对于用户编写自己的UDF非常重要，根据求解过程，用户就能在合适的位置给自定义内存赋值。用户感兴趣的计算数据是通过自定义内存存储的，自定义内存的赋值位置可以放在程序调整宏 、控制方程源项、物料属性UDF，而不能放在迭代循环外边的初始化及边界条件的UDF 程序中。

图2 Fluent中压力基本项求解过程

2 燃料电池问题的基本设置

2.1 调入氢氧燃料电池模块

通过文本框命令行调用质子交换膜燃料电池模块。虽然Fluent提供的质子交换膜燃料电池为氢氧燃料电池，不完全适用甲醇燃料电池问题，但同属于质子交换膜燃料电池，其共同之处是可以借用的，如燃料电池各层材料属性、水合量扩散系数、水合量在催化层固定值等。

2.2 基本模型的选择

自定义标量UDS及内存UDM。在界面中定义自定义标量UDS及自定义内存UDM个数，同时要保证界面中定义的数量大于UDF程序中需要的数量。在质子交换膜燃料电池中至少需要3个标量，即电子电势、离子电势及水含量。自定义内存个数根据用户需要而定，至少包括体电流密度、电流密度矢量、水含量、电渗拖曳系数、反应热等。

选择层流模型、能量方程及组分传输模型。在组分传输模型中不考虑体化学反应，化学反应是通过标量及源项的设置体现出来。

2.3 流体属性的设置

从Fluent物料数据库中选择甲醇空气燃烧混合物,该混合物包括甲醇、氧气、水、二氧化碳、氮气。按顺序设定流体组分，其索引号从0开始，依次为1、2、3、4等，该索引号在UDF中定义物料组分扩散系数时会被使用。流体属性除了密度、比热、黏度等基本参数外，还需要设定2个参数，且这2个参数对计算结果有很大的影响。一个参数是流体组分扩散系数。在多孔介质结构中，气体分子往往被毛细管壁所阻碍，为了计算这些阻碍的影响，修正扩散通量，一般是通过修正扩散系数实现的，该系数主要影响组分传输方程。在质子交换膜燃料电池问题中计算组分扩散系数时可以使用稀释近似模型[5]：

(1)

另一个参数是自定义标量扩散系数(UDS,diffusivity)。对于电子电势和离子电势来说，其扩散系数即在不同流体区域的电导率，交换膜的电导率与电池温度及水含量关系[6]为

(2)

式中，σmem为交换膜的电导率(Ω-1·m-1)；β为膜的电导率系数；λ为膜中水含量；ω为膜的电导率指数。

扩散层、催化层的电导率为常数，集流板固相区域的电导率需要在固体材料属性UDS diffusivity 中设定。

2.4 区域的设置

共有9个计算域，包括阴阳极集流板(2个)、阴阳极流体通道(2个)、阴阳极扩散层(2个)、阴阳极催化层(2个)、质子交换膜(1个)。其中集流板为固体域，流体通道、扩散层、催化层、膜为流体域。区域属性的设置包括多孔介质、源项及固定值。膜、催化层和扩散层均为多孔介质区域，要选择各自的多孔介质材料。多孔介质主要引起动量损失，区域为多孔介质域时，在动量方程中自动引入动量损失源项。计算该源项需要黏性损失系数、惯性损失系数及孔隙度。当流动是层流流动时，惯性损失可以忽略。燃料电池数值计算实际上是在计算流动问题基础上加载电势传输方程，同时由于电化学反应引起能量方程、组分方程、连续方程的源项发生变化，必须要加载源项(UDF 程序)。流体穿过多孔介质区域从阳极到阴极或从阴极到阳极，这种穿透运动有两种驱动力，一是由于动量驱动的，二是由于浓度差驱动的，为了使两极流体不相混合，消除动量驱动作用的影响，在膜域对流体运动的速度设为固定值0。另外膜的水含量影响膜的导电率及水的电渗效果，所以水含量只存在于催化层和膜，对于其他流体域水含量都需设为固定值0。

2.5 边界条件的设置

阴阳极入口的质量流量、流体温度，混合物中各组分质量比例、出口压力、出口回流气体温度、集流板温度、操作压力都是可变量，可根据情况设定。阴阳极集流板电势或电势梯度如下：

阳极电子电势Øan：

Øan=0

(3)

阴极电子电势Øcat：

Øcat=Vcell

(4)

膜离子电势Øm：

(5)

式中n为膜厚。

2.6 求解控制及初始值的设置

压力与速度耦合求解采用SIMPLE算法，采用固定循环模式，即按变量顺序求解，使用BGSTAB方法增加求解过程的稳定性。给自定义标量电子电势和离子电势赋予合理的初始值，是计算是否顺利进行的关键，如果电势初始值不合理，自定义标量残差曲线出现振荡或计算出现错误而停止。合理的电势初始值能够使得阴阳极都有合理的过电势，即电子和离子运动的驱动力，这样才能使电化学反应连续发生。一般离子电势的初始值为负数，如-0.1。

3 编写燃料电池程序

在Fluent软件提供的动量、能量、连续、组分传输等方程的基础上加载电子电势及离子电势标量传输方程，并对原有方程加载源项，并修改混合流体的相关属性，最后数值求解方程组，得到流场、质量组分、电流密度、温度场的分布情况。

3.1 质子交换膜燃料电池基础

质子交换膜燃料电池(PEMFC)结构如图3所示，其阴阳两极各包括集流板(collector)、流道(gas channal)、扩散层(diffusion layer)和催化层(catalyst layer)，阴阳两极中间是电解质膜(electrolyte membrane)[7-8]。燃料从阳极流道入口流入，氧化剂从阴极流道入口流入，阴阳两极入口分别在燃料电池的两个端部，流体相对而行。扩散层和催化层为多孔介质区域，流体扩散到催化层后发生电化学反应，阳极电化学反应产生的氢离子通过电解质膜进入阴极，电子则通过阳极扩散层及集流板和外电路进入到阴极集流板，再经阴极扩散层，进入阴极催化层，在那里和氧化物及氢离子发生阴极电化学反应。阴阳两极不断供应燃料和氧化剂，电化学反应就不断的进行下去，外电路中的负载就得到了电能的供应，这就是简单的燃料电池将燃料化学能转化为电能的过程。

甲醇燃料电池的阳极、阴极催化层电化学反应方程式分别为：

阳极催化层：

CH3OH+H2O=CO2+6e-+6H+

(6)

图3 PEMFC 结构图

阴极催化层：

(7)

阳极催化层电化学反应产生二氧化碳、氢离子和电子，二氧化碳从阳极出口流出，电子从外电路到达阴极，氢离子通过电解质膜到达阴极，这样就形成一个闭合回路[9]。阴极氧气和电子及氢离子反应产生水，水会反向扩散到阳极。在质子交换膜燃料电池中，氢离子的传输离不开水，只有膜的纳米尺度空隙中存在足够的水时，水合氢离子才能够从阳极到达阴极。携带氢离子的水为离子拖曳水，这种现象为电渗拖曳[10]。在甲醇燃料电池中，不仅有水的电渗拖曳，还有甲醇的电渗拖曳和甲醇的扩散，到达阴极的甲醇和阴极的氧反应，产生寄生电流。

3.2 电化学模型

在PEMFC的阴阳极催化层不断发生电化学反应，不断产生和消耗电子和氢离子，反应发生的驱动力是同一网格上固相电势(电子电势)与膜电势(离子电势)之间的电势差，即活化过电势。电子电势传输主要在固相区域中，包括集流板、扩散层及催化层，其传输方程见式(8)。离子电势的传输发生在催化层及膜区域中，其传输方程见式(9)。电势传输方程中的扩散系数即电导率，包括电子在催化层、扩散层、集流板中的电导率，及离子在催化层和膜中的电导率。传输方程中的源项R为体电流密度，即单位时间单位体积产生的电量。体电流的大小与反应物浓度、活化过电势、系统温度、离子交换膜厚度等有关[11-12]。

·(σsolØsol)+Rsol=0

(8)

(9)

式中，σ为电导率(Ω-1·m-1)，下标sol表示固体(阴阳电极)，m代表膜；Ø为电势(V)；R为体电流密度(A·m-3)。

采用简化的Butler-Volmer 公式(即Tafel 公式)来计算阳极和阴极的交换电流密度分别为：

(10)

(11)

式中，jref为参考交换电流密度(A·m-2)；ζ为反应区比表面积(1/m)；[A]和[A]ref分别为反应物浓度及参考浓度 (mol·m-3)；γ为浓度指数；ηan为过电势；α为传递系数；F为法拉第常数(96 485 C/mol)。

在阳极催化层，电子体电流密度等于负的阳极交换电流密度，离子体电流密度等于阳极交换电流密度；在阴极催化层，电子体电流密度等于阴极交换电流密度，离子体电流密度等于负的阴极交换电流密度。阳极阴极的活化过电势不同，有：

阳极催化层：

Rsol=-Ran(<0)，Rm=+Ran(>0)

阴极催化层：

Rsol=+Rcat(>0) ，Rm=-Rcat(<0)

阳极过电势：

ηan=Øsol-Øm

(12)

阴极过电势：

ηcat=Øsol-Øm-Vopen

(13)

3.3 质量源项

流体混合物组分发生变化的地方，需要在连续方程和组分传输方程中增加源项。燃料电池阳极催化层和阴极催化层分别发生电化学反应，从反应方程式(6)和(7)可以看出，阳极催化层消耗掉甲醇和水，产生二氧化碳，阴极催化层消耗掉氧气，产生水，当考虑寄生电流时，阴极还会有二氧化碳产生。阳极催化层消耗的甲醇包括电化学反应消耗、电渗拖曳和甲醇扩散，其质量源项表达式见式(14)；阳极水的变化包括电化学反应和电渗拖曳消耗的水及反向扩散来的水，其质量源项表达式见式(15)；阳极电化学反应产生二氧化碳，其质量源项表达式见式(16)。阴极催化层消耗的氧气包括产生主电流需要的氧气和产生寄生电流需要的氧气，其质量源项表达式见式(17)；阴极水的变化包括主电流反应产生的水、寄生电流反应产生的水、电渗拖曳作用带来的水和反向扩散减少的水，其质量源项表达式见式(18)；阴极产生的二氧化碳由于寄生电流而存在，其表达式见式(19)。在连续方程(mass 方程)中加入变化组分的质量源项。

Dch3ohρ(2Y))

(14)

(15)

(16)

(17)

(18)

(19)

式中，S为质量源项(kg·m-3·s-1)；M为分子量(kg·m-3) ；ρ为混合物密度(kg·m-3) ；Y为组分质量分数；Rp为寄生电流密度(A·m-3)；ndr为电渗拖曳系数；Dλ为水含量扩散系数(mol·m-1·s-1) ；Dch3oh为甲醇扩散系数(m2·s-1) 。

3.4 能量源项

由于电流的产生，使得燃料电池除流道外的各个区域的能量方程中都需要增加源项。在催化层，电化学反应的焓变并没有完全转化成电能，为了克服活化能垒有一部分能量损失掉了，催化层能量方程源项见式(20)，在集流板、扩散层和交换膜中，能量方程源项为导电介质的散热量，见式(21)，式中电流密度为单位时间通过单位面积上的电荷，该量为矢量，在某一方向上分量大小为该方向上电势梯度与电导率的乘积。如果考虑到水蒸气发生相变，那么要在扩散层和催化层的能量源项中增加水的汽化潜热。

Sh=-Δhrect-Ran,catηan,cat

(20)

Sh=I2Rohm

(21)

式中，Δhrect为电化学反应焓变(W·m-3)；I为电流密度(A·m-2)，I=σØ；Rohm为电阻，

3.5 自定义内存

在不同自定义函数之间共用的变量值存储在自定义内存中，因此给自定义内存赋值位置在求解燃料电池问题中极为中重要，要根据求解过程中不同宏的调用顺序，在该内存值在被引用之前为其赋值。在燃料电池问题中，至少需要如下自定义内存：过电势、交换体电流密度、水含量、膜的电导率、水的电渗拖曳系数、电流密度及其在x、y、z方向上的分量，以及在计算过程及结果中感兴趣的变量等。过电势在计算交换体电流密度时使用，在电势传输方程源项中赋值；交换体电流密度在各个方程的源项中都可能用到，所以定义在adjust宏中赋值；水含量在计算膜电导率、水含量扩散系数及水的电渗拖曳系数时都会用到，所以在混合物属性宏中赋值；膜的电导率及水的电渗拖曳系数也在属性宏中赋值；电流密度及其分量只在能量方程源项中引用，且除流道域外各个区域的能量方程均匀源项，所以在能量方程源项中给电流密度赋值即可。在计算过程阳极阴极的电流不断接近，收敛时达到相等，我们既关心电流过程值便于调试程序，也关心计算结果电流的大小，所以在自定义内存中要有电流，电流是催化层所有网格体电流密度的总和，在adjust宏中给电流赋值。

4 程序调试

电化学反应主要通过电势标量传输方程及其他方程增加源项体现出来的。程序调试时，首先不考虑电化学反应，即不求解电势标量传输方程也不增加源项的情况下，确保混合物流场的正确性。要保证阳极入口流入的流体在阳极流道、阴极入口流入的流体在阴极流道，一是要设定膜域中由于压力差引起的宏观流动的速度为零，二是设定合理的混合物扩散系数。待无电化学反应的流场计算正确后，再加载adjust 宏及电势传输方程中源项及质量源项，调试电流密度场，最后加载能量方程中源项，调试温度场。

5 计算结果

开路电压为0.95 V，进料温度353 K，阴阳极参考电流密度为1 000 A/m2，操作压力为101 325 Pa，阳极入口流体质量流量为5×10-6kg/s，阳极入口流体质量流量为6×10-7kg/s，阴阳极电荷传递系数为2，电池扩散层、催化层及膜的气孔率为0.5，阳极浓度指数为0.5，阴极浓度指数为1。流道为单通道时计算得到的电流密度场见图4、温度场见图5、甲醇质量场见图6、水分质量场见图7。

图4 电流密度分布

图5 温度分布

图6 甲醇质量组分分布

图7 水分质量组分分布

6 总结

用Fluent 流体计算软件求解质子交换膜燃料电池问题，关键在于设定正确的自定义电势标量扩散系数(即电导率)、标量传输方程源项(体电流密度)、连续方程及组分传输方程中质量源项、能量方程中能量源项，还要设定合理的自定义内存变量及合适的赋值位置。

[1] 奥黑尔.燃料电池基础[M].王晓红，黄宏，译.北京：电子工业出版社，2007.

[2] 韩占钟，王敬，兰小平.Fluent流体工程仿真计算实例与应用[M].北京：北京理工大学出版社，2004.

[3] 王瑞金，张凯，王刚.Fluent技术基础与应用实例[M].北京：清华大学出版社，2007.

[4] 张建文，杨振亚，张政.流体流动与传热过程的数值模拟基础与应用[M].北京：化学工业出版社，2009.

[5] Sukkee U m,Wang C Y,Chen K S.Computational Fluid Dynamics Modelingof Proton Exchange Membrane Fuel Cells[J].J Electrochemical Society,2000，147(12):4485-4493.

[6] Springer T E,Zawodzinski T A,Gottesfeld S.Polymer Electolyte Fuel CellModel[J].J Electrochemical Soc,1991，138(8):2334-2342.

[7] 熊济时，肖金生，潘牧.质子交换膜燃料电池流场模拟与结构尺寸优化[J].武汉理工大学学报：交通科学与工程版,2009,33(6)：534-537.

[8] 马利军，林才顺，薛方勤，等.质子交换膜燃料电池流场截面设计及三维模拟[J].湿法冶金，2008,27(1):52-55.

[9] 刘桂成，张浩，王一拓，等.液体进料直接甲醇燃料电池工况的匹配及优化[J].电池，2012,42(1):7-10.

[10] 潘科，齐亮，刘洋，等.质子交换膜的电渗系数[J].电池，2008,38(2):75-78.

[11] Kulikovsky A A,Divisek J,Kornyshev A A.Modeling the Cathode Com-partment of Polymer Electrolyte Fuel Cells: Dead and Active Reaction Zones[J].J Electrochemical Society,1999，146(11):3981-3991.

[12] Mazumder S，Cole J V.Rigorous 3-D Mathematical Modeling of PEM Fuel Cells II.Model Predictions with Liquid Water Transport[J].J Electrochemical Soc,2003，150(11):A1510-A1517.