建筑模型
以实际生产的SOFC为原型，研究表明，SOFC周期单元由阴极支持层（CSL）、阴极活化层（CAL）、电解质层（EL）、阳极活化层（AAL）、阳极支持层（ASL）、连接体（）等构成，PEN组件由阳极材料（-onia、NiO-YSZ）构成。
电解质材料为（onia，YSZ），阴极材料由（La0.6Sr0.4Co0.2FeO3-δ，LSCF-GDC）组成，阴极和阳极为双相复合材料，电子导体的体积分数分别为0.4和0.5。图1为SOFC几何模型及PEN组件内物质传递示意图，表1为模型的全局参数。
表1 模型全局参数
当使用碳氢化合物（CH4、H2、H2O）作为燃料时，MSR反应生成的CO也参与电化学反应，阳极侧的电化学反应可利用公式（1）和公式（2）计算。
（1）
（2）
阴极侧的电化学反应可采用公式（3）计算。
（3）
同时镍基阳极中会发生MSR和WGSR反应，当反应温度高于700 ℃时，一氧化碳歧化反应和一氧化碳还原反应均受到抑制，只有甲烷裂解反应（MC）随反应温度的升高而增大，因此本模型仅考虑甲烷裂解反应产生积碳，具体反应方程如（4）~（6）所示。
（4）
（5）
（6）
许多研究者利用碳沉积活性来预测SOFC内部热力学碳沉积的可能性。魏等分析了MSR的基元反应，认为重整反应过程中总的碳沉积活性取决于镍基催化剂表面化学吸附的碳（C*）浓度。碳沉积活性αcar采用公式（7）计算，当αcar大于1时，从热力学上表明会发生碳沉积。
（7）
基于反应动力学阿伦尼乌斯方程，MSR、WGSR和MC的反应速率控制方程如表2所示，SOFC内部物理场的传递过程可以用相应的控制方程和源项表达式来表示。
表2 反应速率控制方程
在应力模型中假设电池各部件均发生弹性变形，力学理论符合胡克定律，控制方程如表3所示，不同温度条件下的热力学参数如表4所示。
3 应力模型控制方程
在SOFC从烧结结束冷却到室温，再由室温加热到高温工作状态的过程中，假定NiO完全还原为Ni，忽略Ni-NiO反复氧化还原引起的体积收缩。
表4 不同温度条件下的热力学参数
基于Chin模型计算正极和负极双相复合材料的物理性能，取烧结温度（1400℃）作为无应力参考温度Tref，选取20℃时E、v、αth作为对应值，σ0＝0，ξ0＝0，计算制造残余应力；取残余应力和残余应变作为预应力σ0和预应变ξ0，取室温（20℃）作为无应力参考温度Tref，选取800℃时E、v、αth作为对应值，计算高温条件下的综合残余应力以及热力学参数变化对热应力分布的影响。
燃料流速为15.4 mL/min,空气流速为61.7 mL/min,入口温度为800 ℃。阳极气道中碳氢燃料H2O/CH4/H2的摩尔分数比为0.675/0.225/0.1,水碳摩尔比为3;阴极气道中O2/N2的摩尔分数比为0.233/0.767。连接器下表面设置为地,上表面设置为工作电位;采用“刚体运动抑制”条件模拟SOFC的夹持力,其他边界条件设置为自由边界。
模型验证
该模型已被验证为网格独立模型，对于纯氢作为燃料，生成了三个高质量的结构化网格，网格数分别为15440、23100和27400，如图2所示。
经计算发现前两个网格的数值解曲线偏差较大，进一步细化网格后，两个网格的数值解曲线基本重合，在允许的误差范围内，认为网格数23 100满足要求。
为验证模型的准确性，操作条件、材料性质及初始条件设置为与参考实验相同的操作条件。在阳极气道中，当使用氢气作为燃料时，H2/H2O摩尔分数比为0.95/0.05；当使用碳氢化合物作为燃料时，H2O/CH4/H2摩尔分数比为0.675/0.225/0.1，水碳摩尔比为3；阴极气道中O2/N2摩尔分数比为0.233/0.767（见图3）。当使用氢气作为燃料时，最大相对误差为6.6%；当使用合成气作为燃料时，最大相对误差为9.1%。
为了验证模型应力性能的准确性，操作条件、材料性质、初始条件均设置为与参考实验相同的操作条件，计算得到的最大压应力值与范等的实验结果和数值模拟结果进行比较，如表5所示，误差均在8%以内。
表5 YSZ层最大压应力对比
结果分析
基于模拟结果，分析了典型稳态条件下SOFC流场设计（并流和逆流）对电池内部各物理量的影响。
图4、图5分别给出了正向和反向流模式下阳极活化层中段中心线处沿y轴正方向的温度分布。从图中可以看出，SOFC内温度场分布受工作电压和进气方式的影响。工作电压越低（电流密度越大），温度越高；相同电压下，反向流模式下的最大温差小于正向流模式下的最大温差，即反向流模式下的温度分布更均匀。
当工作电压小于0.7 V时，两种进气模式的温度均大于1073.2 K（进气温度），表明电化学反应和极化损失释放的总热量大于MSR吸收的热量。
此外，在此条件下，H2的摩尔分数沿y轴正方向逐渐减小，表明电化学反应消耗H2的速率大于MSR和WGSR生成H2的速率；CO的摩尔分数沿y轴正方向先增加后减小，表明在电池前期，MSR生成CO的速率大于电化学反应和WGSR消耗CO的速率，随后电化学反应和WGSR消耗CO的速率大于MSR生成CO的速率。
相反，当工作电压大于0.7V时，意味着MSR吸收的热量大于电化学反应和极化损耗所释放的总热量。
图6为工作电压0.6V、正向和反向流条件下电池表面温度分布。从图6可以发现，正向流时，单电池表面温度沿y轴正方向逐渐上升，最高温为1093K；反向流时，正负极气道入口处温度相对较低，最高温出现在电池中部，最高温为1081.6K。
高温主要是由于放热的电化学反应引起的，而低温则是多孔阳极内发生的强烈吸热的MSR反应和流道内气体的对流冷却共同作用的结果。
顺流模式下，沿y轴正方向，气道内气体的对流冷却及吸热反应MSR使电池前期温度上升减缓，随后温度上升加速了电荷交换过程，促使电化学反应释放更多的热量，使y=40mm处温度达到最高。
逆流模式下，由于空气与燃料气体流动方向相反，燃料入口和空气入口处的电化学反应较弱，主要受到气道内气体冷却对流的影响，温度较低；而在电池中部，电化学反应加强，释放出更多的热量，温度达到最高。
图7、图8分别给出了工作电压0.6V下正向和反向进气模式下CO2和CH4组分分布情况，从图中可以看出两种进气模式下气体组分分布趋势基本相同，由于CO2是电化学反应和WGSR反应的双重产物，因此CO2的摩尔分数分布沿y轴正方向逐渐增大；CH4的摩尔分数沿y轴正方向逐渐减小。
这主要是因为初始H2含量较低，有更多的CH4参与重整反应，补充电化学反应所消耗的氢气；在xz截面（y=20mm）处，顺流模式下CH4摩尔分数消耗大于逆流模式下的，且该截面处的浓度梯度也大于逆流模式下的；翅片下方多孔电极中CH4摩尔分数明显低于流道正上方对应电极中CH4摩尔分数，主要是因为沿y轴正方向不断进行内部重整反应，CH4可以直接从流道扩散到流道正上方对应的多孔电极中，而翅片下方的燃料气体需要较长的扩散路径，气体补充缓慢。
图9和图10分别给出了工作电压0.6 V、正向和反向流动模式下阳极支撑层中部横截面的MSR和WGSR反应速率分布。
从图中可以发现，两种进气模式下MSR反应速率沿y轴正方向逐渐减小，这是因为两种进气模式下y=0mm处CH4浓度最高，因此随着MSR反应的继续进行，反应速率沿y轴正方向逐渐减小；两种进气模式下WGSR反应速率沿y轴正方向先减小后增大。
这是因为在y=0 mm附近，MSR反应强烈，生成的CO浓度较高，因此该区域的WGSR反应速率沿y轴正方向逐渐减小，而随后WGSR与电化学反应CO-O2继续消耗大量的CO，导致该位置的WGSR进行逆反应试图提高CO含量，即WGSR的反应速率变为负值；翅片下方的多孔电极的反应速率小于对应的流道正上方的多孔电极，且顺流模式下的反应速率比逆流模式下的分布更均匀。
图11为工作电压0.6 V、正向和反向进风模式下阳极支撑层中段MC反应速率分布。从图中可以看出，MC反应受到抑制，最大反应速率仅为9.97×10-4mol/(m3·s)；两种进风模式下MCR反应速率沿y轴正方向逐渐减小，正向进风模式下MCR反应速率比反向进风模式下分布更均匀。
图12为工作电压为0.6 V时不同水碳比在阳极支撑层中段的积碳活性分布。从图中可以看出，两种进气模式下积碳活性都在y=0 mm（燃油进气端）附近达到最大值，随后沿y轴正方向逐渐减小。
这是因为y=0 mm处CH4浓度最高，有利于甲烷裂解反应向正方向发展，但总体来看，两种进气方式的碳活性都远小于1。另外可以发现，调整水碳比之后，高积碳活性区域由燃料入口处扩展到SOFC中部，最大积碳活性也随之增大。这说明合理提高水碳比可以有效防止热力学积碳，但水碳比不宜过高，因为高比例的水蒸气会造成Ni基失活，影响电池性能，降低系统效率。
第一主应力适用于陶瓷、金属体系的分析，第一主应力可以区分最大拉应力和压应力（忽略垂直于平面的剪应力），因此本文主要比较第一主应力，正负值分别代表拉应力和压应力。图13给出了工作电压0.6V及正向和反向流动模式下第一主应力的分布情况。
从图中可以发现，顺流模式下，SOFC整体第一主应力变化范围为-24.4~299.2MPa，大于逆流模式下第一主应力变化范围-23.9~296.9MPa；多孔阳极支撑层在xz截面（y=20mm）处拉应力最小，连接件对应的阳极支撑层由于受到下连接件的约束，受到压应力；阳极支撑层厚度方向存在应力梯度，且由于几何形状和材料的不连续性，在相邻部件接触的角部区域出现应力极值。
图14为顺流、逆流模式下各部件第一主应力最大值分布情况，从图中可以发现，两种进气模式下，阴极支撑层拉应力最大，阳极活化层拉应力最小；从电解质层到阴极活化层存在应力突变，其值由32MPa突然增大至。
首先，电解质在所有组分中热膨胀系数最小，杨氏模量最大，对外界影响较为敏感；其次，电解质和正极材料性能不连续，正极较薄，热膨胀系数较大，导致正极侧拉应力最大，从而导致电解质和正极活性层发生突变。
综上所述
通过三维多物理场全耦合模拟，讨论和分析了顺流、逆流两种进气方式下SOFC内部温度、气体成分、化学反应速率、积碳活性及热应力的分布，得到以下结论。
SOFC内部温度场分布受工作电压和气流结构影响，工作电压越低，温度越高；相同电压下，逆流电池内部最大温差小于顺流电池内部最大温差。
气流配置通过影响同一位置气体组分的分布来改变SOFC内部的反应速率。典型条件下，在xz截面（y=20 mm）处，顺流模式的CH4摩尔分数消耗大于逆流模式，且该截面处的浓度梯度也大于逆流模式；顺流模式下反应速率分布更加均匀；甲烷裂解反应受到抑制，积碳活性远小于1。
气流配置对SOFC热应力分布影响不大。典型工况下，顺流模式热应力略大于逆流模式；阴极支撑层拉应力最大，阳极活化层拉应力最小；从电解质层到阴极活化层存在应力突变，其值由32MPa突然增大到40MPa。此外，极端应力常常出现在与其他部件接触的拐角区域。