多孔电极中的气水两相微观渗流1)周 平 * 吴承伟 *，2)*（工业装备与结构分析国家重点实验室，大连理工大学运载工程于力学学部，大连，116024）摘要本文基于微流体流动的基本理论, 采用Lattice Boltzmann 方法，对燃料电池气体扩散层多孔介质内的液滴变化过程进行了数值模拟，揭示了微液滴在多孔电极中的形成、长大、聚集和传输机理与规律, 提出仅仅控制疏水纤维所占的比例的传统方法不能从根本上改变液态水在气体扩散层中的传输能力，研究表明，疏水化处理方式也就是疏水纤维的分布形式对液态水传输有决定性影响。

液态水的分布形态主要取决于气体扩散层的微观结构，但在不同润湿性能纤维区域的连接处，其形态一定程度受到纤维润湿性能的影响。

这些结论对研制高性能的多孔电极具有指导意义。

关键词燃料电池，多孔电极，两相流引 言燃料电池是直接将氢气和氧气的化学能转换为电能的电化学装置。

由于最终产物只有水，燃料电池是非常环保的发电装置。

尤其是质子交换膜燃料电池由于结构紧凑、启动快、效率高、无噪音、工作温度低等优点，不仅被认为是现代环保汽车的最佳动力源，而且是潜艇等需要高隐蔽性能等军事武器装备的理想动力源之一。

多孔电极是质子交换膜燃料电池的重要部件之一，它由气体扩散层和催化层组成。

其中气体扩散层是最关键的功能结构材料，它一般是由随机分布的碳纤维或正交编织的碳纤维束组成的高孔隙度材料，厚度在100-300μm左右，具有较高的电导率和渗透性能。

气体扩散层不但可以和催化剂有效结合，而且更重要的作用是将“燃料”（阴极中的氧气或空气和阳极中的氢气）从外部储气设备输送至催化层发生电化学反应并将反应产物（水）和未反应完的气体排出电极。

所有质子交换膜燃料电池多孔电极的理论分析模型都涉及到多孔介质内两相流传输问题。

两相流模型主要用来分析气体扩散层内的液态水平衡问题，也就是通常所说的水管理技术——一项影响燃料电池性能和寿命最关键的技术难题。

过多的含水量会阻碍燃料气体的传输，而当含水量太少，质子交换膜的质子传导能力下降，会导致电极的电化学反应无法进行，如何平衡水的分布是燃料电池正常运转的关键。

传统的宏观均匀化结合非饱和流理论的分析方法高效、简单，过去已经被大量使用。

但是这种基于宏观尺度的分析方法无法揭示两相流在多孔介质内部传输的本质，因而不能指导气体扩散层微结构设计。

由排水困难所经常产生的水淹问题长期得不到很好解决，这一问题已成为燃料电池电极设计最棘手的问题之一。

Paganin等[1]在他们的实验研究中发现15wt% PTFE 含量的气体扩散层可以获得最佳的性能。

Jordan 等[2]通过对燃料电池性能的实验测试比较提出气体扩散层的形态也是决定其性能的重要因素。

气体扩散层内通常都是不同尺寸孔隙共存，并且有着不同的润湿性能[3]。

这些特点在水的传输过程中扮演什么样的角色，如何控制这些参数已达到最佳的性能是亟待研究的问题。

了解扩散层内水传输机理可以用来指导设计气体扩散层的微结构和疏水处理方法。

Lee等[4] 提出扩散层内的孔径分布特性要比孔隙度更为重要，因为微孔的体积控制着气液传输模式。

传统的宏观均匀化方法借助平均孔隙度和渗透系数来分析流体在多孔介质内部的运动[6][7]。

这1)国家自然科学基金（10672035,10721062, 10802019）资助2) E-mail: cwwu@类方法虽然计算效率高，但不考虑微结构，无法分析诸如孔径分布，混合湿润性能等因素的影响。

并且实际上，气体扩散层也不是均匀介质。

对碳布而言，纤维束的特征尺寸大约是200μm ，和整体厚度（约~300μm ）非常接近。

对碳纸而言，最大孔径超过100μm ，也比较接近气体扩散层的厚度[4][3][8]。

要了解非饱和流在气体扩散层内的传输机理，分析模型必须包含微结构有关参数[5]，如：局部纤维表面的润湿性能，孔径分布和纤维直径。

LBM(Lattice Boltzmann method ，格子波尔兹曼方法)是一种新的计算两相流在复杂几何空间流动过程的方法，并且已经尝试在气体扩散层内水分布模拟中应用[9][10][11][12]。

LBM 的基础是细观运动论方法，基本思想是构造简化运动论模型，在微观上反映微观物理本质，并且宏观平均上满足所遵循的宏观方程。

本文尝试采用LBM 来分析气液两相流在气体扩散层微结构内部的传输过程。

LBM 两相流方法包括很多模型，本文所采用的是基于非理想气体状态方程的多相单组分模型，可以模拟相变过程，最早由Shan 和 Chen (SC)[13]提出。

主要研究内容包括液态水从扩散层传输到双极板流道的过程，分析了碳纤维的分布形式、混合润湿性能对液态水分布的影响。

1 LBM 两相流方法1.1 基本理论LBM 基于细观动力学方程，适用于模拟边界条件复杂的介质内的单相和多相流，并且能得到更细微的流动特性。

LBM 从微观粒子尺度出发建立离散的速度模型，在满足质量、动量和能量守恒的条件下，得出粒子的分布函数，然后对分布函数进行计算平均得到压力、速度、密度等宏观量。

其动力学演化方程是：()()()()()b i t f t f t t t f i i i i ......,1,0,,,=Ω+=++x x e x δδ (1)其中，f i 表示沿速度方向i 的粒子分布函数；Ωi (f (x ,t ))是碰撞因子，表示由于碰撞过程引起f i (x ,t )的改变率。

本文采用D2Q9（二维9个离散速度）模型，微观粒子的离散运动方向为：)0,0(0=e (2a)4,3,2,1),2)1(sin ,2)1((cos=−−=i i i i ππe (2b) 8,7,6,5),2)1(sin ,2)1((cos 2=−−=i i i i ππe (2c) 格点处的宏观密度和速度为：∑=ii f ρ (3a)∑=ii i f (3b)c e u ρ在求解过程中，一般把动力方程（1）分成碰撞然后再流动过程中把值传递给相邻格点的对应方向。

Shan 和 Chen 提出的多相LBM 方法采用相邻流体分布函数(PDF) f i 的演化过程：过程与流动过程。

在碰撞过程中计算方程右端项，粒子间微观作用力的概念来实现两相流模拟，作用力的大小由非理想气体状态方程来控制，用这种作用力作为外力项修正速度场。

采用二维SC-LBM 模型来分析单组分多相流系统，在D2Q9模型中，采用BGK (Bhatnagar-Gross-Krook)碰撞项描述粒子()()()()[]t i i i i i τ, i =0,1,2,…,8 (4)其中f i eq 是平衡粒子分布函数，δt 是时间步长，下标i 表示离散速度方向。

有关：. f t f t f t t t f eq ,,1,,x x x e x −−=++δδτ 是松弛时间，与运动粘性t c s δτυ2)5.0(−= 在Shan 和Chen 提出的单组分多相格子波尔兹曼模型中，引入了位于x 和x ’的粒子间的长程作用力： ()()()()()x x x x x x x F −′′′=x ∑′其中ψψ,int G (5)()x ψ被称作有效质量或者势函数，是只跟当前点密度有关的函数：()()()()()x x x ρρψp c s −=22. (6)非理想气体压力p 由状态方程得到。

在S-C 模型中，虚假流量（数值模拟引起的非真实物理速度）是普遍采用不同状态方程时产生的虚假流量，并提出使用P-R 状态方程可存在的，并且无法与真实流量区分。

Peng and Schaefer [14] 讨论了以使虚假速度下降，以保证计算的收敛性。

用P-R 状态方程得到的非理想气体压力为： ()2221ραρT a RT p −= 21ρρρb b b −+− (7) ()()()T T ,12−[]c c c c cp RT b p T R a T /08778.0 /445724.0with 26992.054226.137464.01222==×−++=ωωα (8) 其中T c 使临界温度，p c 是临界压力，ω 是偏心因子，对于水取值0.344[14]。

对于流固界面的润湿现象可以通过给固壁面指定一个赝势来模拟流固界面的相互作用力，相当赝势ψw 从气相有效质量间存在的相似性。

就本文的分析而言，最终要分析的是多孔介法对气体扩散层内两相流流动模拟的可行性，考虑到液态水在气体扩散层内的节用LBM 两相流模型分析了单管毛细驱动流，并与解析解作了比对。

于指定了界面张力[15]。

固壁面从完全疏水到完全亲水的变化可以将固壁面的的ψg 到液相有效质量ψl 调整来实现。

通常的LBM 模型都采用无量纲参数进行分析，本文的所有参数都是无量纲的。

用LBM 来分析实际流场时要借助于模拟模型和实际流场质内部液态水在表面张力作用下的准静态流动过程，其中固体表面润湿性能和多孔介质微结构是决定这一过程的主要参数。

LBM 两相流模拟流动和实际物理流动的参数很难做到完全相等，如在LBM 模拟中采用的密度比和气态水的扩散速度和实际物理量是有差别的。

在分析中对不真实的密度比可能产生的误差和如何消除气态水扩散速度的影响都作了分析说明。

分析中的边界条件都采用了无滑移边界条件。

1.2毛细驱动微流体流动模拟为了验证以上方主要驱动力是毛细压力，因此本根据Washburn 方程，毛细流的动力方程描述的是五种力的平衡，包括惯性力、表面张力、静水压力、粘性力和重力。

动力方程为[16][17]：()()()()()()()()()()−−+⎣t h t h d t h t h g d l l c 2out 22ρ (9) ′+=⎥⎦⎤⎢⎡′+″+p t h m t h x t h l 23)cos (2μθσρ其中，μl 是液体的动力粘度，x 和m 分别是惯性力和的动能压头的修正因子。

Szekely 等[16]和Batten[17]建议x 取7d/12, m 分别为2.45和3.41。

在本文研究中，采用Batten 的建议，也就是x =7D m=3.41。

/12 ，图1 液体毛细上升示意图毛细流动LBM 分析模型如图1所示，整个方形计算区域划分成301×400个格子，中间两垂直平行板分别被划分成5×240个格子，这些被进行计算的，两平行板间的的距离为51个格子。

计算温度设为0.8T c ，静态接触角35º。

动态接触角是毛细驱动流的的一个重要特征。

当平行板占据的格子是不考虑动态接触角时，如图2所示，模拟结果与Batten 修正的Washburn 方程所求得的解析解几乎完全一致。

图2 LBM 数值模拟结果和解析解的比较2 气体扩散层微液滴的形成、长大与传输过程模拟上面一小节是关于微流体多相流动数值模拟方法的研究，本小节将这种方法用于燃料电池气体扩散层微流体流动的数值模拟。