A catalytic microkinetic analysis software

典型的INPUT文件见第5部分模拟示例,下面我们介绍输入参数具体含义,许多参数都有可靠的默认值,当你在输入文件中遗漏掉它们时,这些默认值就变得有用(这就允许用户可以有非常短的输入文件!)。那些没有默认值选择的参数就需要用户根据具体的体系进行设置,包括反应机理、反应能量和反应条件。文件中所有的内容均为应为字符,可以采用百分号(%)或感叹号(!)作为注释标记符号,分号(;)作为分隔符,除了基元反应,所有的参数设置均采用MATLAB赋值的方式进行,即格式为:变量 = 参数值(variable = value),也允许包含一些简单的MATLAB代码、函数和命令。

4.1 反应机理

微动力学的模拟,反应的基元步骤肯定是必不可少的,它定义了微动力学的反应模型,可以包括吸附、解离、结合和脱附反应,物种包含反应物、产物、表面中间体和反应位点。

反应物种的名称需要满足变量名的命名规则,最简单的可以采用化学式书写,比如水可以用H2O表示,推荐采用标准的化学式写法,有助于反应机理的可阅读性,以及程序自动识别匹配,实现自动分子量和自带的热力学数据的计算。分别采用(gas)和(liq)来区别气相和液相物种(或分别简化为(g)和(l)或(p)和(c)),例如H2O(gas)和H2O(liq)分别表示气相和液相水,如果没有,默认物种是气相的;用井号(#)来作为反应位点,对于多个反应位点可以采用#i来表示,其中i代表第i个反应位点,例如#1和#2分别表示1号和2号反应位点;对于表面物种,采用物种名连接反应位点表示,比如吸附在1号活性位点的水为H2O#1。基元反应中,采用加号(+)分隔反应物种,采用箭头(由小于号、减号和大于号构成<->)分隔正反应逆反应,同时它也作为基元反应输入的标志。

%注释 : 基元反应                          能量输入 [能垒     焓变]   

(1) : H2 + 2#1   <->  2H#1                           [ 0.27  -0.20 ]

(2)  : O2 + 2#1   <->  2O#1                          [ 0.85   0.47 ]

(3) : H#1 + O#2  <->  OH#2 + #1              [ 1.18  -1.09 ]

(4) : H#1 + HO#2 <->  H2O + #1 + #2    [ 1.43  -1.30 ]

以上就是氢气在金属氧化物上燃烧生产水(总反应2H2 + O2 <-> 2H2O)的4步基元反应,其中位点#1代表金属位点,位点#2代表O空位(即O#2代表表面氧位点)。基元反应中冒号(:)前面的内容将只作为注释内容,方括号中的内容为对应的能量输入,具体说明见4.2部分。

4.2 反应能量

有了基元反应步骤,需要进一步提供对应的吉布斯能垒(Ea)和反应自由能(G0),可以输入的能量也可以是反应能垒(Eb)和焓变(H0),然后设置相应的参数来自动进行处理和校正,程序内置了常见分子的热力学焓变、热容、熵(基于Shomate方程)、零点能(基于振动频率)校正的数据,对于复杂的分子,您可以自己手动校正好能量数据直接输入,或者提供相应的热力学参数输入;对于吸附/脱附的能垒,可以选择采用过渡态、碰撞理论、反应热力学等处理方式。

4.2.1 能量输入方式

1 通过设置关键字fun2GetGE= @funName提供自定义的MATLAB函数来计算对应的反应能量信息,其中funName为自定义的MATLAB函数名,可以允许自定义复杂的函数形式,包括考虑覆盖度效应等;

2 直接在基元反应的右侧提供反应能量信息,放在方括号[]中。对于单催化剂模拟,直接提供吉布斯能垒和反应自由能([Ea G0]);对于多催化剂模拟,可以提供基于能量描述符(E1/E2)的能垒和焓变的BEP或线性关系([Ea(scaling) G0(scaling)]/[Ea(BEP) G0(scaling)]),对应的需要提供(E1/E2)的模拟值;

3 直接以变量赋值的方式提供反应能量信息,比如Ea = [Ea1, …],G0 = …。

4.2.2 相关控制参数

1 E1/E2:能量描述符,最多可以允许用两个描述符(E1和E2)来模拟催化剂,设置方式采用赋值方式,比如E1 = [-1:0.1:1],E2 = linspace(-1,1,21);无默认值。

2 Fun2GetGE: 能量输入函数句柄的设置参数,fun2GetGE= @funName,用于(根据描述符)计算对应的反应能量信息,即对应函数形式为[Ea, G0] = funName(E1, E2);默认值@initGE、@initGE1和@initGE2分别对应单催化剂、1个能量和2个能量描述符的情况,需要自己提供对应的函数。

3 fun2E1/fun2E2:基于原能量值E1计算新能量值E1的自定义MATLAB函数(用于考虑覆盖效果等),比如:fun2E1 = @funName,对应函数E1 = funName(E1),或直接使用fun2E1 = @(E1) E1 + 0.2*Q1_I;无默认值。

4 BEPMode: 能垒输入的模式是比例关系还是BEP关系,即Ea =α*E1 + b或α*H0 + b 分别对应0、1;默认值1。

5 ThermoMode:反应物或产物种的热力学校正模式,包括焓变校正(dH)、熵校正(-TdS)、零点能校正(ZPE)及其各种组合,即不校正、dH、dH-TdS、dH-TdS+ZPE、-TdS、-TdS+ZPE、ZPE和dH+ZPE分别对应0、1、2、3、4、5、6和7;默认值0,即不矫正。

6 BarrierMode:吸脱附反应能垒的校正模式,包括不处理、热力学(TD)、过渡态理论(TST)和碰撞理论(CT)处理,分别对应-2、-1、0~1、2;过渡态理论处理时允许设置成0到1的值来考虑不同熵贡献的能垒值;默认值0,即不矫正。

7 ConsMode:能垒限制的模式,包括能垒需要大于0、总能(E0)或自由能(G0)及其组合,即不处理、大于0、E0、G0、(0, E0)、(0, E0)、(E0,G0)或者所有分别对应0、1、2、3、4、5、6、7;默认值0,即不处理。

8 As:吸脱附碰撞理论(CT)处理中必须设置的活性位点面积,单位是埃平方(Å2),对于多位点反应,可以设置成行向量;无默认值。

9 Mr_Species: 吸脱附碰撞理论计算时需要的物种(Species)的相对分子质量,当物种的化学式不是标准写法时,需要自己设定;因为所有元素的相对分子质量都大于1,将其值设为-2、-1和0~1可以将该物种对应的吸脱附过程改为不处理、TD、TST处理;默认值:标准化学式为对应的相对分子质量,非标准时则为-2。

10 自定义或覆盖内置的热力学参数,可采用参数Thermo_Species直接设置物种Species对应的热力学校正参数,可以包括Thermo_Species =

= [dGc];其中,dGc为物种的校正自由能;

= [G, H0, ZPE];其中G = H0 + CpT – TdS; dGc = G – H0;

= [H, H0, CpT, TdS, ZPE, Ef]; 其中Ef 为物种的生成能;

= [A, B, C, D, E, F, G, H, freqi, Ef];其中A-H为shomate参数,freqi为频率。

除内置的常见分子,其余默认值均为0;变温情况下,可增加参数T0作为对应的参考温度;此外还可以新建一个Thermodynamics.data文件,通过关键词CF、TR、AH、FQ和EF设置(分别对应化学式、温度范围、shomate参数、频率和生成能),比如H2的热力学数据设置如下:

CF = H2;                                                                 

EF = [-0.272839]; FQ = [4401];                                          

TR = [ 298, 1000]; AH = [33.07,-11.36,11.43,-2.77, -0.16, -9.98,172.71,0];

TR = [1000, 2500]; AH = [18.56, 12.26,-2.86, 0.27,  1.98, -1.15,156.29,0];

TR = [2500, 6000]; AH = [43.41, -4.29, 1.27,-0.10,-20.53,-38.52,162.08,0];

4.3 反应条件

在动力学的速率表达式中,气相、液相和表面物种(Species)的量(概率)分别以压力(P)、浓度(C)和覆盖度(Q)表示,对于压力和浓度,是无量纲化处理的,即分别除以标准压力(P0: 0.1MPa)和标准浓度(C0: 1mol/L)。物种条件的设置采用X_ Species_Type = value的格式,其中X可以是P(压力)、C(浓度)或Qi(覆盖度)(注:位点#、位点#i和多位点#i#j#k上的物种分别以QQiQixjxk表示(多位点中以x隔开));Species气相、液相和表面物种的名称;Type为需要设置的类型,包括初始化(INIT)、固定化(FROZ)、平衡化(EQUI)和比值化(RATE)处理;value即为设定的参数值。

4.3.1 反应条件初始化(INIT)设置

设置气相、液相和表面物种Species初始压力(P)、浓度(C)和覆盖度(Qi)的初始值,比如P_CO_INIT = 1、C_H2SO4_INIT = 0.1、Q_v_INIT = 1; 默认为0。

4.3.2 反应条件固定化(FROZ)设置

固定气相、液相物种Species压力(P)、浓度(C)的值,被固定的物种,其对应的物料守恒方程dP/dt将从稳态方程中移除,设定的参数值即为稳态时被固定的物种的压力(P)或浓度(C)值。通常可以根据实验条件中反应物及转化率的值,固定所有反应物、产物的压力和浓度,来近似或简化模拟连续流动搅拌反应器(CSTR);否则需要提供实验的条件,包括初始化条件、空速、总压力/浓度等信息(见以下参数Pn、Fgt、Cn、Flt),根据其物料守恒方程来精确求解稳态方程获得。比如P_CO_FROZ = 1、C_H2SO4_FROZ = 0.1;无默认值。

4.3.3 反应步骤平衡化(EQUI)设置

假设某一反应步骤是平衡的来近似处理,即其速率为零;若实际稳态时该反应远离平衡,则该平衡处理结果很可能是无意义的,但可分别对不同反应进行平衡处理,来对活性火山型曲线/面进行分解和分析处理;因为中间体变量与稳态方程是一一对应的,因此新加的平衡处理条件需要替换原来的某稳态方程,所以需要指定替换的稳态方程所对应的中间体,即Qi_Species_EQUI = IndexOfReaction,比如Q­2_OH_EQUI = 4表示第四步基元反应平衡化处理,以其速率等零r(4) = 0替换位点#2上OH表面物种所对应的物料守恒方程dQ2_OH/dt = 0;无默认值。

4.3.4 反应温度(T)设置

单位是开尔文(K),对于单一的反应温度,直接设置其值,如T = 500 K(默认值为室温条件273.15 K);如果是模拟不同反应温度下的情况,可以将其设置为行向量,如T = [300:50:1000]表示温度从300到1000 K每隔50 K布一个点。

4.3.5 连续流动搅拌器(CSTR)模拟

在含非固定化处理的物种的压力(P)或浓度(C)时,需要进一步提供实验的条件,包括空速、总压力、总浓度等信息,否则无法进行求解,包括:

Pn:气相标准化相对总压力,Pn = Ps*Ns/Ng,单位bar

Cn:液相标准化相对总浓度,Cn = Cs*Ns/Nl,单位mol/L

Fgt:气相标准化空速,Fgt = Vgt/Ns,单位m3/s/mol

Flt:液相标准化空速,Flt = Vlt/Ns,单位m3/s/mol

其中,Ns:活性位点数,单位mol

Ps:气相相对总压力,单位bar

Cs:液相相对总浓度,单位mol/L

Ng:气体分子总摩尔数:单位mol

Nl:液体分子总摩尔数;单位mol

Vgt:气相空速,单位m3/s

Vlt:液相空速,单位m3/s

这里采用标准化打包的变量数据Pn、Cn、Fgt、Flt,是因为只有这些标准化的数据改变,才会对计算结果产生影响,因此并不需要全部设置Ns、Ps、Cs、Ng、Nl、Vgt和Vlt所有的参数,而且Pn和Cn只会影响体系到达稳态的(积分)时间(其本质是用来决定反应器的体积的,与下面的布点压力PC相独立),只有Fgt和Flt变化才会影响体系稳态的解;当只有气相或液相时只需要设置对应相的参数即可,不需要全部设置;若空速设为零,则可以计算出反应达到平衡时表面的覆盖度、反应物产物的浓度等;全无默认值。

4.3.6 不同压力或浓度模拟

模拟气相或液相某一或多物种(Species)随不同压力(P)或浓度(C)下的情况,可以通过设置压力P和浓度C参数的行向量值实现,通常情况下,压力和浓度的布点是采用指数方式,例如P = 10.^[-6:1:3],即,压力P从10-6每隔一个数量级布点到103;通过设置对应物种(Species)比值化(RATE)参数,控制哪个或哪些物种参与到压力或浓度的对应改变,例如P_H2_RATE = 2, P_O2_RATE = 1,即H2和O2将随布点压力按2:1对应改变,其余未设置的物种的压力或浓度按原来的设置保持不变;若需要设置的所有比值与初始化(INIT)设置的比值一致,可以不用再单独设置,将自动采用初始化的比值;默认值是0。

4.3.7 敏感度计算(CalcDRC

控制是否进行敏感度(DRC)的计算,包括反应能垒(B)、中间体(I)、压力/浓度(P C)和能量描述符E1/E2(G)的敏感度;其结果矩阵的变量名在datai.mat文件中分别对应DRC_B、DRC_I、DRC_PC和DRC_G;计算敏感度通常需花费要更长的求解时间,特别是基元反应步骤、中间体个数很多的情况下;默认值0;

4.4 控制参数

4.4.1 计算求解相关的参数

1 Ndigits:用于基于可逆度迭代的改进牛顿法计算的有效数字位数,位数越大,求解的稳定性和成功的概率也大,但是耗时可能也越长;默认值是1000。

2 MaxTime:每次改进的牛顿尝试的最大时间,单位秒(s);参数值越大,单次牛顿法求解成功的概率也要高,但是耗时可能更长,相同时间内可以搜索的空间也小;默认值是3。

3 MaxOdeTime:每次时间积分模拟尝试的最大时间,单位秒(s);最大时间越大,求解成功的概率也大,但是耗时可能也越长;默认值是60。

4 TimeMode:是否自动增加每代粒子群算法(PSO)中牛顿迭代的最大时间,分别对应 > 0和0,增加最大时间MaxTime有助于提高求解成功的概率;若大于0,则每代粒子群算法中牛顿迭代的最大时间按以下公式更新:MaxTime = max(MaxTime + 3,(1 + TimeMode)* MaxTime),即每次时间增加TimeMode* MaxTime且每次至少增加3 s;默认值是0.02。

4.4.2 结果分析及可视化的相关参数

1 Tspan:时间积分仿真的时间跨度参数,单位秒(s),时间跨度越大,模拟达到稳态的概率越大,但是耗时可能也越长;默认值是[0 1]。

2 PlotType:是否绘制反应曲线/图,能量分布图,流程图和时间积分模拟,不绘制、绘制曲线/面、势能面图、流程图、时间积分模拟图、所有的以及指定的,分别对应0、1、2、3、4、5 和指定行向量[1 … 4];默认值[1:3]。

3 PlotMode:多反应条件模拟时,是否绘制反应势能面图、流程图和时间积分模拟图,不绘制、第一次、全部绘制或者制定绘制分别对应0、1、inf和指定行向量;默认值是1。

4 ProfMode:以反应总能、标准自由吉布斯能、自由吉布斯能、稳定自由吉布斯能或全部形式绘制反应势能面图,分别对应0、1、2、3和4;默认值是1。

5 PlotList:指定绘制反应结果,可以包括反应速率、选择性、可逆度、覆盖度、表观活化能、以及能垒、中间体、压力/浓度、能量描述符E1/E2的敏感度:Ri、Si_j、Zi、Yi、X[s]i_j、I[s]i_j 、PC[s]i_j、G[s]i_j,其中i /j为指定的下标/活性TOF 下标或它们的组合([s]代表对应参数基于反应物或产物的结果),例如PlotList = {‘R1+1/2*2-3′,’S1+2_3+4′,’Z-2-0.3*1+5′,’Y1+2+3′,’E3′,’B-2-3*1+5 _1+2*2-3’, ‘Is5_2′,’PC1_3′,’Gs1_3’};默认值为空({})。

6 PathOrder:指定的势能面图中反应路径的顺序列表,如PathOrder = {[1 2 3],[1 4]},每个行向量代表一个反应路径,每一个数字代表基元反应的序号;默认值为空({})。

7 PathCoord:指定的势能面图中反应路径的坐标列表,如PathCoord = {[1 2 3],[1 3]},反应坐标与上面的PathOrder一一对应,若为空,则在PathOrder非空时,自动根据物料原则对齐生成坐标;默认值为空({})。

8 PathScale:指定的势能面图中反应路径的系数列表,如PathScale = {[1 2 1/2],[1 -1]};反应系数与上面的PathOrder一一对应;若为空,则在PathOrder非空时,全部系数自动设置成1;默认值为空({})。

9 FigMode:将输出图片格式保存为.fig、.png或全部,分别对应1、2或3;.fig格式可以很方便在MATLAB里面打开查看和进一步编辑处理,.png格式可以用图片软件直接打开查看;默认值是3。

10 SimMode:是否含有反应活性位点以简化反应势能面图及流程图,分别对应0和1;通常包含反应位点的流程图网络通常非常复杂,可阅读和分析性很差,特别是基元反应步骤很多时,默认值是1,即不含反应活性位点。

11 SimValue:反应速率的阈值,用于简化反应的势能面图和流程图;即小于最大反应物和产物速率*SimValue的路径将被忽略,有助于复杂反应网络的主路径提取和分析;默认值是0.05,。

12 CompMode:比较模式,在相同或各自的速率尺度范围内绘制反应势能面图和流程图,分别对应1、0;在相图范围内作图,可以方便比较不同模拟条件下的变化趋势;默认值1。

4.4.3 其他相关参数

1 Istart:是否加载之前的参数设置和稳态方程等数据(ORGi.mat)再继续运行,分别对应1、0;加载时仅更新控制参数部分;若有未完成的任务,设为不加载时将根据INPUT所有参数重新初始化、产生基元反应表达式和稳态方程等,并且如果反应机理、反应能量或反应条件改变时,请删除对应任务的输出文件,特别是doneid_i、data*.mat等,否则,前后两次运行整合了反应机理、反应能量或反应条件不一致的结果;默认值是1。

2 npar:用于平行计算的线程数,>0 时为平行计算,临时计算文件在pari下,运行结束后将自动清空并汇总,并行结果将被汇总到datai.mat和logi_parsum文件中;并产生figi/Prediction文件夹,含有根据部分获得计算的结果,利用机器学习预测完整空间的覆盖度和活性图,可以提前预判模拟是否准确及其规律;默认值是0,即1个线程串行计算。

3 runid:指定到某一次运行计算,若该任务不存在,则开始一个新的任务;存在未完成,则继续计算任务;存在且已完成,则加载数据仅更新控制参数重新进行结果分析及可视化任务;默认值是第一个未完成的任务,若无未完成任务则开始一个新的任务。

4 SaveMode:将所有数据保存为双精度(double)或高精度符号(symbolic)类型,对应0和1;保存双精度速度非常快,需要的内存和硬盘非常小;保存成高精度符号结果速度慢,需要的内存和硬盘大,特别是布点数量很多的时候;默认值是0;

5 SaveFreq:保存计算结果的频率(串行模式下),每隔多少个运行样本保存一次,参数越大,保存的频率及耗时越小,若计算速度很快、计算条件好(内存、容量大,即不会出现内存、容量不足的情况)且稳定(不会出现MATLAB程序、计算机卡死、断电等计算异常中断的情况),可以尽量考虑设大一些,特别是将数据保存为double类型时;默认值是1。

6 SkipMode:从不过滤改进牛顿方法找到的解,或者过滤包含及全是零(全是零的解指每一个位点上只有某一物种是1,其余的均为零的情况)的解,分别对应0、1和2;通常情况下,包含零的解是没有物理意义的,所以需要过滤掉,但是如果您设置的反应机理、反应能量或反应条件就太离谱、不切合实际,则该模拟并不能保证存在有物理意义的解;试图进行这样的模拟完全是浪费时间、毫无意义的,所以模拟前、求解时、甚至是求解结束,请再三检查和确认您设置的反应机理、反应能量或反应条件等参数、软件运行的输出信息等(特别是警告信息);默认值是0。

7 CheckMode:是否检查反应的热力学,不检查、第一次或者全部检查分别对应0、1和2;如果您的反应内循环的热力学不守恒(通常是同一物种,在不同基元反应中考虑共吸附方式不同导致的,当然也有可能是您的数据处理出错),则可能出现一些奇怪、不合实际的结果,比如某些反应路径的方向等,特别是他们被包含在主反应路径中时;默认值是2。

8 CalMode:计算模拟的模式,通常情况下不需要单独设置,可以根据模拟的条件自己判断;在目前的版面里,可以允许考虑布点的变量包括能量描述符E1/E2、温度T、压力P、浓度C、气相空速Fgt和液相空速Flt(分别以012345标记),具体对应关系说明如下:

i. CalMode = 0,对应单条件模拟情况;

ii. CalMode >= 1,对应一维条件模拟情况;具体对应关系为

E1: 1.0、T: 1.1、P: 1.2、C: 1.3、Fgt: 1.4、Flt: 1.5

iii. CalMode >= 2,对应二维条件模拟情况;具体对应关系为

(E1, E2): 2.00、(E1, T): 2.01、(E1, P): 2.02、(E1, C): 2.03、(E1, Fgt): 2.04、(E1, Flt): 2.05、(T  , P): 2.12、(T  , C): 2.13、(T , Fgt): 2.14、(T  , Flt): 2.15、(P , C  ): 2.23、(P, Fgt): 2.24、(P, Flt): 2.25、(C , Fgt): 2.34、(C , Flt): 2.35、(Fgt, Flt): 2.45