0 前言
气体扩散层（gas diffusion layer, GDL）位于PEMFC膜电极（membrane electrode assembly, MEA）的最外侧,通常采用石墨化的碳纸或者碳布制成,起到了增加结构强度以及将反应的氢气和氧气传输至催化层（catalyst layer, CL）参加电化学反应的作用。GDL的孔隙率对燃料电池的气体传质有着重要影响,孔隙率过低,气体传输效率也会随之降低。此外,GDL孔隙率的大小也会影响PEMFC的水管理,当GDL孔隙率过低,积水会堵塞孔隙,导致气体的通透性降低;当孔隙率过高,气体可能会排出过快,导致MEA中大量的氧气难以到达反应位置。
PEMFC运行中,MEA和双极板（bipolar plate, BP）通过直接接触产生回路[3],为降低接触电阻并确保PEMFC单电池的气密性,需要通过热压工艺使BP与MEA充分接触[4]。热压过程中GDL在组装压力下会产生较大变形,从而影响PEMFC的性能和可靠性[5]。BP内流道的结构通常由许多沟槽组成,因而在沟槽和GDL冲压连接时,由于高温高压的影响,GDL容易发生一定程度的不均匀变形。这种变形会影响扩散层内气体的传输。针对GDL压缩的问题,JIAO等[6]通过在BP和GDL间插入不同厚度的金属垫片来控制GDL的厚度,并考察GDL厚度和反应物气体流速对排水效率的影响;LEE等[4]采用有限元法（finite element method, FEM）模拟了金属BP单电池的电池堆装配过程,建立其压力分布的轮廓;魏荣强等[7]建立一个单流道PEMFC三维模型,研究了不同装配压力下GDL厚度变化规律及其对孔隙率和电导率的影响。
PEMFC中反应产水对保持质子交换膜（proton exchange membrane, PEM）的离子电导率是必要的,但大量的水蒸气会稀释反应气体的浓度,减缓反应速率。若液态水不能及时排出,将会堵塞气体通道,阻碍气体进入CL,造成浓差极化加剧。此外,液态水还会覆盖在催化剂表面,从而降低电化学活性面积。以上负面影响都会降低燃料电池的输出电压,这种现象被称为“水淹”[8]。针对GDL的水管理问题,ZHAN等[9]基于一维模型分析了不同GDL结构的液态水相饱和度分布,计算了GDL孔隙中残余水对孔隙度和液体饱和度变化对气体扩散的影响;吴伟等[10]借助多相流格子Boltzmann方法（lattice Boltzmann method, LBM）模拟分析了PEMFC碳纸GDL内的气液两相输运过程,主要研究了GDL疏水性对气液两相流的影响。
在PEMFC中,各种传质、催化和电化学反应等过程非常复杂,仅通过实验进行观察和测量存在误差大、效率低等问题,而且在使用铂等贵金属催化剂和实验运行成本上也具有较大的限制[11]。因此,利用数值模型对电池进行模拟是一种高效且经济的方法。LBM提供了一种可行的框架来模拟电池结构内部的流体传输过程。这种方法通过对复杂结构中的宏观速度、密度和分布函数进行离散化来描述流体的分布特性,并通过碰撞和迁移步骤来演化流体的宏观行为[12],可以处理非线性和非平衡的多相流动问题,并且具有较高的计算效率,目前已经被广泛应用于燃料电池复杂结构及内部流动过程的建模研究,用以指导双极板、扩散层的通道设计。目前LBM多基于单纯的组分扩散模型研究内部物质传输以及两相模型研究内部动态积水过程,鲜有研究采用LBM的多组分模型对其中水蒸气组分的反应生成过程以及传输特性进行模拟,而水蒸气的生成传输过程存在多组分气体的组分变化和竞争效应,并受到燃料电池内部实际工况中多方面复杂因素的影响,包括GDL结构特性、气体反应速率以及GDL内部的积水形成情况等。
1 多相多组分模型
$ \begin{array}{c} f_{k}^{i}\left(x+e_{k} \delta t, t+\delta t\right)-f_{k}^{i}(x, t)= \\ -\frac{1}{\tau_{j}}\left[f_{k}^{i}(x, t)-f_{k}^{i(e q)}(x, t)\right] \end{array}$ (1)
式中：k为不同的速度矢量方向;$ e_{k}$为k方向上的单位向量;$ \delta t$为时间步长;$ f_{k}^{i}(x, t)$为i组分流体在t时刻x处沿第k个方向的粒子分布函数;$ \tau_{i}$为组分i的松弛时间,在扩散问题中为扩散率的函数,其平衡态分布函数为：
$ f_{k}^{i(\mathrm{eq})}=\rho_{i} \omega_{k}\left[1+\frac{u_{i}^{\mathrm{eq}} \cdot e_{k}}{c_{\mathrm{s}}^{2}}+\frac{\left(u_{i}^{\mathrm{eq}} \cdot e_{k}\right)^{2}}{2 c_{\mathrm{s}}^{4}}-\frac{\left(u_{i}^{\mathrm{eq}}\right)^{2}}{2 c_{\mathrm{s}}^{2}}\right]$ (2)
式中：$ \omega_{k}$为速度分布函数在k方向上的权重因子;$ u_{i}^{\mathrm{eq}}$为组分i的平衡态速度;$ \mathcal{c}_{\mathrm{s}}$为格子声速。
各个组分方程中内部的平衡态速度与其共同作用的平衡态速度相关,由式（3）计算：
$ u_{i}^{\mathrm{eq}}=u^{\prime}+\frac{\tau_{i} F_{i}}{\rho_{i}}$ (3)
式中：$ \rho_{i}$为组分i的密度;$ u^{\prime}$为共同速度;对于组分i的作用力$ F_{i}$,其包含流体与其他流体组分间的相互作用、组分与固体壁面的相互作用以及其他的外力三个部分。
每个组分的密度公式是其平衡分布函数求和：
$ \rho_{i}=\sum_{k} f_{k}^{i}$ (4)
式中：$ f_{k}^{i}$为组分i在k方向上的分布函数。
假设各气体组分都是理想气体,即其粒子之间没有相互作用力。在这种情况下,气体分子之间的碰撞过程遵循动量守恒原则。对演化方程求一阶矩：
$ \sum_{i} \sum_{k} \frac{e_{k}}{\tau_{i}}\left[f_{k}^{i}-f_{k}^{i(e q)}\left(\rho_{k}, u^{\prime}\right)\right]=0$ (5)
因此可以得到混合组分的共同速度$ u^{\prime}$：
$ u^{\prime}=\frac{\sum_{i} \frac{\rho_{i} u_{i}}{\tau_{i}}}{\sum_{i} \frac{\rho_{i}}{\tau_{i}}}$ (6)
2 模型验证
通过模拟单相双组分的气体混合过程来验证LBM模型固液边界条件设定的有效性。采用双组分的气体扩散模型,设定模拟区域大小为100 × 100,左侧组分浓度设定为6,右侧组分浓度设定为4,四周的边界设定为反弹边界。将孔隙率φ分别设置为0.8和0.6来验证固体障碍对两种气体组分混合过程的影响。图1为不同计算步数n时的两种组分混合浓度云图,其中x、y分别为模型格点的横、纵坐标。可以看到相同的时间步长下,孔隙率低的结构中组分浓度梯度更大。
Fig. 1 Concentration distribution of gas components for different mixing times with porosities of 0.8 and 0.6图1 孔隙率为0.8和0.6时,不同混合时间气体组分的浓度分布示意图 |
经过12 000个时间步长后,模拟结果显示两种气体组分浓度趋近于5,在模型中选中心点来探究在不同模型中该点的浓度变化,中心组分浓度的变化情况如图2所示。孔隙率达到0.6时两种气体的混合速率比孔隙率为0.8时的混合速率慢,经计算,为在模型中点达到5.3浓度,孔隙率0.8的模型需要经过12 960步长,而孔隙率0.6的模型则需要16 435步长。通过比较可以得出在模型中增大孔隙率会阻碍气体的传质与混合,验证了模型的有效性。
Fig. 2 Trends in concentration of gas components at the center of the two porosities图2 两种孔隙率中心点气体组分的浓度变化趋势 |
3 计算模型
为解决GDL变形以及MEA内部多孔结构的水管理问题,选取一个特征流场,包含气体的进出口通道以及下方的GDL,构建了一个二维的特征模型用于模拟PEMFC中GDL与BP流道的连接处,在特征模型中,只有一个脊存在的区域,并且脊的左侧通道为气体进口,右侧通道为气体出口。二维模型很难准确地模拟出GDL内部碳纤维的复杂结构和条状形态,本文简化碳纤维为随机圆柱分布,二维下模拟纵截面下的碳纤维剖面。假设碳纤维之间相互不重合并按照随机函数在指定区域生成碳纤维。模型构造中在通道中每生成一个纤维截面时计算一次孔隙率,当孔隙率达到设定值之后,停止新增结构。
对于所研究的对流问题,松弛时间$ \tau$可以表示为关于动力黏滞系数$ v$、格子长度$ \Delta x$与格子时间$ \Delta t$的函数关系：
$ \tau=\frac{3 v \Delta t}{(\Delta x)^{2}}+\frac{1}{2}$ (7)
所研究气体流动的格子声速$ c_{\mathrm{s}}=\frac{1}{\sqrt{3}} \frac{\Delta x}{\Delta t} \cdot=\frac{1}{\sqrt{3}} e_{k}$,因此模型中松弛时间的计算公式为：
$ \tau=3 v+\frac{1}{2}$ (8)
由于水蒸气在气体中的比例很小,因此将气体的平均雷诺数$ \overline{R e}$规定为：
$ \overrightarrow{R e}=0.78 \times R e_{\mathrm{N}_{2}}+0.21 \times R e_{\mathrm{O}_{2}}$ (9)
式中：$ R e_{\mathrm{N}_{2}}$为氮气雷诺数;$ R e_{\mathrm{O}_{2}}$为氧气雷诺数。
为了保证雷诺数的恒定,需要选择合适的转换参数,以确保松弛时间取值在0.5 ~ 1之间,从而使算法在合理的计算时间内收敛。经过保证雷诺数不变情况下的单位换算,设定氮气、氧气和气态水的松弛时间分别为0.60、0.63和0.58格子时间。
对于三种气体组分,设定空气的密度为1单位,则三种气体组分在进口处的密度即为其在空气中所占的质量比,分别是氮气78%和氧气21%。对于气态水,根据通常湿度下其在空气中的含量设定为0.1%。
现实中GDL的厚度为50 ~ 500 μm,碳纤维的半径为5 ~ 30 μm,因此本模型重构时设定的基本离散单元格尺寸为1.5 μm;GDL的厚度设定为150格,即225 μm;碳纤维剖面的半径统一设定为5格,即7.5 μm。重构的模型如图3所示。
Fig. 3 Characterization model of the joint part图3 连接处的特征模型 |
模拟通道结构包含两种组分的速度进口以及出口。进口边界采用动力学格式,出口采用纽曼边界,表现为出口处的浓度梯度的固定;左右两边的壁面和固相（流道和碳纤维）的边界采用反弹边界;对于氮气组分在通道底部采用反弹边界,对于氧气则在反弹的同时会有一部分的水蒸气被消耗,同时根据氢氧反应的物质的量比值生成2倍被反应氧气单位的气态水。模拟通道内进口速度保持为0.05格子单位,相同雷诺数下对应的现实气体流动速度为7.5 × 10-8 m/s。
为了确保数值模拟结果不会过分依赖于网格分辨率,将改变孔隙率为0.7时特征模型的网格分辨率,并观察模型出口处的氧气特征浓度比来验证网格无关性并计算网格收敛误差[15]。其中特征浓度比的计算方式如下：
$ \varepsilon_{i}=\frac{C_{i} v_{i}}{\sum_{i} C_{i} v_{i}}$ (10)
式中：$ \varepsilon_{i}$为气体组分i在出口处的特征浓度比;$ C_{i}$为组分i的浓度;$ v_{i}$为组分i的速度在y轴方向的分量,所建立特征模型的出口是y轴向的,因此用速度y轴向的分量和浓度的乘积作为组分的特征浓度。网格收敛误差χ定义如下：
$ \chi=\frac{f_{1}-f_{2}}{f_{1}}$ (11)
式中：$ f_{1}$、$ f_{2}$分别为细网格和粗网格的收敛解。
根据上述计算方法绘制不同模型网格数下的氧气特征浓度比与网格收敛误差绝对值,其中虚框标注部分为模型所使用的网格数,如图4所示。
Fig. 4 Oxygen eigenconcentration ratio and grid convergence error versus number of grids in the model图4 模型中氧气特征浓度比和网格收敛误差与网格数的关系 |
由图4中数据可以看出,在网格分辨率低于设定的情况下网格误差相对较大,而随着网格的细化,网格收敛误差的绝对值基本在0.01上下波动,对氧气的特征浓度比没有显著影响,从而证明了网格无关性。
4 模拟结果与讨论
4.1 孔隙率与扩散层压缩的影响
$ \varphi^{\prime}=1-(1-\varphi) \frac{\delta}{\delta^{\prime}}$ (12)
式中：$ \varphi, \delta$分别为形变之前区域的孔隙率和高度;$ \varphi^{\prime}, \delta^{\prime}$分别为形变之后区域的孔隙率和高度。GDL压缩比 ω 的计算方法如下：
$ \omega=\frac{\delta-\delta^{\prime}}{\delta} \times 100 \%$ (13)
首先考虑没有压缩情况下的孔隙率对组分传递的影响。图5对比了不同孔隙率下三种气体组分的浓度场,可以发现在进口速度和其他初始条件一定的情况下,孔隙率越小,各气体组分在通道内的浓度越大,表明气体被堵塞在通道内,从而影响了传质效率。氧气在右下角的浓度低于其他区域,这一方面是由于空气在模型的右下角浓度偏低,另外一方面是由于氧气作为氢氧反应的反应物在接触底边后会发生反应并且产生消耗。水蒸气是氢氧反应的生成物,其浓度和氧气接触模型底边的浓度有关,因此水蒸气浓度在底边,尤其是模型右下角的部分基本和氧气浓度成反比。综合图6来看,气体出口处靠近双极板壁面一侧的速度要明显高于远离壁面的一侧,这是由于流体总是趋向于走更短的距离通过通道。
Fig. 5 Concentration fields of nitrogen (a-c), oxygen (d-f), and water vapor (g-i) with different porosities图5 氮气（a ~ c）、氧气（d ~ f）、水蒸气（g ~ i）在不同孔隙率下的浓度场 |
Fig. 6 Velocity fields of nitrogen, oxygen and water vapor mixtures with different porosities图6 氮气、氧气、水蒸气混合气体在不同孔隙率下的速度场 |
Fig. 7 Concentration fields of nitrogen (a-d), oxygen (e-h), and water vapor (i-l) at different compression ratios图7 氮气（a ~ d）、氧气（e ~ h）、水蒸气（i ~ l）在不同压缩比下的浓度场 |
Fig. 8 Variation curves of characteristic concentration ratio of oxygen and water vapor at different compression ratios图8 不同压缩比下氧气和水蒸气的特征浓度比变化曲线 |
观察特征浓度比的变化可以发现,氧气的浓度比在压缩比增大之后呈现下降趋势,说明氧气被堵塞在了入口处,无法从出口处排出。水蒸气的特征浓度变化随着压缩比的增加反而有增加的趋势,可能是由于压缩导致气体通道变得狭窄,使得更多的氧气能接触到下部的催化层发生反应。但是由于挤压会降低材料的强度,且由挤压导致的不均匀的扩散层内部分布会导致气体的堵塞。因此,在燃料电池的制备过程中应尽量避免对GDL的挤压和损坏。
4.2 积水液滴的影响
针对存在液滴堵塞现象,面临水管理问题时,将在保证积水液滴面积不变的情况下考察其在GDL底边和在中部两种状况下对气体组分传质的影响。在所建立的特征模型中,尚未涉及双相模型的模拟,因此将液滴设定为球形固体,其边界和碳纤维一样设定为反弹边界。
Fig. 9 Concentration fields of oxygen (a-c) and water vapor (d-f) at different droplet positions图9 氧气（a ~ c）和水蒸气（d ~ f）在不同液滴位置条件下的浓度场 |
Fig. 10 Velocity fields of nitrogen, oxygen and water vapor mixtures at different droplet positions图10 氮气、氧气、水蒸气混合气体在不同液滴位置条件下的速度场 |
通过图11中水蒸气的特征浓度比同样可以发现,位于中部的液滴能够将部分的反应气体引导到催化层处,增加反应气体的浓度,位于底部的液滴则会将催化剂覆盖从而阻碍反应。
Fig. 11 Variation curves of characteristic concentration ratio of water vapor at different droplet positions图11 不同液滴位置下水蒸气的特征浓度比 |
5 结论
通过LBM采用Shan-Chen多组分模型对PEMFC内的GDL与BP流道交界处的物理模型进行模拟,探究了不同孔隙率、不同压缩比对GDL传质的影响,并分析在GDL内部产生积水液滴对电池内部气体的传质和反应的影响,主要结论如下：
（1）在GDL中,孔隙率越小,气体在通道内的浓度越大,造成了气体的堵塞,影响传质效率。而水蒸气的浓度和氧气接触底边时的浓度有关。
（2）扩散层受到压缩变形之后会导致GDL孔隙率减小,气体在流道进口处的堵塞影响气体的传质。但由于压缩导致气体通道变得狭窄,使更多的氧气能接触到下部的催化层而发生反应,在流道附近的反应强度会随着压缩比的增加而增加
（3）位于中部的液滴能够将部分的反应气体引导到催化层处,增加反应气体的浓度,位于底部的液滴则会将催化剂覆盖从而阻碍反应。
