一种基于时间积分和牛顿法自动联用高效求解稳态微观动力学方程组的方法技术

本发明专利技术公开了一种基于时间积分和牛顿法自动联用高效求解稳态微观动力学方程组的方法，采用基于低精度数值的时间积分法对初始值进行时间积分，结合覆盖度敏感度自动积分到粗精度的稳态覆盖度，再联用高精度数值的牛顿法，获得稳态覆盖度；若单次求解失败，可在指数空间随机初始化覆盖度来进行新一轮的求解，用于自动化求解催化剂表面化学反应网络的微观动力学方程组。本发明专利技术的方法对求解稳态微观动力学方程组具有更高的成功率、速度和精度；同时本发明专利技术自动化程度高，只需要设定初始参数，就可以实现自动化迭代求解，过程不需要人机互动，而且算法流程简单，容易上手。

【技术实现步骤摘要】
一种基于时间积分和牛顿法自动联用高效求解稳态微观动力学方程组的方法
本专利技术属于多相催化反应微观动力学分析研究领域，具体涉及一种基于时间积分和牛顿法自动联用高效求解稳态微观动力学方程组的方法，可用于复杂催化剂表面反应网络微观动力学方程的高精度、快速、自动化求解。
技术介绍
多相催化在现代化工中有着不可或缺的作用，尤其是在能源转化和环境领域，其中催化剂扮演着关键角色。现如今，高效催化剂的科学设计已成为研究人员的重要目标，但也是挑战难题。相比于传统实验上相对低效的试错法，近二十多年来密度泛函理论(DFT)和微观动力学分析的发展，可以促使研究人员从理论上评估催化剂的活性及变化规律，进而有助于催化剂的科学筛选和整体设计。对于复杂的催化反应体系而言，微观化学基元反应众多，之间互相链接构成了复杂的动力学体系。DFT可以从微观层面计算得到每一步基元反应的反应能垒和焓变，以及中间物种的吸附能。基于平均场近似的微观动力学则在微观化学基元反应和宏观催化剂反应活性之间扮演了桥梁作用。通过微观动力学分析，可以得到反应的宏观性质(覆盖度、反应速率等)以及决速因素等，进而来描述催化剂的活性趋势。基于DFT计算数据和微观动力学概念，根据质量作用定律写出基元反应的反应速率表达式，可以推导出化学反应在催化剂表面的速率表达式，然后针对每一个吸附物种，依据物料守恒推导出催化剂表面吸附物种覆盖度随时间变化的一组常微分方程组，最后依据稳态近似(即当反应体系达到稳态时，覆盖度随时间的变化率为0)得到稳态方程组，通过迭代求解可以得到稳态条件下反应的反应速率(催化活性)以及表面中间物种覆盖度等数据，具体流程见图1所示。然而，由于微动力学方程组的复杂性(基元反应多、中间物种多、反应级数不同)和基元反应速率常数之间数量级的差异等引发的刚性问题(stiffproblem)，一般很难快速地得到其精确解。目前，微观动力学模拟程序CHEMKIN、CatMAP等采用改进牛顿法进行求解。牛顿法因其收敛速度快而被广泛应用于微观动力学方程组的求解，但需预先给定一个初始值进行迭代求解；当达到给定收敛标准时(可认为体系达到稳态)便终止迭代。求解结果的质量(求解速度、收敛精度和解的合理性)高度依赖于给定的初始值。因此，该方法只有选用合理的初始值，才能快速地得到精确值，否则很难收敛到高精度的解、或者有效解(具有物理意义)；而且由于刚性问题，需要采用高精度数值的计算，才能保证其迭代方向(雅可比矩阵)的有效性和准确性，使得牛顿法迭代通常需要消耗大量的计算机内存和计算时间，特别是在牛顿法求解失败的时候。除此之外，微观动力学模拟程序Cantera、MKMCXX等采用高精度数值的时间积分到足够长的时间尺度作为稳态解，高精度数值的时间积分法虽然可以通过迭代达到相对稳定的解，但通常需要很大的时间步长和计算机内存，并且反复迭代，导致求解时间很长；低精度数值的时间积分法虽然更快速，但通常只能得到覆盖度的大概分布趋势，其数量级跟高精度数值结果差别较大，特别是在刚性程度较高的情况下。若简单地结合时间积分和牛顿法交替迭代，虽然有可能实现结合两种方法特点的可能，但采用未达到稳态的覆盖度作为牛顿法的初始值，导致求解结果具有很大的失败概率，同时往返的多次迭代调牛顿法致使需要更长的求解时间，因此，在实际的应用中，该方法依然存在求解时间长、效率低和经常失败等问题。
技术实现思路
为了克服现有的求解稳态微观动力学方程组求解方法的不足，本专利技术的目的是提供一种新型的基于时间积分和牛顿法自动联用的方法来实现对稳态动力学方程组的高效求解的方法。其简要思想为：本专利技术结合了时间积分法可稳定迭代以及牛顿法快速收敛的优点，采用基于低精度数值(16位有效数字)的时间积分法对初始值进行时间积分，结合覆盖度敏感度可快速、自动积分到粗精度的稳态覆盖度，再联用高精度数值(200位有效数字以上)的牛顿法，获得准确、高精度的稳态覆盖度；若单次求解失败，采用指数空间随机初始化覆盖度来进行新一轮的求解，保证搜索空间的广度与多样性，以提高其求解成功率。为了实现上述目的，本专利技术采用的技术方案如下：为了求解稳态动力学方程组，首先通过低精度数值(16位有效数字)的时间积分法对反应初始覆盖度θ0进行时间积分，并快速得到粗精度下的稳态覆盖度θi，使其作为高精度数值(200位有效数字)牛顿法的初始值求解微观动力学稳态方程；当达到收敛标准时(稳态条件)，可得到反应速率以及表面物种覆盖度等宏观性质，进而评价催化剂的反应催化活性或为进一步改进催化剂的催化活性作参考。本专利技术的第一方面提供了一种基于时间积分和牛顿法自动联用高效求解稳态微观动力学方程组的方法，采用基于低精度数值(16位有效数字)的时间积分法对初始值进行时间积分，结合覆盖度敏感度自动积分到粗精度的稳态覆盖度，再联用高精度数值(200位有效数字以上)的牛顿法，获得稳态覆盖度；若单次求解失败，可在指数空间随机初始化覆盖度来进行新一轮的求解，用于自动化求解催化剂表面化学反应网络的微观动力学方程组，包括以下步骤：步骤0：输入催化反应体系的微观动力学方程组和反应表面物种(1，…，m)的初始覆盖度行向量θ0＝[θ0,1；…；θ0,m]；步骤1：利用低精度数值(16位有效数字)的时间积分法，将常微分方程组在时间尺度上从0积分到ti，对应的覆盖度行向量从θ0积分到θi，得到过程中覆盖度随时间的变化，时间变化列向量为t＝[t0,…,ti]，对应的覆盖度变化矩阵为θ＝[θ0,…,θi]；步骤2：求表面物种m在时间ti的覆盖度敏感度即在对数空间覆盖度随时间的变化斜率，实际求解时采用其差分形式近似根据覆盖度敏感度矩阵X的范数norm(X)是否小于给定的误差(精度)ε，判断其是否趋于稳态；如果是，转到步骤3；否则，更新覆盖度θ0＝θi，提高时间尺度ti＝k×ti，再重复步骤1、2；k表示时间提高系数；步骤3：利用趋于稳态的粗精度覆盖度θi作为高精度数值(200位有效数字以上)牛顿法的初始值，求解稳态微观动力学方程组高精度的覆盖度θs；步骤4：判断牛顿法求解是否成功收敛，如果是，转到步骤5；否则，初始化时间尺度ti，随机初始化覆盖度θ0，再重复步骤1、2、3、4。步骤5：结束。所述步骤1中常微分方程组由表面各吸附物种覆盖度随时间变化的表达式组成，可由步骤0中微观动力学方程组推导得到。所述步骤1中覆盖度随时间的变化表示为：dθi/dt＝∑vjrj，其中rj表示基元反应j的反应速率，vj表示反应速率rj前面的系数；对应的时间积分形式表示为：所述步骤1中的初始时间尺度ti的最优值在10-15到10-8之间。所述步骤1中的时间积分法采用低精度数值(16位有效数字)的刚性常微分方程组求解，具体如MATLAB刚性常微分方程组求解器ode15s、ode23s、ode23t、ode23tb等中的一种。所述步骤2中的覆盖度敏感度这一概念由本专利技术提出；式中，t表示时间，θ表示覆盖度。所述步骤2中的log()为取对数函数，其底数可以是本文档来自技高网...

【技术保护点】
1.一种基于时间积分和牛顿法自动联用高效求解稳态微观动力学方程组的方法，其特征在于，包括以下步骤：/n步骤0：输入催化反应体系的微观动力学方程组和反应表面物种(1，…，m)的初始覆盖度行向量θ

【技术特征摘要】
1.一种基于时间积分和牛顿法自动联用高效求解稳态微观动力学方程组的方法，其特征在于，包括以下步骤：
步骤0：输入催化反应体系的微观动力学方程组和反应表面物种(1，…，m)的初始覆盖度行向量θ0＝[θ0,1；…；θ0,m]；
步骤1：利用低精度数值的时间积分法，将常微分方程组在时间尺度上从0积分到ti，对应的覆盖度行向量从θ0积分到θi，得到过程中覆盖度随时间的变化，时间变化列向量为t＝[t0,…,ti]，对应的覆盖度变化矩阵为θ＝[θ0,…,θi]；
步骤2：求表面物种m在时间ti的覆盖度敏感度即在对数空间覆盖度随时间的变化斜率，实际求解时采用其差分形式近似根据覆盖度敏感度矩阵X的范数norm(X)是否小于给定的误差ε，判断其是否趋于稳态；如果是，转到步骤3；否则，更新覆盖度θ0＝θi，提高时间尺度ti＝k×ti，再重复步骤1、2；k表示时间提高系数；
步骤3：利用趋于稳态的粗精度覆盖度θi作为高精度数值牛顿法的初始值，求解稳态微观动力学方程组高精度的覆盖度θs；
步骤4：判断牛顿法求解是否成功收敛，如果是，转到步骤5；否则，初始化时间尺度ti，随机初始化覆盖度θ0，再重复步骤1、2、3、4。
步骤5：结束。


2.根据权利要求1所述的基于时间积分和牛顿法自动联用高效求解稳态微观动力学方程组的方法，其特征在于，所述步骤1中常微分方程组由表面各吸附物种覆盖度随时间变化的表达式组成，可由步骤0中微观动力学方程组推导得到；
所述步骤1中覆盖度随时间的变化表示为：dθi/dt＝∑vjrj，其中rj表示基元反应j的反应速率，vj表示反应速率rj前面的化学计量系数；对应的时间积分形式表示为：


3.根据权利要求1所述的基于时间积分和牛顿法自动联用高效求解稳态微观动力学方程组的方法，其特征在于，所述步骤1中的初始时间尺度ti在10-15到...

【专利技术属性】
技术研发人员：王海丰，陈建富，来壮壮，胡培君，
申请(专利权)人：华东理工大学，
类型：发明
国别省市：上海;31