PEMFC带沟槽气体扩散层内传输特性孔隙网络模拟

黎方菊，吴伟，汪双凤

（华南理工大学传热强化与过程节能教育部重点实验室，广东广州510640）

引 言

质子交换膜燃料电池（proton exchange membrane fuel cell，PEMFC)比功率高、稳定性好、零排放、能在低温下快速启动，是清洁高效的能源转换装置，但水管理问题仍然是限制其发展的技术难题之一[1-5]。气体扩散层（gas diffusion layer，GDL)是燃料电池的重要多孔部件之一，发挥着传输反应气体、排除多余液态水、传导电子和热量以及支撑膜电极的重要作用[6-7]。GDL 材料结构对燃料电池水管理有重要影响，因此，优化GDL 以改善燃料电池性能受到广泛关注，孔结构改性（射孔，开沟槽）是常用方法之一[8-9]。Gerteisen 等[10]测试了激光射孔GDL和完整GDL构成的燃料电池性能，发现GDL改性后，燃料电池性能提高。Haußmann 等[11]在实验中发现GDL 射孔尺寸与燃料电池功率密度之间存在相关性，并找到一个最佳射孔直径。Wang 等[12]建立三维数值模型研究了射孔直径和间距对GDL 传水性能的影响，结果显示，射孔直径为100 μm，间距为2 mm 时，GDL 传水性能最优。Fang 等[13]采用LBM方法预测了射孔GDL 厚度方向和平面方向的渗透率、扩散率和热导率以及这些性质随饱和度的变化。Niu 等[14]采用VOF 模型研究了GDL 射孔深度和直径对传水性能的影响，发现射孔能降低GDL 内的液位，且直径为100 μm、深度为100 μm 时效果最显著。目前，对带沟槽GDL 的研究较少，Nishida 等[15]在实验中发现使用带沟槽的GDL 时，燃料电池电压比使用完整GDL 时更稳定。Lee 等[16]采用孔隙网络模型研究了GDL 的工程孔道对液态水分布的影响，结果表明，垂直和水平工程孔道均可降低GDL 平均饱和度，当两者都存在时，GDL 平均饱和度降低约30%，但未研究工程孔道参数对氧气和液态水传输性能的影响。沟槽相关参数的研究结果是优化沟槽设计的重要依据，实验研究需借助昂贵的仪器且操作复杂，因此，数值模型研究GDL 沟槽参数对氧气和液态水传输的影响至关重要。本文通过Matlab编程建立带沟槽GDL 的三维孔隙网络模型，探究沟槽位置、沟道深度、沟槽宽度对液相渗透率和氧气扩散率的影响，并从介观尺度分析沟槽的作用机制。

1 数学模型

1.1 孔隙网络模型

孔隙网络模型采用“喉道”表示小空隙，小空隙交汇处的大空隙用“孔隙”表示[17-18]。三维成像技术和数值随机重构是构建孔隙网络的两种方法，实验方法成本高且操作复杂，从随机重构GDL 中提取拓扑孔隙网络需占用大量计算资源和时间[19]，因此，本文采用既包含空间微观信息又高效快捷的规则孔隙网络来代表气体扩散层。

如图1 所示，将气体扩散层的空隙抽像为正方体表示的“孔隙”和长方体表示的“喉道”相互连接。GDL 计算域包括一个气体通道和两个1/2 肋，气体通道宽（Lr）800 μm，肋宽（Lc）800 μm，总宽度L为1600 μm，展向尺寸Ly为1000 μm，GDL 厚（δc）250 μm。GDL入口采用均匀压力边界，四周及肋下采用无渗透边界。模型中孔隙和喉道的尺寸遵从截断正态分布R～N(μ，σ2），且R∈(a，b)，具体参数如表1 所示，孔间距为25 μm。考虑到GDL材料的各向异性，垂直方向喉道径比平面方向稍大，构建的完整GDL孔隙率为63.7%。此外，本文GDL 中沟槽处孔隙直径设置为与孔间距相等，以保证该区域全部为空隙。

表1 孔隙和喉道尺寸分布Table 1 Pore and throat size distribution

1.2 氧气和液态水传输理论

（1）液态水入侵路径预测 实际工作中GDL的毛细数较低(～10-8)，液相侵入GDL 的过程为毛细力主导的指进过程[16,20]，因此，采用渝渗侵入算法[21]模拟准静态水入侵是合理的。入口喉道与液态水库相连，液态水侵入气液界面处具有最小毛细压力阈值的元素，直至液态水到达出口边界（突破）或达到设定的水饱和度。毛细压力阈值由Yang-Laplace公式[22]表示

式中，γ为表面张力，N·m-1；θ为接触角，（°）；Rc为矩形毛细管内切圆半径，m。

图1 带沟槽气体扩散层的规则孔隙网络模型Fig.1 Regular pore network model of GDL with groove

（2）氧气扩散率的计算 氧气传输过程中任意孔隙i遵守质量守恒

式中，qij表示从孔隙j到孔隙i的扩散通量，mol·s-1；gij表示连接孔隙i和孔隙j的扩散传导率，m3·s-1；Ci、Cj为相邻孔隙i和孔隙j的氧气浓度，mol·m-3。

考虑到Knudsen 效应，氧气扩散传导率采用式(3)计算[23]

式中，ld为毛细管长度，m；d为孔径，m；Db、DK,O2和DO2分别表示氧气在空气中的扩散系数、Knudsen扩散系数和氧气在孔隙中的扩散系数，m2·s-1；Rg为气体常数，J·mol-1·K-1；T为气体温度，K；MO2为氧气摩尔质量，kg·mol-1。

孔隙i和孔隙j之间的扩散传导率为两半孔隙与喉道的串联扩散传导率[24]

给定进口氧气浓度Cin和出口氧气浓度Cout后计算孔隙网络中的浓度分布，可得总氧气通量J（mol·s-1），氧气有效扩散系数采用Fick定律计算

式中，A为与物质传递方向垂直的平面面积，m2；De为氧气有效扩散系数，m2·s-1；采用De/Db表示氧气有效扩散率。水侵条件下，气相在液相中的传输系数非常小，所以当孔隙被液态水侵占时，将与该孔隙相连喉道的扩散传导率设为零。设置不同水侵结束条件，可得液态水突破GDL 时的氧气有效扩散率以及氧气有效扩散率随饱和度的变化趋势。

（3）液态水渗透率的计算 液相绝对渗透率的计算过程与氧气有效扩散系数的计算过程类似[25]，此时的液相传导率为

式中，gh为孔隙或喉道的渗透传导率，m3·Pa-1·s-1；η为液态水黏度，Pa·s。

液相绝对渗透率采用Darcy公式计算

式中，Pin和Pout为GDL 入口压力和出口压力，Pa；Q为液态水总流量，m3·s-1。

液相相对渗透率的计算需先对GDL 执行水入侵，当液态水饱和度达到设定值时，将被液态水侵占的孔网络提取出来，可计算不同水饱和度下的有效渗透率Kl，并获得相对渗透率（Kr=Kl/Ka）随饱和度的变化曲线。

单个孔隙和喉道的尺寸大小是随机的，为消除随机性的影响，每个案例循环计算20 次，取其平均值。此外，本文在计算中还做了如下假设：①水侵开始时，气相充满孔隙网络；②孔隙网络均匀润湿，且接触角为110°；③GDL 入口处液态水压均匀分布；④不考虑电池中的热效应，电池温度为40℃；⑤水和气相为不可压缩流体，且两者不互溶。

1.3 模型验证

毛细压力曲线反映了多孔介质的孔隙分布，对氧气扩散和液态水传输路径具有重要影响。本文计算了完整GDL 的毛细压力分布曲线（具体算法见文献[26]），并与实验结果[27]比较。如图2 所示，在较大饱和度范围内模拟结果与实验结果具有较好一致性，饱和度较低时，误差较大，这受规则孔隙网络模型的孔喉尺寸分布限制，但总体误差在可接受的范围内。

图2 孔隙网络模型验证Fig.2 Validation of pore network model

厚度方向的有效扩散率反映了液态水堵塞孔隙对氧气传输的影响。在常用的氧气传输模型中，De/Db可以表示为[28]

本文模拟的De/Db比式（10）的值低，比文献[29]中的实验数据高，但与早期Penman[30]的实验数据和Shou等[31]的模拟结果吻合较好。本文只关注沟槽参数对氧气有效扩散率趋势的影响，误差对该变化趋势基本没影响，因此采用该模型是合理的。

1.4 计算案例说明

研究沟槽位置对气体扩散层传输性能的影响时，共设置了6个案例，其相关参数如表2所示。

研究沟槽深度对GDL 传输性能的影响时，沟槽宽度为100 μm，沟槽深度为距GDL 表面100～250 μm（沟槽穿透GDL）。

研究沟槽宽度对GDL 传输性能的影响时，保持沟槽穿透GDL，沟槽宽度从50 μm到250 μm。

表2 不同沟槽位置的GDL参数Table 2 Parameters of GDL with different groove positions

2 结果与讨论

2.1 沟槽位置

2.1.1 沟槽位置对气液相界面的影响 由图3（a）、（d）可知，厚度减小，突破阻力减小，饱和度相同时突破点增多，即贯穿GDL 的水簇增多。有沟槽时，液态水簇在沟槽附近更紧密，其他区域提供了更多干燥的氧气传输通道。图3 中（b）、（c）、（f）虽然沟槽使液态水集中分布于沟槽附近，但仍未增加贯穿GDL 入口和出口的水簇量。图3（e）中GDL 入口和出口侧沟槽相对排布，液态水需突破的GDL 基质层厚度较薄，形成了贯穿GDL 的优先传水路径，液态水分布最为集中。

2.1.2 沟槽位置对突破参数的影响 如图4 所示，沟槽减小了局部突破距离，案例2 和案例3 沟槽尺寸相同，相比同厚度的GDL（案例1），突破压力稍降。案例4 相比于有一个沟槽的GDL，突破压力更低，因为液态水的传输遵从指进模型[25]，沟槽只减少了局部厚度，液态水突破至沟槽或从沟槽突破至GDL 出口的过程中沟槽以外较厚区域的孔隙会被搜索，因此，整体厚度减少时突破路径更短。案例5和案例6 虽然在GDL 入口和出口边界都有沟槽，但是案例5沟槽位于相对位置，突破路径最短，突破压力最低。

图3 液相饱和度为0.2时不同沟槽位置的GDL气液相界面Fig.3 Gas-liquid interface in GDL with different groove positions when water saturation is 0.2

图4 突破和干燥状态下不同沟槽位置的GDL性能参数Fig.4 Performance parameters of GDL with different groove positions under breakthrough and dry conditions

沟槽增加了局部孔隙率，与完整GDL 相比，有沟槽GDL 在干燥状态下氧气有效扩散率更高，如图4 所示。案例2、案例3 孔隙率相同，故干燥状态下氧气有效扩散率相同。案例5 和案例6 孔隙率相同，但两者孔喉分布不同，案例6氧气传输路径更曲折，故干燥状态下氧气有效扩散率比案例5 稍低。案例4 和案例1 的孔隙率和孔喉分布都相同，所以干燥状态下氧气有效扩散率相同。

图5 不同沟槽位置的GDL突破时水饱和度沿厚度的分布Fig.5 Water distribution in through-plane direction of GDL with different groove positions at breakthrough

图4中液态水突破状态下的氧气有效扩散率需结合图5 进行分析。案例2 沟槽位于GDL 入口侧，液态水优先侵入毛细压力低的沟槽，增加了GDL 入口侧与液态水的接触面积，突破状态下，GDL 入口侧的水饱和度升高，氧气有效扩散率降低。案例3沟槽位于GDL出口侧，根据指进模型[25]，液态水一旦突破至沟槽，就不再侵入GDL 入口侧的孔隙，所以靠近GDL 入口侧的水饱和度降低，突破时氧气有效扩散率比案例1高8.9%，比案例2高10.8%。沟槽增大了孔隙率，故突破状态下案例3 比案例4 氧气有效扩散率稍高。案例5 局部突破路径最短，GDL 入口界面水饱和度最低，突破状态下氧气有效扩散率最高。案例6 突破路径更曲折，GDL 入口侧的水饱和度明显高于案例5，氧气传输通道减少，故突破时氧气有效扩散率比案例5低10.6%。

2.1.3 沟槽位置对氧气和液态水传输参数的影响

沟槽位置不同，孔喉分布不同，氧气扩散率随饱和度的变化曲线也不同，如图6 所示。比较案例1～案例4 氧气有效扩散率随饱和度的分布，在较低的饱和度下，案例2 的氧气有效扩散率最大。随着液态水饱和度继续升高，案例2 的氧气有效扩散率低于案例4，而案例3却显出优势，氧气有效扩散率最高，故有沟槽的GDL 比减小整体厚度的GDL 更具有优势。案例5 液态水分布最集中，干燥的氧气传输通道最多，氧气有效扩散率最高。

由液态水相对渗透率的计算方法可知，贯穿GDL 入口和出口的液态水通量才对液相相对渗透率有贡献。沟槽位置对液相相对渗透率的影响见图7，案例2、案例3、案例6 突破阻力稍降，不能增加贯穿GDL 的液态水簇量，故液相相对渗透率变化并不明显。案例5 突破阻力明显降低，形成了贯穿GDL 的液态水通道，液相相对渗透率提高。案例4整体厚度降低，液态水更容易突破GDL，形成更多贯穿GDL的水簇，因此液相相对渗透率提高。

2.2 沟槽深度

图6 不同沟槽位置的GDL氧气有效扩散率随水饱和度变化曲线Fig.6 Oxgen effective diffusivity as function of water saturation in GDL with different groove positions

图7 不同沟槽位置的GDL液相相对渗透率随水饱和度变化曲线Fig.7 Liquid relative permeability as function of water saturation in GDL with different groove positions

2.2.1 沟槽深度对突破参数的影响 由图8 可知，随着深度增加，突破压力逐渐降低，当沟槽深度为250 μm 时（沟槽穿透GDL），突破压力骤降至490 Pa。沟槽未穿透GDL时，液态水突破GDL 过程中需搜索尺寸小的孔喉，故深度增加，突破压力降低，但不会低于GDL 最低入侵压力。沟槽穿透GDL时，沟槽界面与液态水直接接触，沟槽孔径用当量直径表示，液态水直接从毛细压力低的沟槽内突破GDL，故突破压力骤降。沟槽未穿透GDL 时，液态水突破至沟槽后，从沟槽排出GDL，沟槽长度远大于厚度，液态水在完成填充沟槽前突破GDL，随着沟槽加深，局部孔隙率增大，突破时氧气有效扩散率增大。但当沟槽穿透GDL 时，整个沟槽界面与液态水接触，突破时液态水充满沟槽，故突破时氧气有效扩散率与沟槽深度为200 μm时相比稍降。

图8 不同沟槽深度和宽度的GDL突破参数Fig.8 Performance parameters of GDL with different groove depths and widths under breakthrough conditions

图9 不同沟槽深度的GDL水侵毛细压力曲线Fig.9 Water invasion capillary pressure curve in GDL with different groove depths

2.2.2 沟槽深度对毛细压力分布曲线的影响 如图9所示，液态水填充沟槽时饱和度上升，毛细压力不变，有沟槽GDL 毛细压力曲线出现平台，随着沟槽深度加深，平台出现时的饱和度越低，平台所对应的压力越低。特殊情况出现在沟槽深度为250 μm 时（沟槽穿透GDL），此时，采用当量直径表示沟槽孔径，液态水首先侵入沟槽，所以在水饱和度为0时出现平台，且平台对应的压力只有490 Pa。沟槽穿透GDL 时，GDL 入口界面液态水在较低的压力下从沟槽排出GDL，防止液态水在GDL 入口界面聚集，降低反应物传输阻碍。

2.2.3 沟槽深度对氧气和液态水传输参数的影响

由图10可见，沟槽深度增加，局部孔隙率增大，氧气有效扩散率也增加。氧气有效扩散率随饱和度的分布，可分为三段：第一段液态水还未到达沟槽，氧气有效扩散率随饱和度的下降趋势与无沟槽的下降趋势一样；第二段液态水集中填充沟槽，氧气有效扩散率随饱和度增加而降低缓慢；第三段液态水侵入沟槽以外的区域，液态水分布较分散，阻碍氧气传输，氧气有效扩散率随水饱和度的增大快速降低。

图10 不同沟槽深度的GDL氧气有效扩散率和液相相对渗透率随液相饱和度变化曲线Fig.10 Oxygen effective diffusivity and liquid relative permeability as function of water saturation in GDL with different groove depths

如图10 所示，沟槽深度为100 μm 和150 μm时，没有形成贯穿GDL 的水通道，对液相相对渗透率的影响并不明显。沟槽深度为200 μm时，沟槽底部与GDL 入口的距离较近，液态水很容易侵入沟槽，增加了贯穿GDL 的水簇量，在一定程度上提高了液相相对渗透率。当沟槽深度为250 μm(沟槽穿透GDL）时，沟槽形成了直接贯通GDL 入口和出口的传水通道，故液相相对渗透率显著提高。

2.3 沟槽宽度

2.3.1 沟槽宽度对突破参数的影响 沟槽宽度是另一个重要参数，如图8所示，完整GDL突破压力为6.82 kPa，沟槽穿透GDL 且宽度只有50 μm 时，液态水从沟槽突破GDL，突破压力下降至907 Pa。突破压力随沟槽宽度增加而下降，但当宽度大于150 μm后，突破压力下降并不明显。没有沟槽时，突破时氧气有效扩散率为0.36，沟槽穿透GDL 且宽度只有50 μm 时，液态水从沟槽内突破，其他区域提供了干燥的氧气传输路径，突破时氧气有效扩散率提高到0.45。沟槽宽度增大，突破时液态水占的体积越大，供氧气传输的干燥GDL 相对减小，所以突破时氧气有效扩散率随宽度增加而降低。

图11 不同沟槽宽度的GDL水侵毛细压力曲线Fig.11 Water invasion capillary pressure curve of GDL with different groove widths

2.3.2 沟槽宽度对毛细压力分布曲线的影响 由图11可知，沟槽宽度增大，当量直径增大，毛细压力降低，毛细压力曲线平台所对应的压力越低，平台结束时的液态水饱和度越大。液态水突破GDL 后，被气体吹扫出燃料电池，根据文献[32]的论述，界面液态水压与气相压力之差低于5.39 kPa 时，沟槽为主要传水通道，GDL 其他区域提供干燥的氧气传输通道。在中、高电流密度下，由于沟槽具有排水能力，界面液态水压并不会快速升高，所以沟槽提高了燃料电池的稳定性，与文献[15]中的结论吻合。

2.3.3 沟槽宽度对氧气和液态水传输参数的影响

如图12 所示，与完整GDL 相比，有沟槽时氧气有效扩散率明显提高，同时氧气有效扩散率随饱和度的变化趋势发生变化。开始时GDL 入口界面的水进入沟槽，阻碍了氧气直接传输通道，氧气有效扩散率快速下降；此后液态水填充沟槽，氧气有效扩散率下降趋势变缓；填充完沟槽后，液态水侵入沟槽以外的区域阻塞了更多氧气传输通道，氧气有效扩散率随饱和度的升高下降趋势变陡。填充沟槽时，沟槽宽度为150 μm 的GDL 氧气有效扩散率最高。虽然在高饱和度下，沟槽越宽氧气有效扩散率越大，但宽度增加，液态水在气体通道中的覆盖率增大，导致氧气与GDL 的接触面积减小，阻碍氧气传输。而且沟槽宽度增大，结构缺失变严重，热阻增大，GDL 内会出现局部热点。考虑到以上因素和本文的结果，沟槽宽度应不大于150 μm。

由图13 可见，与完整GDL 相比，有沟槽的GDL液相相对渗透率明显升高。可将液相相对渗透率的曲线分为三个区域：区域一表示液态水填充沟槽，形成液态水突破路径，液相相对渗透率从零发生突变；区域二表示液态水侵入沟槽以外的大孔隙，液态水分布较分散，液相相对渗透率随饱和度缓慢上升；区域三表示液态水侵入GDL 的较小孔喉，将分散的液态水连通成许多贯通GDL 的水簇，液相相对渗透率随饱和度迅速上升。区域三中有沟槽和完整GDL 的液相相对渗透率在水饱和度接近1时发生重合，此时，除少数具有较小尺寸的喉道没有被水侵外，其他较大孔隙和喉道均被水淹，所以此时较小喉道的分布对液相相对渗透率起决定作用，而沟槽未改变较小喉道的分布，故沟槽对液相相对渗透率的影响失效。

图12 不同沟槽宽度的GDL氧气有效扩散率随饱和度变化曲线Fig.12 Oxygen effective diffusivity as function of water saturation in GDL with different groove widths

图13 不同沟槽宽度的GDL液态水相对渗透率随液相饱和度变化曲线Fig.13 Relative liquid permeability as function of water saturation in GDL with different groove widths

3 结 论

本文通过建立三维孔隙网络模型研究了带沟槽GDL的传输性能，得到如下主要结论。

（1）沟槽改变了毛细压力分布，降低了局部突破阻力，液态水更易排至沟槽并由沟槽排出GDL，所以液态水集中于沟槽附近，降低了其他区域的液态水饱和度，为氧气提供了干燥的扩散路径，氧气扩散性能提高。

（2）沟槽对厚度方向液相相对渗透率的影响取决于是否形成贯穿GDL 入口和出口的传水通道，沟槽深度较浅时，不能增加贯通GDL 入口和出口的水簇，液相相对渗透率基本不变。沟槽穿透GDL 时，提供了优先的传水通道，液相相对渗透率明显提高。

（3）GDL 入口和出口都有沟槽，且沟槽相对分布时，突破阻力最小，液态水分布最集中，氧气和液态水传输性能最优。沟槽加深，传输性能增强，沟槽穿透气体扩散层时氧气和液态水传输性能达到最佳。沟槽变宽液相传输性能增强，氧气传输性能在低饱和度范围内先增强后减弱。综合模拟结果和实际情况，沟槽穿透GDL 且宽度为150 μm 时GDL氧气和液态水综合传输性能最优。

符 号 说 明

A——与物质传输方向垂直的平面面积，m2

a,b——分别为孔喉最小半径和最大半径，μm

Cin,Cout,Ci,Cj——分别为GDL 进口和出口氧气浓度，相邻孔隙i和孔隙j的氧气浓度，mol·m-3

De,Db,DK,O2,DO2——分别为氧气有效扩散系数，氧气在空气中的扩散系数，Knudsen 扩散系数和氧气在孔喉中的扩散系数，m2·s-1

d——毛细管直径，m

dp——沟槽深度，m

gij,gd,gd,pi,gd,pj,gd,t——分别为孔隙i和孔隙j之间的扩散传导率，孔喉扩散传导率，半孔隙i和j的扩散传导率，连接孔隙i和j的喉道的扩散传导率，m3·s-1

gh——孔喉渗透传导率，m3·Pa-1·s-1

J——氧气通量，mol·s-1

Ka,Kl——分别为绝对渗透率、有效渗透率，m2

Kr——液相相对渗透率

Lr,Lc,ld——分别为气体通道宽、肋宽、毛细管长度，m

MO2——氧气摩尔质量，kg·mol-1

n——孔隙相连的喉道数

Pin,Pout,Pc——分别为GDL 入口和出口压力，毛细压力，Pa

Q——液态水流量，m3·s-1

qij——孔隙j到孔隙i的扩散通量，mol·s-1

R,Rc——分别为孔喉半径、毛细管内切圆半径，m

Rg——气体常数，J·mol-1·K-1

S——液态水饱和度

T——气体温度，K

w——沟槽宽度，m

γ——表面张力，N·m-1

δc——GDL厚度，m

ε——孔隙率

η——液态水黏度，Pa·s

θ——接触角，(°)

μ——孔喉平均半径，μm

σ——标准偏差，μm