质子交换膜燃料电池催化层局部氧气传输过程仿真方法

技术领域

本发明属于电化学燃料电池领域，具体涉及一种对质子交换膜燃料电池催化层氧气传输过程建立计算模型的方法。

背景技术

质子交换膜燃料电池可将氢能转换成电能而不产生污染物，并具有功率密度高、转换效率高和快速启动等优势，被广泛认为是一项极具前景的清洁能源技术。近年来，质子交换膜燃料电池的功率密度得到了显著提升，推动了其快速发展，然而在其大规模商业化应用前，仍有一些技术需要突破。质子交换膜燃料电池的催化层是其电化学反应发生的场所，催化层是由催化剂聚团和全氟磺酸聚合物电解质薄膜构成的多孔结构，其内部孔隙中的氧气分子需要穿越覆盖在催化剂上的电解质薄膜，才能到达催化剂表面发生反应，由此产生了很高的氧气传输阻力，严重限制了高性能低铂燃料电池技术的发展。

质子交换膜燃料电池催化层中催化剂颗粒的尺寸仅有数纳米，很难利用实验的方法对其周围的传输现象进行观测。分子动力学模拟则是探究微纳尺度传输现象最合适的方法之一，能够帮助探究氧气穿越电解质薄膜到达催化剂颗粒表面的现象，指导质子交换膜燃料电池催化层电极设计，从而可大幅降低实验成本并缩短研发周期。

催化层结构复杂尤其是催化剂载体(无定形碳)属于无序结构，传统分子动力学建模研究可对其结构进行简化，比如将催化剂载体简化为有序的石墨烯层，石墨烯层内部碳原子排布与无定形碳不同，同时石墨烯层是平面结构，无定形碳为颗粒结构，该简化将影响局部结构和氧气传输过程，导致物质密度分布和氧气传输流量预测不准。为此，本发明提供了一种考虑催化层中催化剂、催化剂载体和电解质真实结构的分子动力学模拟方法。特别地，该模型详细考虑了作为催化剂载体的无定形碳真实结构以及多个催化剂颗粒在载体上的分布。所以该模型的建立能够对催化剂周围的局部结构和氧气传输过程做出良好预测。

发明内容

本发明的目的是，提出一种基于分子动力学模拟为基础的质子交换膜燃料电池催化层局部氧气传输过程仿真方法。模型涉及到催化剂颗粒、催化剂载体和电解质薄膜的真实结构，以及电解质水合、水分子运动、质子运动和氧气跨膜传输等过程。特别是考虑了催化剂载体的无定形结构以及多个催化剂颗粒在载体上的分布，且模型能够根据实际情况调整相关参数，从而展示质子交换膜燃料电池催化层催化剂周围的氧气传输过程。

基于分子动力学模拟的燃料电池催化层氧气传输过程仿真方法的建立，包括了6个部分：构建催化剂载体模型；构建催化剂模型；构建全氟磺酸聚合物、水分子、水合氢离子和氧气分子模型；构建电解质薄膜覆盖催化剂及其载体的初始模型；构建电解质薄膜覆盖催化剂及其载体的平衡模型；以及构建氧气穿越电解质薄膜传输过程模型。其具体步骤如下：

(1)构建催化剂载体模型

以无定形碳作为催化剂载体，采用混合逆向蒙特卡洛方法(Hybrid reverseMonte Carlo，HRMC)拟合实验衍射结构信息并最小化系统能量来生成无定形碳结构。其构建过程为：首先，在设定体积的模拟盒子中随机分布碳原子，根据无定形碳的密度来确定碳原子的数目。然后，随机选择一个碳原子,把这个原子的位置标为Old；之后，模拟尝试将这个碳原子移动到一个新的随机位置(与原来位置之间的距离不能超过设定的最大步长距离)，将其标记为New。随后，对本次移动进行判定，其依据是使误差函数χ不断最小化，该误差函数包括了对整体结构和能量的约束，旧位置Old和新位置New的误差函数分别为：

S为结构因子函数，g为径向分布函数，B为键角分布函数。S、g和B函数均需要离散取值，将函数在特定区间内分别线性离散成i、j和k个数据点，分别计算S(i)、g(j)和B(k)；下标old表示原子移动前(旧位置Old)的数值，下标new表示原子移动后(新位置New)的数值，下标exp表示实验数值；W

整个模拟过程从高温开始，模拟线性退火过程，温度逐渐降低到室温，从而得到稳定的无定形碳结构。计算得到原子移动前和移动后的误差函数，如果χ

(2)构建催化剂模型

以铂作为催化剂，依据铂的面心立方晶体结构，其晶格常数为0.39239nm。在给定的空间中填充铂原子，形成不同大小、不同形态的催化剂颗粒。

(3)构建全氟磺酸聚合物、水分子、水合氢离子和氧气分子模型

以全氟磺酸聚合物为电解质材料，根据全氟磺酸聚合物、水分子、水合氢离子和氧气分子的分子式，确实原子之间的成键关系，构建其分子模型。

(4)构建电解质薄膜覆盖催化剂及其载体的初始模型

以催化剂载体(即无定形碳)为基础，在其表面随机放置若干个催化剂颗粒(铂颗粒)，催化剂颗粒的数量通过实际的铂/碳质量比确定。同时在催化剂上方随机插入若干个全氟磺酸聚合物分子、水合氢离子和水分子，其中全氟磺酸聚合物分子通过实际的电解质/碳质量比确定，水合氢离子数量通过模拟系统中负电荷数量确定，从而保证模拟系统的电中和。水分子数量通过实际电解质的水含量确定，电解质水含量表示水分子和水合氢离子数量之和与全氟磺酸聚合物中磺酸集团数量的比值。

(5)构建电解质薄膜覆盖催化剂及其载体的平衡模型

首先，对系统进行能量最小化处理。然后，在初始模型的最上方放置一个虚拟的势能墙，势能墙与原子之间的势能关系式为：

将催化剂载体和催化剂颗粒视为固定不动，沿着电解质厚度方向移动势能墙，压缩电解质直到其密度达到一定数值。然后，保持势能墙的位置不变，将全氟磺酸聚合物之间的势能以及催化剂/载体与其它分子之间的势能降低到原来的0.01。在恒定粒子数、恒定体积和恒定温度(以下简称NVT系统)下，维持系统温度为高温(如1000K)，持续设定时间。紧接着再进行设定时间的模拟，在此过程中，将势能逐渐恢复到初始水平，将温度逐渐降低到模拟温度。之后，在NVT系统下，使系统温度在高温和模拟温度之间来回变化多次。最后，将势能墙放置到不影响电解质形貌的位置。仍将催化剂载体视为固定不动，将催化剂颗粒看作刚体(即将颗粒看成一个整体进行移动或转动)，在NVT系统下继续模拟设定时间，直到系统的能量和Z方向的密度分布不再发生变化。

(6)构建氧气穿越电解质薄膜传输过程模型

在电解质薄膜上方插入设定数量的氧气分子，该数量可依据实际情况来确定,并将势能墙放置在氧气上方。在NVT系统下，维持系统中氧气、电解质和催化剂的温度为模拟温度，持续模拟设定时间，从而使系统达到氧气稳定传输的状态。再进行设定时间的模拟，从而对氧气传输过程进行分析。在整个模拟过程中，如果氧气分子中任一原子与催化剂颗粒任一原子之间的距离小于0.31nm，则认为氧气分子到达了反应位点，从而将其从系统中去除，并在系统的最上端随机插入氧气分子，从而维持系统中氧气分子数量不变。为了提升计算效率，可以每隔设定时间对系统中氧气分子的位置进行一次核查。

根据上述步骤，可建立基于分子动力学模拟的燃料电池催化层氧气传输过程仿真方法。经过数值仿真，即可获得氧气传输流量、氧气密度分布、全氟磺酸聚合物密度分布、水密度分布等。

本发明的特点以及有益之处在于：

(1)综合考虑质子交换膜燃料电池催化层中，电解质薄膜包裹催化剂的微观结构以及电解质薄膜水合、水运动、水合氢离子运动、以及氧气传输现象，建立了基于分子动力学模拟的催化层氧气传输过程仿真方法，包含对如催化剂周围电解质薄膜形貌和氧气传输路径等关键性问题的合理解释。提供了一个较完整的以实验方法难以观测的催化剂周围微观传输现象的方法，为进一步研究催化层局部氧气传输机理和局部结构优化设计提供了有效的工具。

(2)目前所用模型大都忽略了催化剂载体和催化剂颗粒的真实三维结构，本模型利用逆向蒙特卡洛方法实现对催化剂载体：无定形碳的重构，以及依据铂的面心立方晶格结构实现对催化剂颗粒的重构，从而得到了碳载铂的真实结构。

附图说明

图1模型搭建过程及无定形碳、催化剂颗粒、全氟磺酸聚合物、水分子、水合氢离子和氧气分子的结构示意图。

图2无定形碳模型与实验测得的结构因子S(Q)之间的比较。

图3全氟磺酸聚合物、水分子和水合氢离子的原子局部电荷。

图4不同模拟时间下穿越电解质薄膜的氧气分子数量以及氧气传输流量。

图5两个催化剂颗粒周围全氟磺酸聚合物的密度分布图。

图6两个催化剂颗粒周围水的密度分布图。

图7两个催化剂颗粒周围氧气的密度分布图。

图8 X-Z平面内的全氟磺酸聚合物密度分布云图。

图9 X-Z平面内的水密度分布云图。

图10 X-Z平面内的氧气密度分布云图。

表1结构因子函数实验数据。

表4径向分布函数实验数据。

表3键角分布函数实验数据。

表4分子动力学模拟的关联势能参数。

具体实施方式

以下通过实施例进一步说明本发明的设计方案，需要说明的是本实例是为了清楚解释建模步骤而进行的叙述性说明，并不以此限定本发明的保护范围。

构建催化层氧气传输过程仿真方法包括：构建催化剂载体模型；构建催化剂模型；构建电解质材料、水分子、水合氢离子和氧气分子模型；构建电解质薄膜覆盖催化剂及其载体的初始模型；构建电解质薄膜覆盖催化剂及其载体的平衡模型；以及构建氧气穿越电解质薄膜传输过程模型6个部分。

催化剂载体、催化剂、电解质材料、水分子、水合氢离子和氧气分子的结构以及模型搭建过程如图1所示。

计算步骤如下：

(1)构建催化剂载体模型

采用混合逆向蒙特卡洛方法生成无定形碳模型，碳原子之间的势能利用环境相关原子间势(environment-dependent interatomic potential)来表示，模拟盒子X方向长度为

结构因子函数S(q)在q∈[0,20]区间内被离散为j＝401数据点，径向分布函数g(r)在r∈[0,19.9878]区间内被离散为j＝1304数据点，键角分布函数B(θ)在θ∈[1,180]区间内被离散为j＝180数据点，其对应的实验值被列表1、表2和表3中。

模拟分为了四个阶段，初始温度设为6000K，然后温度逐渐降低至3000K、1500K和300K，每个降温过程运行1000万步，温度降低到300K后继续运行1000万步。

在第一阶段，即温度从6000K降低至3000K，成本函数S(q)、g(r)和B(θ)权重因子W分别设置为0.016、0.016和0.00008。

第二阶段，即温度从2000K降低至1500K，对应的权重因子分别设置为0.008，0.008，0.00004。

第三阶段，即温度从1500K降低至300K，对应的权重因子分别设置为0.004、0.004和0.00001。

第四阶段，温度维持在300K，对应的权重因子分别设置为0.001、0.001和0.000005。

经过以上四个阶段之后，在权重因子不变的条件下，再运行100万步进行数据收集，从而得到无定形碳的构型。为了节约计算时间，本实施算例选用半球形的无定形碳，即选中一个半径为4nm的半球形区域，然后删除多余的碳原子，从而形成一个半球形的碳颗粒，如图1所示(左上)。无定形碳模型结构因子与实验测得的结构因子之间的比较如图2所示。

(2)构建催化剂颗粒模型

依据铂的面心立方晶体结构，在半径为2nm的球形中填充铂原子，从而形成球形的铂颗粒，如图1(右上)所示。

(3)构建电解质材料、水分子、水合氢离子、氧气分子模型

全氟磺酸聚合物包括10个重复的单体，每个单体包括了16个CF

(4)构建电解质薄膜覆盖催化剂及其载体的初始模型

设置模拟盒子的X方向长度为9.2nm，Y方向的长度为9.2nm，Z方向的长度为50nm。将半球形无定形碳的球心位置设置为X＝4.1nm、Y＝4.1nm、Z＝0nm。在无定形碳表面随机放置2个催化剂颗粒，同时在催化剂颗粒上方随机插入30个全氟磺酸聚合物分子、300个水合氢离子和3000个水分子。

(5)构建电解质薄膜覆盖催化剂及其载体的平衡模型

首先，对系统进行能量最小化处理。然后，在初始模型的最上方放置一个虚拟的势能墙，势能墙与原子之间的关联参数ε

将催化剂载体和催化剂颗粒视为固定不动，沿着电解质厚度方向移动势能墙，压缩电解质直到其密度达到1.739g/cm

在NVT系统下，维持系统温度为1000K，持续模拟1ns；紧接着再模拟1ns，在此过程中，将势能逐渐恢复到初始水平，将温度逐渐降低到353K。

之后，在NVT系统下，使系统温度在1000K和353K之间来回变化8次，总共持续1ns。

最后，将势能墙放置到7nm处，仍将催化剂载体视为固定不动，将催化剂颗粒看成一个整体进行移动或转动的刚体，在NVT系统下继续模拟20ns，从而得到平衡结构，如图1所示(右下)。

(6)构建氧气穿越电解质薄膜传输过程模型

为了确保有足够多的氧气到达催化剂表面，从而增加样本，在平衡的电解质薄膜上方插入1800个氧气分子，并将势能墙放置到氧气上方，即Z＝13nm，如图1所示(左下)。在NVT系统下，维持系统中氧气、电解质和催化剂的温度为353K，持续10ns，从而使系统达到氧气稳定传输的状态，再进行15ns的模拟，从而对氧气传输过程进行分析。

在整个模拟过程中，如果氧气分子中任一原子与催化剂颗粒任一原子之间的距离小于0.31nm，则认为氧气分子到达了反应位点，从而将其从系统中去除。并在系统的最上端随机插入相应数量的氧气分子，从而维持系统中氧气分子数量不变。为了提升计算效率，可以每100fs对系统中氧气分子的位置进行一次核查。

在整个分子动力学模拟过程中，原子之间的势能方程为：

方程式右侧的五项分别表示键伸缩势能、键角弯曲势能、二面角扭曲势能、分子间作用势能、静电作用势能。

K

此外，第四项中不同类型原子之间的参数ε

原子之间的势能参数可见表4，极性分子全氟磺酸聚合物、水分子和水合氢离子的原子局部电荷可见图3。

整个仿真过程利用开源平台LAMMPS执行，时间步长为1fs，利用Nose-Hoover扩展系统法实现温度的控制，X和Y方向采用周期性边界条件。每2ps输出一次原子位置信息，从而计算平均密度分布；每3ns输出一次穿越电解质薄膜氧气分子平均数目，用于计算氧气传输流量。

附图4是穿越电解质薄膜氧气分子数量随时间的变化以及计算得到的氧气传输流量。

附图5-附图7是两个催化剂周围的全氟磺酸聚合物、水和氧气密度分布。从图5和图6中可以看出，催化剂附近存在全氟磺酸聚合物和水密度峰值，即存在致密超薄结构，该致密结构导致了很低的氧气密度，如图7所示。

附图8-10是X-Z平面内的全氟磺酸聚合物、水和氧气密度分布云图，从图8和图9可以看出，电解质薄膜中的水分子形成多个聚团结构，被全氟磺酸聚合物包围。

从附图10可以看出，氧气分子在电解质薄膜中的密度分布与电解质薄膜的结构密切相关，具有分散分布的特点，因此可以优化电解质薄膜的结构，尤其是全幅磺酸聚合物和水分子聚团的结构，来提升氧气密度和促进氧气传输过程。

表1结构因子函数实验数据

/>

/>

/>

表2.径向分布函数实验数据

/>

/>

/>

/>

/>

/>

/>

/>

/>

/>

/>

表3.键角分布函数实验数据

/>

表4分子动力学模拟的关联势能参数

/>

/>