CN113724801A

CN113724801A - 分解反应的分析方法

Info

Publication number: CN113724801A
Application number: CN202111010863.3A
Authority: CN
Inventors: 李艳平; 严大洲; 杨涛; 温国胜; 韩治成; 姚心; 刘诚; 万烨; 赵雄; 司文学; 孙强; 张升学
Original assignee: China ENFI Engineering Corp
Current assignee: China ENFI Engineering Corp
Priority date: 2021-08-31
Filing date: 2021-08-31
Publication date: 2021-11-30

Abstract

本公开涉及一种分解反应的分析方法，包括：建立分解反应的模拟体系，并设置模拟体系的反应条件；确定模拟体系中各粒子的属性参数；根据各粒子的属性参数和模拟体系的反应条件，确定各粒子的初始速度和初始位置；对模拟体系进行划分形成多个进程，对各进程进行划分，形成硬球模型区域和拟颗粒模型区域，并建立硬球模型区域的未来事件列表；基于未来事件列表处理硬球模型区域内粒子间的碰撞和反应，基于反应的时间步处理拟颗粒模型区域内粒子间的碰撞和反应；统计反应物粒子、中间产物粒子和生成物粒子信息，生成统计结果，当模拟体系中的反应物粒子的数量为零时，统计结果为分解反应的分析结果。该分析方法能显著提高对分解反应的分析速度和效率。

Description

分解反应的分析方法

技术领域

本公开涉及分子动力学领域，尤其涉及一种分解反应的分析方法。

背景技术

目前，在分子动力学领域，通常采用高斯原理和第一性原理等量化计算方法对电子层面上的微观化学反应进行计算，但是采用该方法对全模拟体系的演化进行模拟和分析时，通常需要巨大的计算量。并且，采用该方法进行分析时，需要对反应过程中化学键的形成和断裂进行建模，并且在每个时间步均需要使模拟体系达到电荷的动态平衡，从而导致其计算速度较慢。其次，采用该方法进行分析的时候，其能够计算的原子总数量较小，从而适得其计算的规模较小，进而导致分析的准确性较低。最后，采用该方式进行分析时，只能够针对特定的反应过程，从而导致其分析的通用性较差。

需要说明的是，在上述背景技术部分公开的信息仅用于加强对本公开的背景的理解，因此可包括不构成对本领域普通技术人员已知的现有技术的信息。

发明内容

本公开目的在于提供一种分解反应的分析方法，该分析方法能够显著提高对分解反应的分析速度和分析效率，并且具有较好的通用性。

本公开提供了一种分解反应的分析方法，包括：

建立所述分解反应的模拟体系，并设置所述模拟体系的反应条件；

确定所述分解反应的模拟体系中各粒子的属性参数；

根据所述各粒子的属性参数和所述模拟体系的反应条件，确定所述各粒子的初始速度和初始位置；

对所述模拟体系进行划分以形成多个进程，并对每个所述进程进行划分，以形成硬球模型区域和拟颗粒模型区域，并建立所述硬球模型区域的未来事件列表；

基于所述未来事件列表处理所述硬球模型区域内粒子间的碰撞和反应，并基于反应的时间步处理所述拟颗粒模型区域内粒子间的碰撞和反应粒；

统计所述模拟体系中的反应物粒子、中间产物粒子和生成物粒子信息，并生成统计结果，当所述模拟体系中的反应物粒子的数量为零时，所述统计结果为所述分解反应的分析结果。

在本公开的一种示例性实施例中，所述确定所述分解反应的模拟体系中各粒子的属性参数包括：

对各所述粒子进行编号，设置各所述粒子的质量、计算各所述粒子硬球半径，并对各所述粒子的质量和硬球半径进行无量纲化处理，以得到各所述粒子的属性参数。

在本公开的一种示例性实施例中，所述建立所述分解反应的模拟体系，并设置所述模拟体系的反应条件，包括：

根据所述分解反应的类型，确定所述模拟体系在各方向上的尺寸；

根据所述模拟体系在各方向上的尺寸，建立所述模拟体系；

设置所述模拟体系的反应温度、反应活化能和每种编号的粒子的数量，以形成所述模拟体系的反应条件。

在本公开的一种示例性实施例中，所述根据所述各粒子的属性参数和所述模拟体系的反应条件，确定所述各粒子的初始速度，包括：

根据所述反应温度和所述各粒子的质量，计算所述各粒子的热速度；

设定所述各粒子的流动速度和所述各粒子的运动角度；

根据所述各粒子的热速度、流动速度和运动角度，确定所述各粒子的初始速度。

在本公开的一种示例性实施例中，所述对所述模拟体系进行划分以形成多个进程，并对每个所述进程进行划分，以形成硬球模型区域和拟颗粒模型区域，包括：

根据所述模拟体系，将所述模拟体系划分为多个进程；

确定每个所述进程的边界位置和中间位置；

将所述边界位置进行划分，以将所述边界位置分隔为多个边界网格区；

将所述中间位置进行划分，以将所述中间位置分隔为多个中间网格区；设置所述边界网格区为所述拟颗粒模型区域，并且设置所述中间网格区为所述硬球模型区域。

在本公开的一种示例性实施例中，基于未来事件列表处理所述硬球模型区域内粒子间的碰撞和反应，并基于反应的时间步处理所述拟颗粒模型区域内粒子间的碰撞和反应，包括：

根据所述未来事件列表里的事件发生时间的先后顺序，判断所述硬球模型区域中的两个即将发生碰撞的各粒子是否发生反应；

当所述硬球模型区域中的两个即将发生碰撞的粒子不发生反应时，则按照正常的碰撞事件进行处理，并更新所述未来事件列表；

当所述硬球模型区域中的两个即将发生碰撞的粒子发生反应时，根据所述反应物粒子发生的反应，更新所述未来事件列表和所述硬球模型区域中的反应物粒子、中间产物粒子和生成物粒子的种类和数量；

对所述拟颗粒模型区域和所述硬球模型区域内的所有粒子均添加第一标记；

在当前反应的时间步中，处理所述拟颗粒模型区域内各粒子间的碰撞，并将即将分解的粒子的第一标记修改为第二标记；

在所述当前反应的时间步后，并在下一个反应的时间步前，首先检测所述拟颗粒模型区域内各个粒子的标记，将添加有所述第二标记的粒子分解，以形成所述中间产物粒子和所述生成物粒子。

在本公开的一种示例性实施例中，所述硬球模型区域中的各粒子中至少包括所述反应物粒子，所述根据所述各粒子发生的反应，更新所述未来事件列表和所述硬球模型区域中的反应物粒子、中间产物粒子和生成物粒子的种类和数量，包括：

根据所述反应物粒子发生的反应，将模拟体系中发生反应的所述反应物粒子转化为一个中间产物粒子，并生成另一个中间产物粒子；

删除所述未来事件列表中关于所述反应物粒子的未来事件，并增加关于所述中间产物粒子的未来事件，以更新所述未来事件列表；

根据更新的所述未来事件列表，对所述中间产物粒子的碰撞进行预测和处理，并更新所述中间产物粒子的半径；

当所述中间产物粒子的半径与所述生成物粒子的半径相同时，所述中间产物粒子为新生成的所述生成物粒子；

根据转化的所述反应物粒子和所述新生成的生成物粒子更新所述硬球模型区域中的反应物粒子、中间产物粒子和生成物粒子种类和数量。

在本公开的一种示例性实施例中，所述建立所述硬球模型区域的未来事件列表，包括：

设定第一时间步，并对所述硬球模型区域在每个所述第一时间步中的未来事件进行预测；

将预测的所述未来事件按照事件发生时间的先后顺序进行排序，以形成所述硬球模型区域的未来事件列表。

在本公开的一种示例性实施例中，对所述硬球模型区域在每个所述第一时间步中的未来事件进行预测，包括：

在每个所述第一时间步中，对所述硬球模型区域的两个即将发生碰撞的粒子的碰撞时刻进行预测；

对所述粒子穿越网格边界的时刻、所述中间产物粒子的半径变化为所述生成物粒子半径的时刻进行预测；

其中，对所述硬球模型区域的两个即将发生碰撞的粒子的碰撞时刻进行预测包括：

根据所述两个即将发生碰撞的粒子的初始速度、初始粒子半径、初始位置和粒子半径增长速度，计算所述两个即将发生碰撞的粒子发生碰撞的时间；

根据所述两个即将发生碰撞的粒子发生碰撞的时间，预测所述两个即将发生碰撞的粒子的碰撞时刻。

在本公开的一种示例性实施例中，所述两个即将发生碰撞的粒子发生碰撞的时间为：

其中，

△t为所述两个即将发生碰撞的粒子发生碰撞的时间，

为第一个粒子的初始粒子半径，

为第二个粒子的初始粒子半径，v_r1为第一个粒子的粒子半径增长速度，v_r2为第二个粒子的粒子半径增长速度，r₁ ⁰为第一个粒子的初始位置，

为第二个粒子的初始位置，V₁为第一个粒子的初始速度，V₂为第二个粒子的初始速度。

本公开提供的技术方案可以达到以下有益效果：

本公开所提供的分解反应的分析方法，能够将模拟体系划分为多个进程，并针对多个进程分别进行计算，从而能够实现大规模、多进程的并行模拟计算，从而使得本公开所提供的分析方法具有较大的计算规模，进而使得使用该分析方法分析的准确性较高。

其次，本公开所提供的分析方法，将每个进程划分为硬球模型区域和拟颗粒模型区域，从而能够对不同的区域采用不同的计算方法，因此该分析方法不需要同步演化模拟体系中的所有粒子，也就能够节省大量的计算时间，大幅度提高计算效率。

另外，本公开所提供的分析方法，将模拟体系中的各个粒子假设为无内部结构的硬球，从而本公开不需要考虑各个粒子的动态电荷信息，也就因此减少了大量的计算量。并且，本公开所提供的分析方法在分析的过程中，只需要各个粒子的硬球半径、速度和位置等信息，从而使得该分析方法具有较好的通用性。

应当理解的是，以上的一般描述和后文的细节描述仅是示例性和解释性的，并不能限制本公开。

附图说明

此处的附图被并入说明书中并构成本说明书的一部分，示出了符合本公开的实施例，并与说明书一起用于解释本公开的原理。显而易见地，下面描述中的附图仅仅是本公开的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。

图1示出了根据本公开一示例性实施例的分解反应的分析方法的流程示意图；

图2示出了根据本公开一示例性实施例的两个生成物粒子与反应物粒子之间的位置关系的示意图。

具体实施方式

现在将参考附图更全面地描述示例实施方式。然而，示例实施方式能够以多种形式实施，且不应被理解为限于在此阐述的实施方式；相反，提供这些实施方式使得本公开将全面和完整，并将示例实施方式的构思全面地传达给本领域的技术人员。图中相同的附图标记表示相同或类似的结构，因而将省略它们的详细描述。

用语“一个”、“一”、“该”、“所述”用以表示存在一个或多个要素/组成部分/等；用语“包括”和“具有”用以表示开放式的包括在内的意思并且是指除了列出的要素/组成部分/等之外还可存在另外的要素/组成部分/等；用语“第一”和“第二”等仅作为标记使用，不是对其对象的数量限制。

本公开首先提供了一种分解反应的分析方法，该分解反应的分析方法相比于现有技术来说，具有较大的计算规模、较高的准确性、较高的计算效率、较少的计算量和较好的通用性。具体的，如图1所示，该分解反应的分析方法可以包括：

步骤S10、建立分解反应的模拟体系，并设置模拟体系的反应条件；

步骤S20、确定所述分解反应的模拟体系中各粒子的属性参数；

步骤S30、根据各粒子的属性参数和模拟体系的反应条件，确定各粒子的初始速度和初始位置；

步骤S40、对模拟体系进行划分以形成多个进程，并对每个进程进行划分，以形成硬球模型区域和拟颗粒模型区域，并建立硬球模型区域的未来事件列表；

步骤S50、基于未来事件列表处理硬球模型区域内粒子间的碰撞和反应，并基于反应的时间步处理拟颗粒模型区域内粒子间的碰撞和反应；

步骤S60、统计模拟体系中的反应物粒子、中间产物粒子和生成物粒子信息，并生成统计结果，当模拟体系中的反应物粒子的数量为零时，统计结果为分解反应的分析结果。

下面对上述各个步骤进行详细说明：

在上述步骤S10中，可以建立分解反应的模拟体系，并设置模拟体系的反应条件。具体的，可以根据分解反应的类型，确定模拟体系在各方向上的尺寸，再根据模拟体系在各方向上的尺寸，建立模拟体系。举例而言，该模拟体系可以为沸石ZSM-5体系，该沸石ZSM-5体系具有多个沸石晶胞和孔道，反应物粒子可以分布在孔道内。但不限于此，也可以建立其他模拟体系，可以根据实际需要设置。

进一步的，可以设置模拟体系的反应温度、反应活化能和每种类型的粒子的数量，以形成模拟体系的反应条件。同时，也可以对反应活化能进行无量纲化处理。并且，可以将模拟体系的反应温度、无量纲化处理后的反应活化能、每种类型的粒子的数量写入参数文件中。

另外，还可以设置模拟体系的进程数。由于在后续计算的过程中，一个进程负责一个大的区域的计算。从而，当进程数较多时，能够显著提高计算的速度和效率。同时，由于相邻的进程间可以通过消息传递接口交换边界信息，从而能够实现大规模、多进程的并行计算，使得计算规模显著增加，计算规模可以达到10000数量级。当本公开对模拟体系的进程数设置完成之后，也可以将该进程数写入参数文件中。

在步骤S20中，可以确定分解反应的模拟体系中各粒子的属性参数。具体的，该分解反应的模拟体系中可以具有反应物粒子、中间产物粒子和生成物粒子，也可以仅仅具有反应物粒子，或者具有反应物粒子和中间产物粒子或反应物粒子和生成物粒子。

进一步的，可以对各粒子进行编号，设置各粒子的质量，计算各粒子的硬球半径，并对各粒子的质量和硬球半径进行无量纲化处理，以得到各粒子的属性参数。

具体地，可以对反应物粒子进行编号，设置反应物粒子的质量，计算反应物粒子的硬球半径，并对反应物粒子的质量和硬球半径进行无量纲化处理，以得到反应物粒子的属性参数。同时，也可以对生成物粒子进行编号，设置生成物粒子的质量，计算生成物粒子的硬球半径，并对生成物粒子的质量和硬球半径进行无量纲化处理，以得到生成物粒子的属性参数。当模拟体系中存在中间产物粒子的时候，也可以对中间产物粒子进行编号，设置中间产物粒子的质量，计算中间产物粒子的硬球半径，并对中间产物粒子的质量和硬球半径进行无量纲化处理，以得到中间产物粒子的属性参数。

在本公开的一个实施例中，当分解反应的反应物只有一个时，且生成物粒子为两个时，反应物粒子的编号可以为A，生成物粒子的编号可以为B和C，中间产物粒子的编号可以为B’和C’。但不限于此，当分解反应的生成物粒子为三个时，生成物粒子的编号可以为B、C和D，以此类推。进一步的，为了区分模拟体系中的各个反应物粒子，可以对各个反应物粒子进行编号，将其编为1、2、3…n。同时，为了更进一步的对同种生成物粒子中的各个生成物粒子进行区分，可以对同种生成物粒子中的各个生成物粒子进行编号，例如：可以对编号为B的生成物粒子中的各个生成物粒子进行编号，可以将其编为1、2、3…n等。需要说明的是，本公开对各种编号的形式不做限定，也可以不使用英文字母或者阿拉伯数字进行编号，例如：可以使用罗马数字对其进行标号等，这也在本公开的保护范围之内。

进一步的，本公开通过对反应物粒子的质量、生成物粒子的质量、中间产物粒子的质量、反应物粒子的硬球半径、生成物粒子的硬球半径和中间产物粒子的硬球半径进行无量纲化处理，能够便于后续的分析和计算，使得后续在计算的过程中所有的物理量均为无量纲的单位，进而能够在后续计算中大幅度减少单位换算的计算量，也就因此能够使得本公开的计算量大幅度减少，计算效率大幅度提升。

在本公开的一个实施例中，可以采用极值化法、标准化法、均值化法、标准差化法对上述各个物理量进行无量纲化处理。但不限于此，也可以使用其他方法对上述各个物理量进行无量纲化处理，这均在本公开的保护范围之内。

另外，本公开还可以对粒子的特征时间进行确定，并对该特征时间也同样进行无量纲化处理。

当对本公开对各粒子的属性参数通过上述方式进行确定了以后，可以将确定好的各粒子的属性参数整体写入分析计算的参数文件中，以作为后续计算的基本计算参数，并可以在后续计算中按照需要在参数文件中找到和调取计算所需的任意参数。

在步骤S30中，可以根据各粒子的属性参数和模拟体系的反应条件，确定各粒子的初始速度和初始位置。具体的，可以根据模拟体系的反应温度和各粒子的质量，计算各粒子的热速度。其中，各粒子在立体模拟体系中满足：

其中，v₀为粒子的热速度，k_B为玻尔兹曼常数，T为反应温度，m为粒子的质量。

另外，各粒子在平面模拟体系中满足：

进一步的，可以设定各粒子的流动速度和各粒子的运动角度。需要说明的是，各粒子的运动角度是通过随机数生成的，本公开对各粒子的流动速度和各粒子的运动角度不做限定，可以根据分析的实际需要进行设定，这均在本公开的保护范围之内。

在本公开的一个实施例中，可以根据各粒子的热速度、流动速度和运动角度，确定各粒子的初始速度。在立体模拟体系中，可以设置相互垂直的三个方向X、Y和Z方向。通过各向同性的球形分布可以得到：

其中，v_x为粒子在X方向上的初始速度，v_y为粒子在Y方向上的初始速度，v_z为粒子在Z方向的初始速度，v_rx为粒子在X方向上的流动速度，v_ry为粒子在Y方向上的流动速度，v_rz为粒子在Z方向上的流动速度，v₀为粒子的热速度，

为粒子的速度方向与X方向之间的方位角，θ为粒子的速度方向与Z方向之间的仰角，其中，

与θ由设置的时间种子随机地生成。

在平面模拟体系中，可以设置相互垂直的X和Y方向，即可得到：

其中，v_x为粒子在X方向上的初始速度，v_y为粒子在Y方向上的初始速度，v_rx为粒子在X方向上的流动速度，v_ry为粒子在Y方向上的流动速度，v₀为粒子的热速度，

为粒子的速度方向与X方向之间的夹角。

需要说明的是，上述所涉及的方位角

的取值范围为[0，2π]，仰角θ的取值范围为[0，π]。

在本公开的一个实施例中，可以根据模拟体系中的孔道形状确定反应物粒子的初始位置。举例而言，可以将各粒子随机、均匀的分布在模拟体系的孔道内，并且使得任意两个粒子之间不发生重叠即可，但不限于此，也可以采用其他的方式来设置各粒子的初始位置，这均在本公开的保护范围之内。

在步骤S40中，可以对模拟体系进行划分以形成多个进程，并对每个进程进行划分，以形成硬球模型区域和拟颗粒模型区域，并建立硬球模型区域的未来事件列表。

具体的，可以根据模拟体系，将模拟体系划分为多个进程。当本公开在上述步骤中设置有进程数时，可以根据上述设置的进程的数量将模拟体系划分为多个进程。并且，可以确定每个进程的边界位置和中间位置。可以将边界位置进行划分，以将边界位置分隔为多个边界网格区，并且可以将中间位置进行划分，以将中间位置分隔为多个中间网格区。进一步的，可以设置边界网格区为拟颗粒模型区域，并且可以设置中间网格区为硬球模型区域。

在本公开的一个实施例中，可以设定第一时间步，并对硬球模型区域在每个第一时间步中的未来事件进行预测。具体地，可以在每个第一时间步中，对硬球模型区域的两个即将发生碰撞的粒子的碰撞时刻、所述粒子穿越网格边界的时刻、所述中间产物粒子的半径变化为所述生成物粒子半径的时刻进行预测。其中，对硬球模型区域的两个即将发生碰撞的粒子的碰撞时刻进行预测可以包括：可以根据两个即将发生碰撞的粒子的初始速度、初始粒子半径、初始位置和粒子半径增长速度，计算两个即将发生碰撞的粒子发生碰撞的时间。并且可以根据两个即将发生碰撞的粒子发生碰撞的时间，预测两个即将发生碰撞的粒子的碰撞时刻。

在本公开的一个实施例中，由于模拟体系中硬球模型区域内任意两个粒子的位置没有发生重叠，从而两个粒子需要满足的碰撞条件是两个粒子的质心之间的距离等于两个粒子碰撞时的半径之和。即可得到：

其中，

为第一个粒子碰撞时的位置，

为第二个粒子碰撞时的位置，

为第一个粒子碰撞时的半径，

为第二个粒子碰撞时的半径。其中，需要说明的是，上述各个参量均为向量。

在本公开的一个实施例中，第一个粒子碰撞时的位置和第二个粒子碰撞时的位置可以为：

其中，r₁ ⁰为第一个粒子的初始位置，

为第二个粒子的初始位置，V₁为第一个粒子的初始速度，V₂为第二个粒子的初始速度，△t为所述两个即将发生碰撞的粒子发生碰撞的时间。

进一步的可以设置粒子半径增长参数，通过设置粒子半径增长参数能够避免新生成的粒子与周围粒子的位置发生重叠。因此，上述第一个粒子碰撞时的半径可以为：

其中，

为第一个粒子的初始粒子半径，v_r1为第一个粒子的粒子半径增长速度。

上述第二个粒子碰撞时的半径可以为：

其中，

为第二个粒子的初始粒子半径，v_r2为第二个粒子的粒子半径增长速度。

由于两个粒子发生碰撞的条件是两个粒子必须存在相互碰撞的趋势，即要求满足上述的计算公式有解。从而，根据求根公式可以得到：

需要说明的是，为了保证上述计算公式有解，需要满足△/4>0。

进一步的，可以根据上述公式以及△/4>0的条件，即可以得到两个即将发生碰撞的粒子发生碰撞的时间为：

其中，△t为两个即将发生碰撞的粒子发生碰撞的时间，

为第一个粒子的初始粒子半径，

当得到两个即将发生碰撞的粒子发生碰撞的时间之后，可以根据即将发生碰撞的粒子发生碰撞的时间，预测两个即将发生碰撞的粒子的碰撞时刻。

需要说明的是，上述未来事件并不仅仅包括粒子间的碰撞事件，还可以包括其他事件，例如：粒子穿越网格边界的事件、中间产物的粒子半径恰好变化为生成物的粒子半径事件等。

进一步的，可以将预测的未来事件按照时间发生的先后顺序进行整理，以形成硬球模型区域的未来事件列表。

在步骤S50中，可以基于未来事件列表处理所述硬球模型区域内粒子间的碰撞和反应，并基于反应的时间步处理所述拟颗粒模型区域内粒子间的碰撞和反应。

具体的，可以根据初始未来事件列表里的事件发生时间的先后顺序，判断硬球模型区域中的两个即将发生碰撞的粒子是否发生反应。需要说明的是，当两个即将发生碰撞的粒子在质心连线方向的相对平动能之和大于设定的反应临界能量时，则认为该粒子能够发生分解反应。并且，无论是粒子与模拟体系的活性点位发生碰撞，还是粒子与模拟体系中的其他粒子发生碰撞时，都有可能发生分解反应。反之，当两个即将发生碰撞的粒子在质心连线方向的相对平动能之和小于或等于设定的反应临界能量时，则认为该粒子一般不会发生分解反应。

当硬球模型区域中的两个即将发生碰撞的粒子发生反应时，可以根据粒子的发生的反应，更新未来事件列表和硬球模型区域中反应物粒子、中间产物粒子和生成物粒子的种类和数量。

进一步的，在本公开的一个实施例中，针对拟颗粒模型区域。由于拟颗粒模型区域位于进程的边界处，并且在拟颗粒模型区域发生的反应受到时间驱动，因此，通常情况下需要分不同的情况处理拟颗粒模型区域中的分解反应过程，例如：需要分别考虑两个反应物粒子处在一个进程中和处在不同的进程中的情况，以及分别考虑产生的生成物粒子处在同一个进程中和处在不同进程中的情况。此外，产生的生成物粒子有可能都处在拟颗粒模型区域或者硬球颗粒模型区域，也有可能部分位于拟颗粒模型区域，部分位于硬球颗粒模型区域。因此，为了防止进程间通讯时有可能发生发送的粒子数与收回的粒子数不相等的问题，并且简化计算过程提高计算效率，本公开提供了一种新的处理方式。

可以对拟颗粒模型区域和硬球模型区域内的所有粒子均添加第一标记。在当前反应的时间步中，可以处理拟颗粒模型区域内各粒子间的碰撞，并即将分解的粒子的第一标记修改为第二标记；在当前反应的时间步后，并在下一个反应的时间步前，可以首先检测拟颗粒模型区域内各个粒子的标记，将添加有第二标记的粒子分解，以形成中间产物粒子。

举例而言：可以对拟颗粒模型区域和硬球模型区域中的每一个粒子存储的物理属性中均添加一个表示粒子是否需要分解的标记Decflag，并且设置粒子的第一标记为Decflag＝0。只有在粒子即将要分解的情况下将该粒子的第一标记修改为第二标记，该第二标记可以为Decflag＝1。但是需要说明的是，该粒子并不是在当前反应的时间步中发生分解。当前反应的时间步仅用于处理拟颗粒模型区域内各粒子的碰撞，在当前反应的时间步处理完碰撞后，该粒子被送回其相应的进程中，并将粒子的标记修改为第二标记，即Decflag＝1。当下一个反应的时间步开始之前，粒子进入拟颗粒模型区域后，首先需要检测粒子的Decflag。若该粒子的Decflag＝1，则将该粒子分解为中间产物粒子。可以理解的是：此处所述的下一个反应的时间步开始之前指的是：当下一个反应的时间步同步更新拟颗粒模型区域内各个粒子的状态之前。

另外，需要注意的是，当粒子分解为中间产物粒子时，需要更新该粒子和中间产物粒子周围处在中间网格区中的未来事件列表。产生的中间产物粒子是否在同一个进程并不会对具体的计算操作产生影响，这是因为当同步拟颗粒模型区域内所有粒子的状态后，下一步操作是将粒子发送到其高级进程中。即，拟颗粒模型区域并不会因为产生的中间产物粒子不在同一个进程而增加很多计算量。由于拟颗粒模型算法是时间驱动，粒子在当前的时间步分解与在下一个时间步开始之前进行分解在计算操作上并没有什么区别。由此可见，这样的处理方式既简化了计算，又可以得到合理的模拟结果。

进一步的，拟颗粒模型区域允许粒子之间发生重叠，在同步更新所有的粒子的状态时，需要根据每个粒子的粒子半径增长速度更新粒子的半径。只有当粒子半径增长速度不为0时，对粒子的半径进行更新；当粒子半径等于0时，不对粒子半径进行更新。需要注意的是，由于拟颗粒模型受时间驱动，从而有可能发生中间产物粒子的半径小于生成物的半径或者中间产物粒子的半径大于生成物的半径的情况。但是，无论遇到哪一种情况，均需要在当前时间步设置中间产物粒子的半径等于生成物粒子的半径，并修改粒子类型参数为生成物粒子的类型参数，以及设置粒子半径增长速度为0。

在本公开的一个实施例中，硬球模型区域中的各粒子中至少包括反应物粒子。可以根据反应物粒子发生的反应，将模拟体系中发生反应的反应物粒子转化为一个中间产物粒子，并生成另一个中间产物粒子。可以删除未来事件列表中关于反应物粒子的未来事件，并增加关于中间产物粒子的未来事件，以更新未来事件列表。

进一步的，在增加关于中间产物粒子的未来事件时，需要考虑中间产物粒子的初始速度和初始位置。中间产物粒子的初始速度可以依据反应物粒子碰撞后更新的反应物粒子的速度确定中间产物粒子的初始速度。基于反应前后动量守恒和能量守恒的原则，可以得到与生成物粒子的速度相关的标量方程，当生成物粒子具有两种时，标量方程的数量为4个。

需要说明的是，由于速度为矢量，从而当生成物粒子具有两种时，其速度可以表达为6个标量，此时，求解4个标量方程无法得到唯一的、确定的6个标量。因此，本公开选择了满足4个标量方程的一组简单解，即两个中间产物粒子的初速度与粒子碰撞后更新的反应物粒子的速度相等时的解。

依据动量守恒和能量守恒，可以获得任意两个粒子碰撞后的速度更新公式：

其中，V₁'为第一个粒子碰撞后的速度，V₂'为第二个粒子碰撞后的速度，e为碰撞恢复系数，m₁为第一个粒子的质量，m₂为第二个粒子的质量，V₁为第一个粒子碰撞前的速度，V₂为第二个粒子碰撞前的速度，r₁为第一个粒子碰撞时的位置，r₂为第二个粒子碰撞时的位置。

更进一步的，在确定中间产物粒子的初始位置时，为了节省搜索合适的粒子位置和检测粒子之间是否发生堆叠而产生的计算量，可以将中间产物粒子的质心放置在反应物粒子的某一条直径上。

举例而言，如图2所示，当生成物粒子具有两个时，可以将两个中间产物粒子的质心放置在反应物粒子的某一条直径的两端。其位置可以由下式进行计算：

其中，B'为第一个中间产物粒子，C'为第二个中间产物粒子，r_B'X为第一个中间产物粒子在X方向上的位置，r_B'Y为第一个中间产物粒子在Y方向上的位置，r_B'Z为第一个中间产物粒子在Z方向上的位置，r_X为中间产物粒子在X方向上的位置，r_Y为中间产物粒子在Y方向上的位置，r_Z为中间产物粒子在Z方向上的位置，d_AB'为第一个中间产物粒子的质心与反应物粒子的质心之间的距离，

为反应物粒子的直径方向与X方向之间的夹角，θ为反应物粒子的直径方向与Z方向之间的夹角，r_C'X为第二个中间产物粒子在X方向上的位置，r_C'Y为第二个中间产物粒子在Y方向上的位置，r_C'Z为第二个中间产物粒子在Y方向上的位置，d_AC'为第二个中间产物粒子的质心与反应物粒子的质心之间的距离。

需要说明的是，当生成物为其他数量时，也可以使用上述公式进行计算，只需要增加或者减少上述公式的数量即可，这均在本公开的保护范围之内。

在本公开的一个实施例中，当对中间产物粒子的初始位置确定了以后，可以将中间产物粒子插入到相应的网格中。需要注意的是，中间产物的质心位置与反应物粒子的质心位置不一定在同一个网格中。由于模拟体系中的粒子数量发生了变化，从而需要构建反应物粒子和中间产物粒子所在的网格和周围相邻的网格的未来事件列表。当发生反应的反应物粒子生成中间产物粒子后，与该发生反应的反应物粒子有关的未来事件均作为失效事件进行处理。

另外，当本公开采用计算机进行计算的时候，可以将存储反应物粒子的物理属性的内存块转换为存储中间产物粒子的物理属性的内存块，即：只将内存块中的反应物粒子的物理属性修改为中间产物粒子的物理属性，而不需要彻底删除反应物粒子。通过此方式，能够节省计算开销，并提高计算效率。

在本公开的一个实施例中，可以根据更新的未来事件列表，对中间产物粒子的碰撞进行预测和处理，并更新中间产物粒子的半径。具体的，可以引入粒子半径增长速率，并将反应物粒子的位置坐标标记为中间产物粒子的位置坐标，并且保持中间产物粒子的半径与反应物粒子的半径大小相等，以避免发生粒子的半径突然增大或者减小的情况。从而能够在计算的过程中，不需要检测生成物是否与周围现有的粒子发生重叠，大大节省了因为检索合适的生成物位置而产生的计算量。

因此，中间产物粒子的半径更新公式可以为：

其中，

为中间产物粒子在t₁时刻的粒子半径，

为中间产物粒子在t₀时刻的粒子半径，v_z为中间产物粒子的粒子增长速度，t₁为更新后的当前时刻，t₀为更新前的时刻。另外，当v_z>0时表示中间产物粒子的半径逐渐增大，当v_z<0时表示中间产物粒子的半径逐渐减小。

进一步的，当中间产物的粒子半径与生成物粒子半径相同时，中间产物粒子为新生成的生成物粒子。具体的，在模拟体系的演化过程中，硬球模型区域将根据中间产物粒子的半径和生成物粒子的半径之间的关系，相应的增大或者减小中间产物粒子半径，即：中间产物粒子在与其他粒子发生碰撞时，中间产物粒子的半径将会逐渐发生变化，直到中间产物粒子的半径与生成物粒子的半径相同时为止。同时，可以将该中间产物粒子的类型更新为生成物粒子的类型，以形成新生成的生成物粒子。但是，需要说明的是，当处理中间产物粒子的半径与生成物粒子的半径相等这一未来事件时，需要设置生成物的粒子增长速度为0。

更进一步的，可以根据转化的反应物粒子和新生成的生成物粒子更新硬球模型区域中的反应物粒子、中间产物粒子和生成物粒子的种类和数量。

在步骤S60中，可以统计模拟体系中的反应物粒子、中间产物粒子和生成物粒子信息，并生成统计结果，当模拟体系中的反应物粒子的数量为零时，该统计结果为该分解反应的分析结果。

具体地，可以设置第二时间步，并且每隔第二时间步对反应物粒子、中间产物粒子和生成物粒子信息进行统计。该第二时间步与第一时间步可以相同。另外，此处所说的反应物粒子和生成物粒子信息可以包括：反应物粒子、中间产物粒子和生成物粒子的种类、数量等。但不限于此，还可以实时的统计模拟体系中的反应物粒子、中间产物粒子和生成物粒子信息，这均在本公开的保护范围之内。

除此之外，当模拟体系中的反应物粒子的数量不为零时，可以重复上述步骤S50。

下面通过具体的实施例对上述各个步骤进行进一步的说明，并对上述方法的准确性进行验证。

该实施例中的模拟体系可以为沸石ZSM-5体系，并且选择T12位点为沸石ZSM-5的活性位点。当发生碰撞的两个粒子在质心连线方向的相对平动能之和大于设定的反应的临界能量时，该分解反应可以为：

反应R1：A→B+C

反应R2：A→2B，

其中A未反应物粒子，B和C为生成物粒子。

反应R1和R2的临界能量分别标记为E_al和E_a2。只有在反应R1中，反应物粒子A必须与沸石的T12位点发生碰撞后才有可能发生反应。反应R2中，反应物粒子A是在催化剂的表面和孔道中发生反应，反应本身与沸石的活性位点没有关系。因此，设定反应R1为主反应，反应R2可被认为是副反应。模拟体系在X、Y和Z三个方向采用周期性边界条件，即：当某粒子从某个方向离开模拟体系时，将以相同的粒子类型从对称的位置进入模拟体系。

为便于讨论分析，在该实施例中使用反应物粒子A的质量m_A＝2.66×10^-26kg、反应物粒子A的粒子半径为：

特征时间为：

进一步的，可以对反应物粒子A的质量、反应物粒子A的粒子半径及特征时间的单位分别进行无量纲化处理，并且将活化能E_al无量纲化处理为1.0。

在第一个实施例中，反应R1和R2的临界能量分别为E_a1＝1.00和E_a2＝0.85。模拟体系中初始的反应物粒子为A，反应后生成物粒子为B和C。模拟体系初始化时，反应物粒子A随机地分布在沸石Y方向的直孔道中。设置模拟体系的气体负载量为平均每个沸石晶胞单元中含有4个气体分子。模拟体系含有的沸石晶胞数量为8×8×8(16.056×15.936×10.736nm²)晶胞。模拟体系的反应温度设置为398K。不同粒子的质量m和粒子半径σ_HS如表1所示。基于反应前后质量守恒设定粒子的质量，并且粒子的半径根据粒子的质量的变化而等比例发生变化，即：粒子的质量越大，其半径也越大。

在模拟体系初始化时，各粒子的数量如表1中的n₀所示。待完成模拟体系初始化后，发起程序开始计算，当检测到体系中粒子A被消耗完毕时，即认为完成了所有反应，此时模拟体系演化结束，获得的反应后的粒子数量如表1中的n₁所示。模拟结果表明，反应前后粒子数量的变化符合预设的反应关系，验证了本公开提供的分析方法的正确性。

表1

名称	m	σ<sub>HS</sub>	n<sub>0</sub>	n<sub>1</sub>
					A	1.00	1.00	2048	0
B	0.50	0.50	0	2438
					C	0.50	0.50	0	1658

在第二个实施例中，反应R1和R2的临界能量分别为E_a1＝1.00和E_a2＝1.10。模拟体系中初始的反应物粒子为A，反应后生成物粒子为B和C。模拟体系初始化时，反应物粒子A随机地分布在沸石Y方向的直孔道中。设置模拟体系的气体负载量为平均每个沸石晶胞单元中含有3个气体分子。模拟体系含有的沸石晶胞数量为9×9×9(18.063×17.928×12.078nm²)晶胞。模拟体系的反应温度设置为448K。不同粒子的质量m和粒子半径σ_HS如表2所示。基于反应前后质量守恒设定粒子的质量，并且粒子的半径根据粒子的质量的变化而等比例发生变化，即：粒子的质量越大，其半径也越大。

在模拟体系初始化时，各粒子的数量如表2中的n₀所示。待完成模拟体系初始化后，发起程序开始计算，当检测到体系中粒子A被消耗完毕时，即认为完成了所有反应，此时模拟体系演化结束，获得的反应后的粒子数量如表2中的n₁所示。模拟结果表明，反应前后粒子数量的变化符合预设的反应关系，验证了本公开提供的分析方法的正确性。

表2

名称	m	σ<sub>HS</sub>	n<sub>0</sub>	n<sub>1</sub>
					A	1.00	1.00	2187	0
B	0.50	0.50	0	2672
					C	0.50	0.50	0	1702

在第三个实施例中，反应R1和R2的临界能量分别为E_a1＝1.00和E_a2＝1.20。模拟体系中初始的反应物粒子为A和少量的生成物粒子B，反应后生成物粒子为B和C。模拟体系初始化时，反应物粒子A和少量的生成物粒子B随机地分布在沸石Y方向的直孔道中。设置模拟体系的气体负载量为平均每个沸石晶胞单元中含有5个气体分子。模拟体系含有的沸石晶胞数量为11×11×11(22.077×21.912×14.762nm²)晶胞。模拟体系的反应温度设置为418K。不同粒子的质量m和粒子半径σ_HS如表3所示。基于反应前后质量守恒设定粒子的质量，并且粒子的半径根据粒子的质量的变化而等比例发生变化，即：粒子的质量越大，其半径也越大。

在模拟体系初始化时，各粒子的数量如表3中的n₀所示。待完成模拟体系初始化后，发起程序开始计算，当检测到体系中粒子A被消耗完毕时，即认为完成了所有反应，此时模拟体系演化结束，获得的反应后的粒子数量如表3中的n₁所示。模拟结果表明，反应前后粒子数量的变化符合预设的反应关系，验证了本公开提供的分析方法的正确性。

表3

名称	m	σ<sub>HS</sub>	n<sub>0</sub>	n<sub>1</sub>
					A	1.00	1.00	6355	0
B	0.50	0.50	300	7057
					C	0.50	0.50	0	5953

在第四个实施例中，反应R1和R2的临界能量分别为E_a1＝1.00和Ea2＝1.00。模拟体系中初始的反应物粒子为A和少量的生成物粒子B，反应后生成物粒子为B和C。模拟体系初始化时，反应物粒子A和少量的生成物粒子B随机地分布在沸石Y方向的直孔道中。设置模拟体系的气体负载量为平均每个沸石晶胞单元中含有6个气体分子。模拟体系含有的沸石晶胞数量为10×10×10(20.070×19.920×13.420nm²)晶胞。模拟体系的反应温度设置为458K。不同粒子的质量m和粒子半径σ_HS如表4所示。基于反应前后质量守恒设定粒子的质量，并且粒子的半径根据粒子的质量的变化而等比例发生变化，即：粒子的质量越大，其半径也越大。

在模拟体系初始化时，各粒子的数量如表4中的n₀所示。待完成模拟体系初始化后，发起程序开始计算，当检测到体系中粒子A被消耗完毕时，即认为完成了所有反应，此时模拟体系演化结束，获得的反应后的粒子数量如表4中的n₁所示。模拟结果表明，反应前后粒子数量的变化符合预设的反应关系，验证了本公开提供的分析方法的正确性。

表4

名称	m	σ<sub>HS</sub>	n<sub>0</sub>	n<sub>1</sub>
					A	1.00	1.00	5650	0
B	0.50	0.50	350	6513
					C	0.50	0.50	0	5137

因此，综上所述，通过不同的反应温度、模拟体系的尺寸、反应活化能和初始粒子种类和数量的模拟，所进行分析后的结果均满足要求，从而充分验证了本公开所提供的分析方法的准确性。

本领域技术人员在考虑说明书及实践这里公开的发明后，将容易想到本公开的其它实施方案。本公开旨在涵盖本公开的任何变型、用途或者适应性变化，这些变型、用途或者适应性变化遵循本公开的一般性原理并包括本公开未公开的本技术领域中的公知常识或惯用技术手段。说明书和实施例仅被视为示例性的，本公开的真正范围和精神由所附的权利要求指出。