SOFC多孔阳极内微尺度流动的Lattice-Boltzmann模拟

谢 龙，王玉璋，于建国

（上海交通大学动力机械与工程教育部重点实验室，上海 200240）

固体氧化物燃料电池 (SOFC)一次发电效率高，易实现SOFC-GT联合循环，电极材料成本低，燃料适用面广，全固态结构适合进行模块化设计和放大，是中长期最具吸引力的新型、绿色能源解决方案[1]。SOFC的一大发展瓶颈是较高的操作温度要求电池的关键部件必须具有一定的化学和热的相容性，局部热不平衡引起的极大热应力易导致微米级电解质层的密封性问题甚至断裂，使得单电池难以做大[2]。目前主要采用具有低电池运行温度的阳极支撑结构以改善SOFC工作温度高这一技术难点。而在高电流密度或燃料利用率时，阳极支撑SOFC的性能主要受到阳极侧浓度过电势控制[3-5]。阳极中燃料气体以扩散形式供应给电解质层，反应物气体的扩散速率与气体成分和阳极微结构有关[6-7]。已有的对电极微结构研究的局限性在于依据“容积平均”的概念将多孔介质看成一种在大尺度上均匀分布的虚拟连续介质，并采用难以测量（比表面积、迂曲度等）或难以承担特性尺度角色（孔隙率、平均孔隙尺度等）的几何参量来表征其几何结构特性，难以模拟多孔电极内真实的微尺度流动。

本文建立了多孔电极内气体输运的Lattice-Boltzmann(LB)模型，结合多孔介质的数值构造方法，模拟阳极支撑SOFC的微尺度流动过程。在模型有效性验证基础上，进一步分析了不同孔隙率、电流密度、燃料组分及阳极结构对多孔电极内质量传递及浓差极化的影响。模型为多孔电极内气体输运过程的微尺度研究提供了思路，对于制备出性能优良、具有高输出功率密度的电池具有重要的理论及实际工程意义。

1 模型建立

1.1 物理模型

图1以恒定工作温度下的SOFC多孔阳极为研究对象。电极模型的上下侧取为周期性边界；左侧为燃料气体浓度边界(concentration boundary,CB)，燃料气体通过多孔介质流道扩散至右侧三相界面(triple phase boundary,TPB)发生电化学反应；右侧为燃料气体消耗与产物生成的流量边界，固体颗粒表面处理为法向零流量、切向无滑移速度边界。本文燃料选定为H2、H2O与N2的混合气以模拟微尺度多组分流动系统。

1.2 Lattice-Boltzmann模型

LB方法是一种基于统计演化分布函数、研究粒子团碰撞和迁移的介观方法，模拟适用范围从Kn数小于10－3的宏观连续流到101数量级的自由分子区内流动。SOFC的多孔介质电极中Kn≈1，发生分子扩散和Knudsen扩散[9]，同时LB方法具有独特的边界处理方式，使得LB方法对SOFC多孔电极的模拟具有极好的适应性。本文建立的二维等温多组分LB模型由McCracken[10]的不同分子量气体两组分扩散DLS(different lattice speed)模型发展而来。模型演化方程为：

其中，从fαi(x,t)到为演化过程的迁移步，分布函数fαi迁移（或迁移后插值）到临近格点；为演化过程的碰撞步，分为同组分气体内的自我碰撞和异组分气体间的交叉碰撞。

浓度边界与流量边界采用Zou-He格式，固体颗粒表面边界采用非平衡反弹格式[11]。LB模型中的分子扩散由孔隙内分子间的自我碰撞和交叉碰撞体现，而Knudsen扩散则是由分子与固体颗粒的碰撞反映。

1.3 参数选取与多孔介质重构

LB模拟中参数为格子单位，根据文献[11]中物理单位与格子单位的转换规则，模拟1 mm厚的平板式SOFC多孔阳极（孔隙尺度为微米级）时，其厚度方向网格规模将达到103数量级。而根据相似原理，可对计算参数做无量纲处理，通过保证无量纲准则数的一致性来模拟真实流动，从而减小计算规模。选取无量纲准则数为无量纲流量J*=JL/(CTD13)，互扩散系数比D12/D23、D13/D23和孔隙率Φ。研究对象选取为1073 K、0.1 MPa环境下SOFC多孔阳极中H2、H2O与N2多组分理想气体流动系统，物理参数为D12=3.37×10－4m2/s，D13=0.692×10－4m2/s，D23=1.085×10－4m2/s[12]，Φ=0.5，L=0.002 m，TPB 侧J=0.0415 mol/(m2·s)，燃料侧 CT=11.4 mol/m3，燃料气体摩尔比为 H2∶H2O∶N2=0.47∶0.03∶0.5，计算得J*=0.067。物理单位与格子单位转换时需注意：格子单位中孔隙率Φ与物理单位一致；总体网格规模必须足够大才能使多孔结构具有足够的“分辨率”，而计算机有限的计算能力又要求网格规模足够小，本文选取120×120网格以满足两方面的需求，即格子单位中特征长度L=120；当互扩散系数过小时，互扩散松弛时间接近于0.5 s，将会导致数值稳定性问题，因此将计算中格子单位的Dij扩大103数量级；J和CT的取值可调整以满足J*的要求，如J选定为0.069，则此工况下CT值取为1140 mol/m3。由于采用相似原理进行数值模拟，故引入无量纲电流密度i*以进行模拟结果的无量纲分析，同时根据Nernst方程得到浓差极化表达式：

模拟时将多孔介质信息以数字化网格形式读入。首先以二值图形式建立多孔介质，然后将孔隙结构二值图转换为数值0、1形式的数字化网格，即通过标志函数来表示多孔介质的空间分布：

最后将数值矩阵作为数值模拟的网格输入。

需要注意的是，为了实现对多孔电极微尺度流动的完全模拟，还需要引入表征Knudsen扩散与分子扩散比值的无量纲参数D1,k/D1,3（其中D1,k为组分1的Knudsen扩散系数）来保证数值模拟与真实流动的完全相似。一种方法是保持现有网格形式不变（即固定D1,k），通过增大D1,3（即弱化分子间的碰撞，减小连续区域效应）来间接模拟实际多孔电极内的非连续区域效应；另一种方法是修正数值网格（通常是减少孔隙间格子数）以保证网格的Knudsen数与实际一致。本文采用连通情况较好二维多孔介质为研究对象，其内部扩散过程以分子扩散为主，主要关注LB方法模拟多孔电极内微尺度流动的可行性。为了充分模拟多孔电极内的微尺度流动过程，进一步的研究除了引入上述的无量纲扩散系数比之外，还需要以X射线CT扫描等技术[13]重构出的三维多孔介质为研究对象。

2 模型计算与有效性分析

孔隙尺度LB模型的每步迭代计算都需要对网格格点类型进行逻辑判断，并对流体格点与固体边界格点各分布函数进行计算，因此计算量巨大。一种提高计算效率的方法是采用稀疏矩阵，分理出固体格点，计算只在流体格点上进行，不需要进行格点的逻辑判断和搜索，从而减少了存储空间并能提高计算效率[14]；另一种方法是利用LB本身的并行特性，通过并行计算来提高计算效率[15]。本文采用MPI(message passing interface)并行方法计算LB模型，并行算法设计时按一维行分解将计算区域分配给各个处理器，如图2。值得注意的是，若按一维列分解或二维块分解划分计算区域，临近分区互相决定了对方的边界条件，不利于数值稳定性，影响计算精度。采用8核、单核2.5 GHz工作站进行并行计算，图3给出了处理器数量与加速比（串行运算时间与并行运算时间之比），可见随着处理器数的增加，加速效果趋于稳定。

为验证LB多组分流动模型的正确性，采用Yakabe[5]的阳极支撑SOFC性能试验数据进行模拟计算，同时与Suwanwarangkul[16]采用的Dusty-Gas Model(DGM)模型预测进行比较，如图4。可见本模型计算结果与文献所报道结果较为一致，说明LB微尺度模型可有效模拟阳极支撑SOFC的浓差极化。三者之间的差异可能是由于采用的多孔介质以及对多孔介质的描述方式不同，以及计算中未充分考虑Knudsen扩散的影响。

3 结果与分析

图5 稳态工况下多孔阳极内参数分布

图5为稳态工况下多孔电极内燃料气体摩尔分数的典型分布以及系统流线图。H2从燃料侧扩散至TPB参与反应，TPB上H2不断消耗生成H2O，H2O则从TPB扩散至燃料流道；N2在TPB上不参与电化学反应，在多孔电极内表现为梯度较小的浓度差分布；电极的总速度分布则是从TPB指向燃料流道，电极内的电化学反应热与欧姆热不断被带出电极，以达到良好的冷却效果，避免局部过热区的产生。经计算发现各组分流量在电极内各截面上保持恒定，满足质量守恒，因此各组分分布云图中摩尔分数偏小的截面上或狭窄流道中流速较大。

3.1 孔隙率、电流密度对质量传递和浓差极化的影响

图6给出了不同孔隙率与电流密度时TPB上H2的平均摩尔分数分布，可见固定J*时，随着孔隙率的增大，电极内气体扩散通道增加，维持一定摩尔流量所需的浓度差减小；每个孔隙率对应着一定的J*阀值（或极限电流密度i*），该流量阀值表示扩散到TPB处的H2被完全消耗，该阀值随着孔隙率的增大而增加；固定孔隙率时，浓度差随着J*的增加而增大以满足大摩尔流量的需求。

图7给出了不同孔隙率与电流密度情况下多孔阳极的浓差极化特性曲线，可见固定孔隙率时，浓差极化随着i*的增大而增大，固定i*时，浓差极化随着孔隙率的减小而增大。原因是随着孔隙率的减小，H2难以扩散到TPB参与反应且电化学反应生成的H2O难以往外扩散；或者随着i*的增大，TPB上电化学反应速度加快从而H2大量消耗且H2O大量生成，因此TPB上H2的摩尔浓度减小且H2O浓度增大，从而引起较大的极化损失。同时还观察到浓差极化在很大i*范围间呈近似线性分布，而当i*接近该孔隙率对应的极限电流密度时浓差极化出现较大的阶跃，i*恒定时浓差极化随孔隙率减小而增大的趋势也越发明显。因此高i*或低孔隙率情况下浓差极化的变化更加敏感。

3.2 气体组分对浓差极化的影响

图8给出了孔隙率为0.45时，固定燃料气体中N2摩尔分数为50%，变化H2与H2O摩尔比时各i*对应的浓差极化曲线，可见浓差极化随燃料侧H2摩尔比的增大而减小，且除了大的i*或小的H2摩尔体积时表现出的浓差极化陡增以外，其余区域为线性变化。原因是当i*、孔隙率以及N2摩尔分数一定时，多孔电极两侧H2浓度差保持恒定；电极内N2纯扩散不参与反应，在TPB上H2与H2O摩尔浓度之和保持恒定。因此随着燃料侧H2摩尔分数的降低，TPB处H2摩尔分数线性减小，而H2O则线性增加，从而引起浓差极化的近似线性变化。

图9给出了孔隙率为0.45时，固定燃料气体中H2与H2O摩尔比为4∶1，变化H2浓度时各i*对应的浓差极化曲线。i*、孔隙率以及H2/H2O摩尔比一定时，多孔电极两侧H2浓度差随H2摩尔分数的增加而线性减小；而随着N2摩尔分数的减小，TPB上H2与H2O摩尔浓度之和则线性增加。因此燃料侧H2摩尔分数降低会引起TPB处H2摩尔分数的线性减小以及H2O摩尔分数的线性增加，从而引起浓差极化在大H2摩尔体积时的近似线性变化以及小H2摩尔体积时的陡增。

图9 变H2浓度时的浓差极化曲线

3.3 阳极结构对浓差极化的影响及电极结构优化

表1给出了i*=0.064，孔隙率为0.45时不同阳极结构对浓差极化的影响。其中多孔电极结构有二级规则分形Sierpinski地毯结构多孔介质[17]，以等半径圆为基本构造单位的顺排、错排和随机多孔介质，以及真实多孔介质。

顺排、错排与随机三种多孔阳极结构的浓差极化值比较接近说明不同电极结构的主要特性尺度为孔隙率，而浓差极化的细微差异则是由不同迂曲度引起，说明迂曲度为多孔电极的次要特性尺度。分形与真实多孔介质得到了相同的浓差极化值，则说明真实多孔介质具有分形特性。各种形式阳极结构在极化特性上的不同表现说明对多孔结构的合理优化可以减弱电极的浓差极化特性。

表1 不同阳极结构的浓差极化

基于此，对多孔电极结构的优化路线应该从减少浓差极化出发，充分考虑欧姆极化以及活化极化的影响，使极化损失之和达到最小值。i*一定时，孔隙率越大浓差极化越小，而较大的孔隙率会导致电极的结构性能下降，结构整体的欧姆极化损失增加；选择合适孔隙率的同时，多孔结构应满足燃料附近孔隙率较大而TPB附近孔隙率较小，既使得燃料气体可以顺利的往TPB扩散，又保证TPB附近足够的有效反应面积从而减小活化极化损失。

4 结论

(1)LB方法适用于燃料电池多孔电极内多组分气体输运过程的微尺度模拟，通过LB模型可以建立多孔阳极微结构及燃料组分与质量传递和浓差极化的关系。进一步耦合电极内传热过程以及TPB上电化学反应的LB模型将是设计与分析SOFC多孔阳极更为准确有效的工具。

(2)SOFC存在极限电流密度，该阀值随着多孔阳极孔隙率的增大而增加；多孔阳极浓差极化随无量纲电流密度i*或燃料侧H2摩尔比的变化作近似线性变化，而在TPB处的H2被完全消耗时变化剧烈；通过优化阳极多孔结构可以减小浓差极化，以提升单电池性能。

[1]WANG Y Z,YOSHIBA F,WATANABE T,et al.Numerical analysis of electrochemical characteristics and heat/species transport for planar porous-electrode-supported SOFC[J]. Journal of Power Sources,2007,170(1):101-110.

[2]毛宗强，黄建兵，王诚，等.低温固体氧化物燃料电池研究进展[J].电源技术，2008，32(2)：75-79.

[3]王玉璋，于建国，惠宇，等.新型平板式固体氧化物燃料电池的开发和性能分析[J].动力工程，2008，28(6)：949-954.

[4]YI J,ANIL V V.Fuel composition and diluent effect on gas transport and performance of anode-supported SOFCs[J].Journal of theElectrochemical Society,2003,150(7):942-951.

[5]YAKABE H,HISHINUMA M,URATANI M,et al.Evaluation and modeling of performance of anode-supported solid oxide fuel cell[J].Journal of Power Sources,2000,86(1/2):423-431.

[6]VIRKAR A V,CHEN J,TANNER C W,et al.The role of electrode microstructure on activation and concentration polarizations in solid oxide fuel cells[J].Solid State Ionics,2000,131(1/2):189-198.

[7]ZHAO F,VIRKAR A V.Dependence of polarization in anode-supported solid oxide fuel cells on various cell parameters[J].Journal of Power Sources,2005,141(1):79-95.

[8]BIRD R B,STEWART W E,LIGHTFOOT E N.Transport phenomena 2nd EDN[M].New York:Wiley,2002.

[9]于建国，王玉璋，翁史烈.以煤气化合成气为燃料的平板式固体氧化物燃料电池性能[J].中国电机工程学报，2011，30：88-93.

[10]MCCRACKEN M E,ABRAHAM J.Lattice Boltzmann methods for binary mixtures with different molecular weights[J].Phys Rev E,2005,71(4):8.

[11]何雅玲，王勇，李庆，等.格子Boltzmann方法的理论及应用[M].北京：科学出版社，2009.

[12]JOSHI A S,PERACCHIO A A,GREW K N,et al.Lattice Boltzmann method for multi-component,non-continuum mass diffusion[J].Journal of Physics D:Applied Physics,2007,40(23):7593-7600.

[13]LEFEBVRE,PHILIPPE L,BANHART,et al.Porous metals and metallic foams:current status and recent developments[J].Advanced Engineering Materials,2008,10(9):775-787.

[14]郭照立，郑楚光.格子Boltzmann方法的原理及应用[M].北京：科学出版社，2009.

[15]VIDALA D,ROYA R,BERTRAND F.On improving the performance of large parallel Lattice Boltzmann flow simulations in heterogeneous porous media[J].Computers&Fluids,2010,39(2):324-337.

[16]SUWANWARANGKUL R,CROISET E,FOWLER M W,et al.Performance comparison of Fick's,dusty-gas and Stefan-Maxwell models to predict the concentration overpotential of a SOFC anode[J].Journal of Power Sources,2003,122(1):9-18.

[17]王唯威，淮秀兰.分形多孔介质导热数值模拟分析[J].工程热物理学报，2007，28(5)：835-837.