基于多尺度耦合的管式SOFC数值模拟

马 旭，杨 晨

(1. 珠海城市职业技术学院，广东 珠海 519090；2. 重庆大学能源与动力工程学院，重庆 400044)

1 引言

建模与仿真的发展趋势是从计算机宏观计算模拟到粒子团的模拟再到分子级的模拟，即从宏观到介观再到微观。在宏观层次上主要是依据三大守恒定律进行的CFD模拟；在介观层次上近年来研究最多的是基于粒子团的格子玻尔兹曼方法，该方法通过一定的数学变换可以推导出三大守恒定律；在微观层次上主要是依据分子动力学方法开展的模拟。在单尺度仿真中，尽管已经获得了不错的成果，但是仍然存在一些局限性，例如模型无法准确恰当地描述疑难问题，现有算法难以适应不断增长的模型规模和精度要求以及计算机硬件的速度和效率的限制等问题。

多尺度耦合模拟是解决上述问题的有效途径[1]。多尺度模拟在空间尺度上划分出不同的层次，针对问题描述难易程度、模型规模精度、计算速度等问题，使用不同层次的模拟方法和手段[2]，以便能更恰当准确的描述疑难问题，提高模型计算的精度，实现并行计算，提高计算效率。多尺度模拟是仿真发展的趋势，是大型复杂系统有力的描述手段[3]。

本文以管式固体氧化物燃料电池(SOFC)为例，从多尺度模拟的角度，在介观层次上采取格子玻尔兹曼方法模拟电化学反应扩散，在宏观尺度上采取CFD方法模拟管式SOFC内部的传热、流动过程，通过介观层次和宏观层次的耦合实现多尺度模拟。

2 格子玻尔兹曼方法的电化学反应扩散模型

格子玻尔兹曼方法中模拟的是流体粒子团，流体粒子团可以理解为微观充分大、宏观充分小由若干流体分子组成的微团。其具体形式如下[4]

fi(x+eiΔt，t+Δt)=fi(x，t)+Ωi(f(x，t))

(1)

式中，fi是粒子速度分布函数；ei是当地粒子速度；Ωi(f(x，t))是碰撞算子，它表示粒子团发生碰撞后fi的变化率。

格子玻尔兹曼方法中关于电化学反应扩散问题，其基本思想是对碰撞项进行修正，即在碰撞项中加入反应项[5]，具体形式如下

(2)

其中

(3)

(4)

将(2)、(3)和(4)式代入(1)式，利用恰普曼—恩斯科格展开，可得到宏观的反应扩散方程：

(5)

式中，M是格子模型中的速度数；D′为扩散系数；R(ρ)是反应项，形式由具体化学反应决定。

3 管式SOFC的CFD模型

SOFC的基本工作原理如图1。燃料和空气需先压缩并加热到规定温度。燃料被送到SOFC的阳极，空气被送到SOFC的阴极。在阴极室中，随着空气的流动，O2逐渐溶入多孔电极中，在多孔电极中得到电子生成O2-，O2-通过电解质内的离子导体传输至阳极。在阳极室中，CH4与H2O反应生成H2、CO、CO2的混合气体，其中H2、CO与阴极传送过来的O2-发生电化学反应生成H2O、CO2，并产生电子，电子经外部电路被再次送到阴极，如此循环往复，便产生了电能。

图1 SOFC工作原理

过程中涉及到的主要化学反应有[6]

CH4+H2O=3H2+CO

(重整反应)

CO+H2O=H2+CO2

(置换反应)

H2+O2-=H2O+2e

(阳极)

CO+O2-=CO2+2e

(阳极)

0.5O2+2e=O2-

(阴极)

总体反应为

CH4+2H2O=4H2+CO2

CO+0.5O2=CO2

H2+0.5O2=H2O

3.1 质量守恒方程

管式SOFC的质量守恒方程可由下式表示[7]

∇·(ερu)=0

(6)

式中，ε是管式SOFC电极的孔隙率。

燃料电池内部是由CH4、H2O、H2、CO、CO2等气体形成的混合气体，其密度为：

(7)

式中，mi是混合气体中某气体的质量百分比，Mi是混合气体中某气体的摩尔质量。

3.2 动量守恒方程

管式SOFC的动量守恒方程可由下式表示

∇·(ερuu)=-ε∇P+∇·(εμeff∇u)+Su

(8)

式中，Su为源项，在流道中Su数值取零。在多孔电极中发生电化学反应使得气体浓度发生改变，进而引起传质现象的产生。在传热传质过程中，由于流道内的速度分布差异很大，在电极界面将产生很大的速度梯度，因此不能忽视粘性力和惯性力的影响，此时应用布林克曼方程描述动量，源项为

(9)

式中，μeff是气体的有效粘性系数，k是多孔电极的渗透率。

3.3 能量守恒方程

管式SOFC的能力守恒方程可由下式表示[8-9]

∇·(ερCpuT)=∇·(keff∇T)+ST

(10)

式中，ST为热源项，keff为有效导热系数(W/m*K)。

3.4 边界条件

在流道入口处初始流速和温度参数采用定值，即

uin=u0；Tin=T0

在流道出口处假设气体流速和温度已达到稳定状态，即

4 尺度耦合的技术路线

对管式SOFC采用多尺度耦合模拟的技术路线[8]如图2。在宏观尺度上，采用三大守恒定律建立CFD模型模拟流动和传热问题。在介观尺度上，采用格子玻尔兹曼方法模拟反应扩散问题，其原因是SOFC阳极为多孔介质，分布不连续，采用宏观模型难以进行准确描述。

图2 管式SOFC多尺度耦合模拟的技术路线图

本文的多尺度耦合计算流程如图3。假设任意初始温度T，根据式(1)～式(4)计算反应放热量，并以此作为边界条件，输入流动和传热的CFD模型中计算温度场，取其平均值Tc， 将T和Tc进行比较，当相对误差关系式err=(T-Tc)/Tc成立，则退出计算，计算流程如图3。

图3 计算流程图

5 数值模拟结果及分析

5.1 格子玻尔兹曼电化学反应模拟结果及分析

本文采用格子玻尔兹曼D2Q9模型模拟管式SOFC的电化学反应扩散现象，二维网格划分为256×256，H2、O2的初始浓度分别为15.2250mol和9.9167mol，网格初始分布为正态随机分布。管式SOFC的电流密度可由下式表示[10]

(11)

式中，z为反应的H2的量。

图4所示为1000℃下反应达到平衡后H2的浓度分布。图中x轴表示管式SOFC阳极的厚度，y轴表示管式SOFC阳极的长度，z轴表示反应达到平衡后H2的浓度分布。图4结果表明反应达到平衡后阳极内H2浓度呈多个正态随机分布，这种正态分布与所采取的初始分布相关，也与反应进行的程度相关，反应进行的越彻底，各网格节点上浓度的偏差越小，因此随着反应温度的升高，反应进行的程度越高，阳极内各节点达到平衡后，其浓度分布越均匀。

图4 1000℃时反应平衡后H2浓度分布

图6 电流密度与反应温度关系

图7 反应时间与反应温度关系

图5、6、7所示为H2的浓度、电流密度、反应时间与反应温度的关系。结果表明，H2的浓度随着反应温度的增加而逐渐减少；电流密度随着反应温度的增加而逐渐变大，电流密度的大小和H2量有关，H2量越多电流密度就越大；反应达到平衡所需要的时间随着反应温度的增加变得越来越短；本文在模拟时，模拟所需要的相关参数，如指前因子、活化能等参数都是依据有关手册选取的，这可能与实际情况存在偏差，导致800℃时反应所需时间较长。

5.2 耦合数值模拟结果及分析

图8为空气入口流速为1.634m/s时管式SOFC多尺度耦合模拟的温度场分布。结果表明，空气温度和管壁温度在入口处变化较大，随着空气的流动，温度很快达到稳定。管式SOFC内部的温度可以通过控制空气的流动状态来控制，因为空气起着冷却的作用，随着空气流速的增加，被空气带走的电化学反应的放热量越多，温度降低的越多，达到平衡后便处于稳定状态，且温度分布趋于均匀。

图8 管式SOFC多尺度耦合模拟的温度场分布

在流动方向上，随着电化学反应的不断进行，放出的热量不断增多，流道中的气体被不断加热，气体温度逐渐升高，这加速了电化学反应的进行，使得放热量快速增加。另外，在这过程中管壁也被不断加热，温度逐渐升高。管式SOFC整体温度分布情况为入口温度最低，出口温度达到最高值。因此，管式SOFC内部的温度分布与反应放热量相关。

综上所述，可以得出影响管式SOFC内部温度分布的两个重要因素，一个为反应放热量，另一个为空气的流动状态。

通过多尺度耦合模拟，发现随着燃气的流动，反应发生的位置不断向后移动，反应产生的热量不断增加，对流道和管壁加热效果越明显。因此沿着管式SOFC流动方向，流道和管壁的温度越来越高。

图9所示为计算稳定后管式SOFC的电解质壁面的温度分布，图中以1.634m/s和1.0m/s为例，分别以热量和对流换热作为边界条件计算温度分布。结果表明，采用热量作为边界条件的温度整体上比采用对流换热作为边界条件的温度高。当采用热量作为边界条件时，流速高的温度高，其原因在于热量边界对温度分布影响大，流速不足以将热量全部带走。当采用对流换热作为边界条件时，流速低的温度高，其原因在于对流换热边界中对流效应影响大，空气流速对温度分布起主要作用，流速越高，温度越低，因此稳定后流速低的温度高。

图9 管式SOFC电解质壁面温度分布

图10、11所示为空气流道入口流速u=1.634m/s时的速度场分布和管式SOFC中截面的速度分布。结果表明，管式SOFC的中部流速最大，壁面处流速最小。

图10 管式SOFC流场分布

图11 管式SOFC中截面速度分布

图12是管式SOFC管壁中轴线上的速度分布，结果表明多孔介质中的流速分布与多孔介质的准周期性分布一致。

图12 管式SOFC管壁中轴线速度分布

6 结论

本文从介观层次和宏观层次出发，介观层次上采用格子玻尔兹曼方法建立了管式SOFC的电化学反应扩散模型，宏观层次上采用计算流体力学方法建立了管式SOFC的CFD模型，并将格子玻尔兹曼方法计算的反应放热量作为边界条件求解CFD模型，实现了多尺度耦合模拟，得到了管式SOFC的稳态数值模拟结果。从结果中可知，H2的浓度和反应时间随反应温度逐渐降低，电流密度随反应温度逐渐增加。空气流道中的气体温度和管壁温度沿着流动方向逐渐升高。影响温度分布的两个重要因素，一个是反应放热量，另一个是空气的流动状态。反应放热量越多，管式SOFC内部的总体温度越高。另外，空气流速越高，管式SOFC内部温度降低的越多，达到平衡后便处于稳定状态，且温度分布趋于均匀。

结果表明采用多尺度耦合方法模拟管式SOFC是行之有效的，在一定程度上解决了管式SOFC介观尺度和宏观尺度多物理场耦合的问题。