金属支撑型固体氧化物电解池的3D 建模与性能分析

张梦茹，王恩华，胡浩然，欧阳明高，王贺武，卢兰光

（1.北京理工大学 机械与车辆学院，北京 100081；2.清华大学 车辆与运载学院，北京 100084）

为促进我国的碳达峰和碳中和工作顺利进行，加快产业结构的优化[1]，需要积极发展清洁能源[2-3].氢能作为21 世纪的“终极能源”，是一种高效、零碳、应用场景丰富的环境友好型绿色能源[4-5].固体氧化物电解池（solid oxide electrolysis cell, SOEC）采用全固态电池结构[6]，具有效率高[7-8]、成本低、寿命长等诸多优势[9-10].相比低温的碱性电解池和质子交换膜电解池，SOEC 在高温下运行，能有效降低电能消耗，提高电解效率[11].

与低温燃料电池技术相比，SOFC 发电具有很大的发展潜力，但技术尚未成熟.由于研究SOFC 内部的详细过程的实验成本高且难度大，因此仿真建模是SOFC 系统设计过程中非常重要的工具[12].诸多学者开展了SOFC 的流体仿真工作[13-15]，电池的三维结构设计对电堆内的流场和浓度场等的分布有重要影响[16-18]，进而影响到电堆的工作性能.CHEN 等[19]建立了单电池SOFC 堆的三维电化学—多物理场耦合模型，研究了平行圆柱流场对SOFC 性能的影响.采用新型平行圆柱形流场的SOFC 电堆的电流密度可提高38%，反应物分布均匀性提高了60%.GUO等[20]研究了连接体设计对SOFC 性能的影响，采用三维CFD 对比了传统肋、离散肋、交错菱形肋、窄肋、连续波状肋、梯形肋和交错长方体肋的工作性能.与传统的基于平行流道的连接体设计相比，新型的连接体可以实现SOFC 阴极中的O2更均匀分布，SOFC 的峰值功率密度可调高27.86%.

金属支撑SOEC 可在700 °C 以下的中低温环境下工作，有效提高电池的可靠性和耐久性[21].目前关于金属支撑型的数值研究较少，并且基本上采用零维模型进行分析，对电堆内部工作的电流密度和温度场等的分布特性还缺乏了解.ZHANG 等[22]建立了金属支撑型SOEC 的零维数学模型，详细分析了结构参数和运行条件对电池性能的影响.结果表明，电极的孔隙率和弯曲度影响到电极的质量输运，从而影响到浓差损失和电池性能.此外，运行电压越高越有助于电池电解，但是当电压超过0.3 MPa 时，其对电解性能影响的显著性降低.首先，零维模型仅能在有限的变化因素上解释电池性能降低的原因，但是不能直观描述产生电池性能降低的变化关系；其次，低维的模拟需要附加更多的假设条件，因此与实际电池中发生的过程存在较大的出入.因此，针对设计的一种铈基电解质的金属支撑SOEC，本研究采用三维数值仿真方法分析了SOEC 单电池的电化学工作性能和内部多物理场分布特性.首先根据设计的金属支撑SOEC 单电池，建立了单通道的三维多物理场模型.设计的金属支撑SOEC 在650 °C 的高温下运行，内部的热辐射可能会影响SOEC 的工作性能.因此，建立的模型还耦合了辐射传热模型.随后，基于仿真结果分析了该金属支撑SOEC 在650 °C 的电解性能，并进一步详细地分析了单电池内部的物质传输和电化学反应以及温度场的变化特性.研究结果可为下一步金属支撑SOEC 的设计开发提供参考.

1 金属支撑型SOEC 物理模型

GDC 电解质基[23]的金属支撑固体氧化物燃料电池可以在中低温下展现较好的电化学性能.由于SOEC 的工作电流较大，ZHANG 等[22]研究了一种可降低电解质氧分压的金属支撑SOEC，其单电池结构如图1 所示.

图2 表示SOEC 工作时，内部发生传热传质和电化学反应等物理化学过程[24].

图2 SOEC 工作原理图Fig.2 SOEC working principal diagram

用不可压缩的纳维-斯托克斯(N-S, Navier-Stokes)方程描述开放流道中的流动，Brinkman 方程描述多孔电极内的流动.采用ButlerVolmer 方程描述电极反应层内电荷转移的电流密度[25]，氢气侧：

氧气侧有：

基于局部热平衡的假设，介质内的能量方程可表示为

式中：λ为热导率，W/(m·K)；cp为比热容，J/(kg·K).

对于空气和燃料流道壁面，包含热辐射在内的净辐射热量密度为[28]

式中：h为对流传热系数，W/(m2·K)；qr为表面净传热速率，W/m2.

根据图1 设计的金属支撑SOEC 单电池结构，建立了单通道的单电池几何模型，具体如图3 所示.表1列出了具体的SOEC 单电池的几何结构尺寸.

表1 金属支撑SOEC 单电池几何结构尺寸Tab.1 Geometric parameters of designed metal-supported SOEC

图3 建立的单通道金属支撑SOEC 几何模型Fig.3 Geometric model of single-channel metal supported SOEC

建立的金属支撑SOEC 各部分材料的热物理和相应的电池的主要参数可参考文献[17, 29].

建立的多物理场的边界条件设置如表2 所示.

表2 金属支撑SOEC 多物理场模型边界条件Tab.2 Boundary conditions of multi-physical field model for metal supported SOEC

为研究金属支撑SOEC 的工作性能，首先对SOEC模型的计算精度进行了验证.基于文献[17]中的SOFC模型，在COMSOL 6.0 软件中建立了相应的三维仿真模型，仿真得到的结果与文献的对比如图4 所示.仿真结果与实验结果的误差保持在±3.05%以内，具有良好的一致性.

图4 建立的三维多物理场模型预测结果与文献[17,25]的结果对比Fig.4 Comparison of predicted results of the 3D multi-physical field model with those in literature [17,25]

2 结果与讨论

设SOEC 的初始工作温度为923.15 K，运行压力为0.1 MPa, 氢气侧水蒸气流量为3.0 m/s，氧气侧空气流量为5.0 m/s.利用COMSOL 6.0 软件，对设计的金属支撑SOEC 内部的工作过程进行了仿真分析.

2.1 极化曲线

图5(a)为金属支撑SOEC 的极化曲线.随着电流密度的增加，电池工作电压逐渐增大.开路电压为1.05 V，当电流密度为3A/cm2时，输出电压为1.46 V，增加了0.41 V，说明金属支撑SOEC 在大电流密度下运行时存在较大的极化损失.基于COMSOL 软件的计算结果，利用MATLAB 计算了MS-SOEC 的各极化损失的具体数值，计算结果如图5(b)所示.图5(b)显示了金属支撑SOEC 工作时欧姆损失和极化损失随电流密度的变化曲线.欧姆损失随电流密度的增加近似线性增加.当电流密度为0.4 A/cm2时，总电压损失为0.19 V，极化损失和欧姆损失分别占84.26%和15.74%；电流密度为2.4 A/cm2时，总电压损失为0.38 V，极化损失和欧姆损失分别占66.28%和33.72%.不同电流密度工况下，极化损失的占比始终明显大于欧姆损失，这是由于小电流密度下活化极化损失占比较大，而在大电流密度下浓差极化损失占比较大，因此作为二者之和的极化损失占比始终较大.

图5 650℃下设计的金属支撑SOEC 的工作性能Fig.5 Work performance of metal supported SOEC under 650 ℃:

2.2 流场分析

流道内H2O、H2、O2的浓度沿x轴方向的变化曲线如图6 所示，其中的浓度值为在对应x处的流道yz平面内的平均值.沿着流动方向，氢气侧的H2O 的摩尔浓度逐渐下降，对应的H2浓度从入口处的0.1逐渐增加到出口处的0.31，氧气侧的O2浓度从入口的0.21 逐渐增加到0.26.图中还显示了氢气侧和氧气侧的平均温度变化.氢气侧的温度变化范围较大，温度先从949.36 K 快速上升，到1.16 cm 附近达到最大值954.55 K，随后又稍有下降，出口处温度为954.03 K.氧气侧的温度变化呈现相似的趋势，但是变化范围明显小很多，基本维持在931.0 K 附近，这主要是因为氧气侧通入的空气流量较大，与阴极之间的热交换量也更大，避免了温度的快速上升.前半段的温度迅速上升主要受到水蒸气浓度的影响，靠近入口附近，水蒸气浓度较大，电化学反应的活化极化较低，随着H2O、H2和O2浓度的变化，过电势加大，导致产热增加，温度升高.

图6 金属支撑SOEC 内沿流道方向主要组分的平均摩尔浓度及温度变化曲线Fig.6 Mean concentrations of main components and temperature curves along the flow passage in the metal supported SOEC

金属支撑体和流道内的H2浓度沿X轴的分布云图如图7 所示.图7(a)为金属支撑体中间位置截面的H2浓度分布，沿X轴方向的变化趋势与图6 中的结果类似.同时，在接近进出口位置处，连接体接触界面附近的H2浓度与流道附近的H2浓度差别较小，但是在中间位置附近二者之间的差异较大.图7(b)为对应流道与金属支撑体接触界面处的H2浓度分布，其变化趋势与多孔金属支撑体内的变化趋势一致，仅在数值上有0.02 左右的下降.

图7 金属支撑体中间截面的和金属支撑与氢电极接触界面的H2 摩尔浓度分布Fig.7 Distribution of H2 concentration in the middle plane of metal support and the interface between metal support and hydrogen electrode

2.3 电流密度分布

图8 显示了GDC 主电解质内部的氧离子电流密度，在左侧入口处的电流密度为2.67 A/cm2，电流密度沿X轴方向逐渐降低，右侧出口处的电流密度为2.18 A/cm2.这主要是由于参与反应的H2O 的浓度沿轴向方向逐渐减小.在电池运行过程中，在靠近入口附近，电流密度变化较为剧烈，在连接体与单电池接触区域的边缘的电流密度明显高于附近区域.另一方面，在x=0.4 cm 处，对应流道中心点区域出现电流密度突然增大的现象，这主要是因为此处的H2O 浓度较高，同时局部温度也较高，导致此时的电化学反应较为剧烈.

图8 沿厚度方向GDC 电解质中心位置处的xy 截面电流密度分布Fig.8 The xy section current density distribution at the center position along the GDC electrolyte thickness direction

2.4 温度场分析

图9 显示了在x轴向x=0.1 cm 和1.9 cm 的两个不同切面上的温度分布云图.图9(a)中最高温度位于GDC 电解质的中心点，约为974.38 K，沿着GDC电解质中心向两侧逐渐降低，左右两侧的温度降低到947.17 K.在10Sc1CeSZ 电解质层，也呈现类似的变化规律，温度从中心的962.03 K 逐渐降低到两侧的947.16 K.相应地，金属支撑的温度从中心区域的958.67 K 向两侧逐渐降低，边沿处的最低温度为950.28 K.连接体内部的温度分布总体上比较均匀，氢气侧的连接体温度约为953.87 K，氧气侧的连接体温度约为926.19 K.图9(b)显示了类似的变化趋势，但是中心区域的高温区的温度有所升高，且高温区域的面积明显扩大，在氢气侧和氧气侧的流道内也分别出现了局部的高温区.在x=0.1 cm 截面内的GDC电解质的最高温度为974.38 K，在x=1.9 cm 截面内的GDC 电解质的最高温度为976.39 K.在yz截面内的温度分布主要由三个方面的原因：首先，最高温度出现在中心位置，也即是电极和电解质的反应区域，虽然SOEC 是吸热反应，但是由于电化学反应剧烈，电流密度较大，产生了大量的欧姆热；其次，支撑体的导热率为20 W/(m·K)，使得从电池中心位置产生的热量在金属基底和连接体内快速传播；最后，造成进出口不同位置处二者分布的原因，一是氧气侧温度分布相似，氧气侧通入的空气流量较大，与阴极之间的热交换量也更大，避免了温度的快速上升；二是由于氢气侧是相对流量较小的水蒸气，因此对电化学反应热量的反应较为敏感，温度沿轴向方向变化幅度较大.

图9 金属支撑SOEC 在进出口附近的温度分布云图Fig.9 Temperature distribution on two typical positions of YZ section in the metal supported SOEC

3 结 论

为研究金属支撑型固体氧化物电解池的传质传热性能，建立了计算金属支撑型固体氧化物电解池性能的三维多物理仿真模型.详细研究了电解池内部的速度、浓度和温度分布以及电化学反应情况, 对多层电解质设计的金属支撑型固体氧化物电解池的性能进行了评估.主要结论如下.

①沿着平行流道方向，H2的摩尔浓度从0.1 近似线性增加到0.31，O2的摩尔浓度从0.21 近似线性增加到0.26，而H2O 的浓度从0.9 逐渐下降0.64，H2摩尔浓度上升速度大于O2摩尔浓度.在流道的前半段，氢气侧和氧气侧流道内的流体温度出现较为明显的升高，氢气侧的流体温度从949.36 K 上升到954.55 K，氧气侧的流体温度从930.57 K 上升到932.07 K，在流道的后半段，二者的温度基本维持不变，仅有轻微的下降.

②采用金属支撑有利于提高SOEC 内部温度均匀性.在流道入口的x=0.1 cm 截面内的GDC 电解质的平均温度为949.88 K，在流道出口附近的x=1.9 cm截面内的GDC 电解质的平均温度为954.97 K，温度仅增加了5.09 K.