化合反应类型的硬球-拟颗粒耦合模拟方法与流程

化合反应类型的硬球
‑
拟颗粒耦合模拟方法
技术领域
1.本公开涉及经典分子动力学模拟计算技术领域，尤其涉及一种化合反应类型的硬球
‑
拟颗粒耦合模拟方法。


背景技术：

2.研究纳微尺度上的扩散和反应过程耦合机制对于理解宏观、复杂的化学反应体系具有重要意义，例如单分子催化反应过程。处在纳微尺度上的气体分子，其状态可以很好地符合稀薄气体假设，并且表现出了分子级别的离散特性，不再适用于连续介质假设，必须寻求其他的模拟计算方法。针对该研究领域，分子动力学(molecular dynamics，md)方法因其较高的计算精度表现出了巨大的应用潜力。md方法是联系微观世界和宏观世界的一种强有力的计算机模拟方法，能够对研究体系进行分子水平的描述，通过求解牛顿运动方程获得各原子、分子的微观状态随时间的动态变化规律，从而获得分子体系统计的各种物理和化学特性。
3.基于粒子间的相互作用力，经典的md模拟方法包括硬球(hard
‑
sphere,hs)模型、软球(soft
‑
sphere)模型，以及为了解决hs模型难以并行的问题而提出的拟颗粒模型(pseudo
‑
particle modeling,ppm)。hs粒子的碰撞可以等效为一系列ss粒子的相互作用，但是ss模型的计算效率并不如hs模型高。ppm结合了hs模型较高的计算效率和ss模型较好的可扩展性的优点。但是，针对稀薄气体的模拟，由于所有粒子的状态在ppm体系中需要同步更新，使得ppm的计算效率明显低于hs模型。hs
‑
ppm耦合算法能够同时满足较高的计算精度、较大的计算规模(纳米和纳秒的时空尺度)、及较高的计算效率的大规模并行计算需求。
4.然而，目前的文献资料还没有关于在hs
‑
ppm耦合算法中如何实现简单典型反应动力学模型方面的报道，制约了hs
‑
ppm耦合算法在研究纳微尺度上扩散和反应过程耦合机制领域的应用。
5.所述背景技术部分公开的上述信息仅用于加强对本公开的背景的理解，因此它可以包括不构成对本领域普通技术人员已知的现有技术的信息。


技术实现要素：

6.本公开的目的在于提供一种化合反应类型的硬球
‑
拟颗粒耦合模拟方法，该方法具有可并行计算、计算效率高、计算规模大、适用范围广、能够获得全系统演化的信息等优点。
7.为实现上述发明目的，本公开采用如下技术方案：
8.根据本公开的第一个方面，提供一种化合反应类型的硬球
‑
拟颗粒耦合模拟方法，包括：
9.获取粒子a、粒子b和粒子c各自的粒子属性参数，并设置粒子c’的粒子属性参数，a+b
→
c，c’为化合反应过程中的中间产物，所述粒子属性参数包括粒子质量、粒子硬球半径、粒子类型编号和粒子半径增长速率；
10.设置化合反应模型的模型参数，以及所述粒子a、粒子b、粒子c’和粒子c的粒子数量和反应参数，所述模型参数包括模型反应温度、进程数和模型几何参数；
11.基于化合反应模型的所述模型参数，以及所述粒子a和粒子b各自的粒子属性参数，确定所述粒子a和粒子b各自的初始速度和初始位置；
12.根据所述模型参数，对化合反应模型进行进程划分、网格划分，以及硬球模拟子区和拟颗粒模拟子区划分，并对所述硬球模拟子区建立未来事件列表；
13.开始模拟计算过程，检测所述粒子a、粒子b、粒子c和粒子c’各自的反应状态、粒子半径和粒子位置的变化过程，并统计各个粒子的数量信息，当所述粒子a或所述粒子b的数量为0时，输出计算结果。
14.在本公开的一种示例性实施例中，所述获取粒子a、粒子b和粒子c各自的粒子属性参数，并设置粒子c’的粒子属性参数，a+b
→
c，c’为化合反应过程中的中间产物，所述粒子属性参数包括粒子质量、粒子硬球半径、粒子类型编号和粒子半径增长速率之后还包括：
15.对所述粒子属性参数进行无量纲化处理。
16.在本公开的一种示例性实施例中，所述基于化合反应模型的模型参数，以及所述粒子a和粒子b各自的粒子属性参数，获取粒子a和粒子b各自的初始速度和初始位置包括：
17.根据化合反应模型的所述模型反应温度，以及所述粒子a和所述粒子b各自的粒子质量，获取所述粒子a和所述粒子b各自的初始速度；
18.根据化合反应模型的模型几何参数，获取所述粒子a和所述粒子b各自的初始位置。
19.在本公开的一种示例性实施例中，所述硬球模拟子区的数量为多个，多个所述硬球模拟子区之间通过所述拟颗粒模拟子区连接。
20.在本公开的一种示例性实施例中，所述根据所述模型参数，对化合反应模型进行进程划分、网格划分，以及硬球模拟子区和拟颗粒模拟子区划分，并对所述硬球模拟子区建立未来事件列表包括：
21.基于所述进程数将化合反应模型划分为多个子模块，每个进程负责一个所述子模块的模拟计算；
22.对每个所述子模块进行网格划分，使得每个所述子模块的边界区域为所述拟颗粒子区，剩下区域为所述硬球模拟子区，所述拟颗粒子区围绕所述硬球模拟子区；
23.对所述硬球模拟子区建立未来事件列表。
24.在本公开的一种示例性实施例中，所述开始模拟计算过程，检测所述粒子a、粒子b、粒子c和粒子c’各自的反应状态、粒子半径和粒子位置的变化过程，并统计各个粒子的数量信息，当所述粒子a或所述粒子b的数量为0时，输出计算结果包括：
25.检测是否有粒子发生碰撞，当检测到有两个粒子发生碰撞时，判断发生碰撞的粒子的种类；
26.当发生碰撞的粒子为所述粒子a和所述粒子b时，根据预设的反应原则，判断发生碰撞的所述粒子a和所述粒子b是否发生反应，并根据判断结果做出对应的处理，所述反应原则根据所述反应参数进行设定；当发生碰撞的粒子不是所述粒子a和所述粒子b时，则更新粒子碰撞后各自对应的速度、粒子半径和粒子位置；
27.统计各类粒子的数量信息，当所述粒子a或所述粒子b的数量为0时，输出计算结
果。
28.在本公开的一种示例性实施例中，所述判断发生碰撞的所述粒子a和所述粒子b是否发生反应，并根据判断结果做出对应的处理包括：
29.若判断所述粒子a和所述粒子b发生反应，则删除发生碰撞的所述粒子a和所述粒子b中的其中一者，并将另一者转存为所述粒子c’，且更新所述粒子a、所述粒子b和/或所述粒子c’在碰撞后各自对应的速度、粒子半径和粒子位置，以及更新所述硬球模拟子区的未来事件列表；若判断所述粒子a和所述粒子b未发生反应，则更新所述粒子a和所述粒子b碰撞后各自对应的速度、粒子半径和粒子位置。
30.在本公开的一种示例性实施例中，当所述粒子a和所述粒子b发生反应，所述粒子a、所述粒子b和所述粒子c’满足第一关系式和第二关系式
31.其中，e为结余能量，e＞0，m
a
和m
b
分别为粒子a和粒子b在碰撞前的粒子质量，v
a
、v
b
分别为粒子a和粒子b在碰撞前的速度，mc’为发生反应后粒子c’的粒子质量，v
c
为发生反应后粒子c’的速度。
32.在本公开的一种示例性实施例中，当所述粒子a和所述粒子b发生反应，比较所述粒子a和所述粒子b的粒子半径的大小，删除其中粒子半径较大的一者，并将其中粒子半径较小的一者转存为粒子c’。
33.在本公开的一种示例性实施例中，在所述拟颗粒子区，当所述粒子a和所述粒子b发生反应，且所述粒子a和所述粒子b不在同一个进程时，在当前时间步内比较所述粒子a和粒子b的粒子半径的大小，并对其中粒子半径较大的粒子予以标记，在下一个时间步同步演化所述拟颗粒子区所有粒子的状态之前，检测所有粒子的标记信息，删除被予以标记的粒子。
34.在本公开的一种示例性实施例中，所述粒子c’的粒子半径变化满足第三关系式，
[0035][0036]
其中，t0和t1分别表示粒子运动的上一个时刻及当前时刻，和分别对应于粒子在t0和t1时刻的粒子半径，v
r
为所述粒子c’的粒子半径增长速率，v
r
>0和v
r
<0分别表示粒子的粒子半径变化过程是一个逐渐增大和逐渐减小的过程。
[0037]
在本公开的一种示例性实施例中，所述判断发生碰撞的所述粒子a和所述粒子b是否发生反应，并根据判断结果做出对应的处理还包括：
[0038]
监测粒子c’的运动状态和粒子半径，当所述粒子c’的粒子半径等于所述粒子c的硬球半径时，将粒子c’的粒子类型编号更改为所述粒子c对应的粒子类型编号。
[0039]
在本公开的一种示例性实施例中，更新所述硬球模拟子区的未来事件列表包括：
[0040]
删除与所述粒子a和粒子b有关的无效事件，增加与粒子c’有关的事件。
[0041]
本公开提供的化合反应类型的硬球
‑
拟颗粒耦合模拟方法，采用了硬球
‑
拟颗粒耦合算法，该算法忽略了粒子之间的相互吸引力，将粒子假设为无内部结构的硬球，硬球分子之间除碰撞瞬间外不存在相互作用力，因此，无需计算粒子的动态电荷等信息，可以在很大
程度上减少整个反应体系演化时的计算量。本公开将化合反应模型划分为硬球模拟子区和拟颗粒模拟子区，其中，硬球模拟子区采用事件驱动，拟颗粒子区采用时间驱动，因此，只需同步演化拟颗粒子区内粒子，而不必同步演化整个反应体系内所有粒子，大大节省了计算时间，提高了计算效率。此外，本公开提供的硬球
‑
拟颗粒耦合模拟方法能够针对反应过程符合a+b
→
c的化合反应体系进行模拟研究，只需要根据具体的反应物质相应地更新粒子的粒子属性参数和化合反应模型的模型参数，具有较强的通用性。
附图说明
[0042]
通过参照附图详细描述其示例实施方式，本公开的上述和其它特征及优点将变得更加明显。
[0043]
图1是本公开示例性实施例中硬球
‑
拟颗粒耦合模拟方法流程示意图。
具体实施方式
[0044]
现在将参考附图更全面地描述示例实施例。然而，示例实施例能够以多种形式实施，且不应被理解为限于在此阐述的范例；相反，提供这些实施例使得本公开将更加全面和完整，并将示例实施例的构思全面地传达给本领域的技术人员。所描述的特征、结构或特性可以以任何合适的方式结合在一个或更多实施例中。在下面的描述中，提供许多具体细节从而给出对本公开的实施例的充分理解。
[0045]
在图中，为了清晰，可能夸大了区域和层的厚度。在图中相同的附图标记表示相同或类似的结构，因而将省略它们的详细描述。
[0046]
所描述的特征、结构或特性可以以任何合适的方式结合在一个或更多实施例中。在下面的描述中，提供许多具体细节从而给出对本公开的实施例的充分理解。然而，本领域技术人员将意识到，可以实践本公开的技术方案而没有所述特定细节中的一个或更多，或者可以采用其它的方法、组元、材料等。在其它情况下，不详细示出或描述公知结构、材料或者操作以避免模糊本公开的主要技术创意。
[0047]
当某结构在其它结构“上”时，有可能是指某结构一体形成于其它结构上，或指某结构“直接”设置在其它结构上，或指某结构通过另一结构“间接”设置在其它结构上。
[0048]
用语“一个”、“一”、“所述”用以表示存在一个或多个要素/组成部分/等；用语“包括”和“具有”用以表示开放式的包括在内的意思并且是指除了列出的要素/组成部分/等之外还可存在另外的要素/组成部分/等。用语“第一”和“第二”等仅作为标记使用，不是对其对象的数量限制。
[0049]
本公开提供的化合反应类型的硬球
‑
拟颗粒耦合模拟方法，适用于研究纳微尺度上发生的气固微观化学反应过程，能够为模拟研究纳微孔道内发生的扩散和反应耦合过程提供技术支撑，对纳微尺度上扩散和反应耦合机制的研究具有重要意义。
[0050]
针对纳微尺度上发生的气固微观化学反应过程，由于纳微尺度上气体分子之间的距离相距较远，气体分子间的引力作用可被忽略，基于气体碰撞理论，本公开将气体分子假设为无内部结构的硬球且硬球间除碰撞外无相互作用力的假设比较合理。另外，针对扩散和反应过程的耦合机制的研究，模拟研究关注的重点不是反应过程中断键、成键这样具体的变化过程，可以采用经典的md方法进行传质和反应过程的研究，因此，本公开采用hs
‑
ppm
耦合算法模拟研究纳微尺度上的扩散和反应过程的耦合机制。
[0051]
如图1所示，本公开提供一种化合反应类型的硬球
‑
拟颗粒耦合模拟方法，适用于模拟研究包含由多个化合反应组成的化学反应网络的复杂化学反应体系能够进行高效和大规模的并行模拟计算，具有可并行计算、计算效率高、计算规模大、适用范围广、能够获得全反应体系演化的信息等优点。该方法包括以下步骤：
[0052]
步骤s100，获取粒子a、粒子b和粒子c各自的粒子属性参数，并设置粒子c’的粒子属性参数，a+b
→
c，c’为化合反应过程中的中间产物，粒子属性参数包括粒子质量、粒子硬球半径、粒子类型编号和粒子半径增长速率；
[0053]
步骤s200，设置化合反应模型的模型参数，以及粒子a、粒子b、粒子c’和粒子c的粒子数量和反应参数，模型参数包括模型反应温度、进程数和模型几何参数；
[0054]
步骤s300，基于化合反应模型的模型参数，以及粒子a和粒子b各自的粒子属性参数，确定粒子a和粒子b各自的初始速度和初始位置；
[0055]
步骤s400，根据模型参数，对化合反应模型进行进程划分、网格划分，以及硬球模拟子区和拟颗粒模拟子区划分，并对硬球模拟子区建立未来事件列表；
[0056]
步骤s500，开始模拟计算过程，检测粒子a、粒子b、粒子c和粒子c’各自的反应状态、粒子半径和粒子位置的变化过程，并统计各个粒子的数量信息，当粒子a或粒子b的数量为0时，输出计算结果。
[0057]
本公开提供的化合反应类型的硬球
‑
拟颗粒耦合模拟方法，采用了硬球
‑
拟颗粒耦合算法，该算法忽略了粒子之间的相互吸引力，将粒子假设为无内部结构的硬球，硬球分子之间除碰撞瞬间外不存在相互作用力，因此，无需计算粒子的动态电荷等信息，可以在很大程度上减少整个反应体系演化时的计算量。本公开将化合反应模型划分为硬球模拟子区和拟颗粒模拟子区，其中，硬球模拟子区采用事件驱动，拟颗粒子区采用时间驱动，因此，只需同步演化拟颗粒子区内粒子，而不必同步演化整个反应体系内所有粒子，大大节省了计算时间，提高了计算效率。此外，本公开提供的硬球
‑
拟颗粒耦合模拟方法能够针对反应过程符合a+b
→
c的化合反应体系进行模拟研究，只需要根据具体的反应物质相应地更新粒子的粒子属性参数和化合反应模型的模型参数，具有较强的通用性。
[0058]
下面结合附图对本公开实施方式提供的一种化合反应类型的硬球
‑
拟颗粒耦合模拟方法的各步骤进行详细说明。
[0059]
步骤s100，获取粒子a、粒子b和粒子c各自的粒子属性参数，并设置粒子c’的粒子属性参数，a+b
→
c，c’为化合反应过程中的中间产物，粒子属性参数包括粒子质量、粒子硬球半径、粒子类型编号和粒子半径增长速率。
[0060]
在步骤s100中，在化合反应中，反应物a与反应物b进行反应，最终生成生成物c。c’为化合反应过程中的中间产物。在本公开中，将c’定义为存储生成物。
[0061]
粒子属性参数为与化合反应中粒子有关的参数，包括粒子质量、粒子硬球半径、粒子类型编号、粒子半径增长速率等。其中，粒子类型编号是指针对不同的粒子对应给予不同的数字标签。如粒子a的粒子类型编号为1，粒子b的粒子类型编号为2，粒子c的粒子类型编号为3，粒子c’的粒子类型编号为4。粒子半径增长速率用于描述各个粒子的粒子半径变化过程，其中粒子a和粒子b的粒子增长速率为0，粒子c’的粒子增长速率根据实际情况进行设定。
[0062]
在本公开一些实施例中，步骤s100包括：
[0063]
步骤s110，根据粒子c的粒子质量和粒子硬球半径设置粒子c’的粒子质量和粒子半径增长速率。
[0064]
粒子c’为存储生成物，在该步骤中，粒子c’的粒子质量根据粒子c的质量进行设定，具体地，粒子c’的粒子质量等于粒子c的质量，粒子c’的粒子半径增长速率与粒子c的粒子硬球半径有关。粒子a、粒子b和粒子c各自的粒子质量和粒子硬球半径，可根据实验经验自由设定，也可根据实验数据获得。本公开将粒子质量和粒子硬球半径作为粒子属性参数，研究化合反应过程中的粒子动力学性质改变。在此需说明的是，在后续模拟过程中，随着时间的推移，只有粒子c’的粒子半径发生改变。
[0065]
在一些实施例中，步骤s100之前还包括：
[0066]
步骤s100
‑
0，建立化合反应的反应模型。
[0067]
在步骤s100
‑
0中，建立化合反应的反应模型，确定化合反应a+b
→
c中粒子的运动模拟空间。化合反应对应的反应模型可以包含单个模拟空间，也可以包含多个相同或不同的模拟空间。举例而言，sapo
‑
34是具有类菱沸石型结构的分子筛，其骨架主要由硅铝磷酸盐组成，具有双六元环、八元环与四元环构成的椭球形cha笼以及三维交叉孔道结构。sapo
‑
34分子筛既具有特殊的吸水性能和质子酸性，同时又有良好的水热稳定性，使得它可用作吸附剂和催化剂载体。可用于小分子的择形催化反应，比如作为甲醇制烯烃的催化剂，还可以作为汽车尾气净化催化剂载体。因此，针对沸石分子筛sapo
‑
34体系，可模拟sapo
‑
34分子筛的孔道结构，建立包含独特孔道构型的反应模型，以模拟研究化合反应中粒子在纳微孔道内的扩散、吸附等行为。在此需说明的是，该步骤可以是基于真实的孔道结构建立的反应模型，也可以是简化的反应模型，如直接建立直方孔道模型，具体本公开不做限定。
[0068]
在一些实施例中，步骤s100之后还包括：
[0069]
步骤s100
‑
1，对粒子属性参数进行无量纲化处理。
[0070]
在该步骤中，对粒子属性参数进行无量钢化处理，将粒子属性参数中的单位去除。在本公开一具体实施例中，采用粒子a的粒子质量和粒子硬球半径对其他粒子的粒子质量和粒子硬球半径进行无量钢化处理，以便于后续进行分析和讨论。此外，当粒子的属性参数还包括时间、速度等参数时，也对时间、速度等进行无量纲化处理。具体在一实施例中，采用粒子a的粒子质量m
a
、粒子硬球半径及特征时间分别对模拟体系中的质量、长度及时间的单位进行无量纲化处理。
[0071]
在另一些实施例中，也可采用粒子b或粒子c的粒子质量和粒子硬球半径对其他粒子的粒子质量和粒子硬球半径进行无量钢化处理，对此，本公开不做限定。
[0072]
在步骤s200中，设置化合反应模型的模型参数，以及粒子a、粒子b、粒子c’和粒子c的粒子数量和反应参数，模型参数包括模型反应温度、进程数和模型几何参数。
[0073]
在该步骤中，模型参数是指与化合反应模型整体有关的参数，例如模型的反应温度、进程数和模型几何参数。其中，模型反应温度也即反应体系温度，在整个模拟过程中保持不变。进程数可基于计算规模人为进行设定。模型几何参数包括反应模型整体的几何结构和几何尺寸，如反应模型的结构为直方孔道构型，该直方孔道构型在x、y、z方向上的长度等。在此需说明的是，模型几何参数也可以包括模型的体积、表面积等几何数据，具体可根据实际模拟需求进行设定，对此，本公开不做限定。
[0074]
粒子数量是指各个粒子的个数。在该步骤中，粒子数量可根据实验经验自由设定，也可根据实验数据获得。其中，粒子a和粒子b的粒子数量不为0，粒子c’的数量为0，粒子c的粒子数量可以为0，也可以不为0。在此需说明的是，该步骤中设置的粒子数量为模拟反应开始前的粒子数量。在后续模拟过程中，随着时间的推移，化合反应逐步进行，各个粒子的数量也在逐渐改变。
[0075]
反应参数是指与化合反应过程有关的参数，例如反应临界能。本公开引入反应临界能，作为反应物a和b发生反应的条件和规则。在后续模拟计算过程中，当反应物a和b发生碰撞时，基于反应临界能判断a和b是否发生了反应。
[0076]
在步骤s300中，基于化合反应模型的模型参数，以及粒子a和粒子b各自的粒子属性参数，确定粒子a和粒子b各自的初始速度和初始位置。
[0077]
在化合反应中，反应物a和反应物b反应生成生成物c。在该步骤中，根据化合反应模型的模型参数，以及粒子a和粒子b的粒子属性参数，确定反应物粒子a和粒子b的初始速度和初始位置。
[0078]
在一些实施例中，步骤s300包括：
[0079]
步骤s310，基于化合反应模型的模型反应温度，以及粒子a和粒子b各自的粒子质量，确定粒子a和粒子b各自的初始速度；
[0080]
步骤s320，基于化合反应模型的模型几何参数，确定粒子a和粒子b各自的初始位置。
[0081]
在步骤s310中，各个粒子的初始速度包括热速度和流动速度两部分，粒子的热速度v0可由模型反应温度t计算：(3d体系)、(3d体系)、(2d体系)，其中k
b
表示玻尔兹曼常数，m为对应粒子的粒子质量。物理量v
rx
、v
ry
和v
rz
分别表示粒子在x、y和z三个方向的流动速度。基于各向同性的球形分布原则，由计算式：
[0082][0083][0084][0085]
获得粒子在x、y和z三个方向的初始速度v
x
、v
y
和v
z
。速度方向与正x轴间的方位角及速度方向与正z轴间的仰角θ由时间种子随机生成，其变化范围分别为[0,2π)和[0,π)。
[0086]
基于上述内容，可获得粒子a和粒子b的初始速度。
[0087]
在步骤s320中，基于化合反应模型的模型几何参数，确定粒子a和粒子b各自的初始位置。
[0088]
基于反应模型的几何参数，为粒子随机分配初始位置。模型的几何参数包括反应模型整体的几何结构和几何尺寸。几何结构可以是各种三维几何结构，如长方体、圆锥体、球体或正方体等，也可以各种二维几何图形，如长方形、正多边形、圆形、不规则平面图形等，具体根据实际的模型需求进行设定。不同反应模型可设置不同的几何尺寸，如反应模型的结构为直方孔道构型，该直方孔道构型的几何尺寸可包括其在x、y、z方向上的长度。
[0089]
在一些实施例中，将各个粒子均匀地分布在反应模型的模拟空间内，并且要求各个粒子间不能发生重叠。
[0090]
在步骤s400中，根据模型参数，对化合反应模型进行进程划分、网格划分，以及硬球模拟子区和拟颗粒模拟子区划分，并对硬球模拟子区建立未来事件列表。
[0091]
在该步骤中，将整个化合反应模型依据应用的计算模型进行划分，以保证耦合算法较高的计算效率和计算精度及较好的可扩展性，具体地，将化合反应模型划分为硬球模拟子区和拟颗粒模拟子区，硬球模拟子区简称为hs子区，拟颗粒模拟子区简称为pp子区。hs子区采用hs模拟，pp子区采用ppm。
[0092]
相关技术中，hs模拟的核心是依次处理以发生时间排序的事件列表，此过程具有本质的串行性，因此，单独的hs模拟难以并行化。ppm是一种具有强扩展性的时间驱动算法，模型中允许拟颗粒间的轻微重叠。本公开采用hs
‑
ppm方法，通过pp子区将多个hs子区相连接，从而实现耦合模拟。在一些实施例中，硬球模拟子区的数量为多个，多个硬球模拟子区之间通过拟颗粒模拟子区连接。
[0093]
在本公开一些实施例中，步骤s400包括：
[0094]
步骤s410，基于进程数将反应模型划分为多个子模块，每个进程负责一个子模块的模拟计算；
[0095]
步骤s420，对每个子模块进行网格划分，使得每个子模块的边界区域为拟颗粒子区，剩余区域为硬球模拟子区，拟颗粒子区围绕硬球模拟子区；
[0096]
步骤s430，对硬球模拟子区建立未来事件列表。
[0097]
在步骤s410中，基于进程数将反应模型划分为多个子模块，每个进程负责一个子模块的模拟计算。
[0098]
采用的hs
‑
ppm耦合算法通过时间驱动的ppm模拟提高了事件驱动的hs模拟的可扩展性，使得能够基于计算规模设置所需要的进程数，然后根据进程数将整个反应模型划分为多个子模块，每个进程只负责一个子模块的模拟计算，相邻进程间通过消息传递接口(message passing interface,mpi)交换边界信息，进而实现大规模、多进程的并行模拟计算，计算规模可以达到10000数量级。
[0099]
在步骤s420中，对每个子模块进行网格划分，使得每个子模块的中部网格为硬球模拟子区，围绕硬球模拟子区网格的为拟颗粒子区。
[0100]
网格单元是算法中粒子存储和事件列表建立的基本单位。具体在一实施例中，可将靠近子模块边界的两层网格定为pp子区，采用ppm模拟，将pp子区围绕的内部网格定为hs子区，采用hs模拟，即每一个子模块中hs子区占据大部分区域，pp子区只占据边界的两层网格。
[0101]
对于每个进程对应的子模块，在子模块边界外侧构建一层虚拟网格，以此来处理该进程与相邻进程间的数据交换。
[0102]
在步骤s430中，对硬球模拟子区建立未来事件列表。
[0103]
未来事件可包含未来可能发生的所有事件，如各个粒子在网格中的位置变化、粒子的粒子半径变化、粒子间的相互碰撞等，具体本公开不做限定。
[0104]
在一些实施例中，硬球模拟子区建立未来事件列表前，还需对粒子未来发生事件进行预测，例如，预测粒子发生碰撞的时刻。在实际模拟中，预测粒子碰撞时刻对构建硬球模拟子区建的未来事件列表非常重要。
[0105]
某些粒子在整个模拟过程中，其粒子半径会随着反应的进行逐渐发生变化。本公
开在预测粒子的碰撞时刻时，引入粒子半径增长速率v
r
，v
r
>0和v
r
<0分别表示粒子的粒子半径变化过程是一个逐渐增大和逐渐减小的过程。若某个粒子的v
r
≠0，如粒子c’，则需要考虑粒子的σ
r
从当前运动时刻到预测碰撞的时刻这段时间间隔内发生的变化。因此，在预测粒子碰撞时刻的公式中引入粒子的v
r
。
[0106]
具体地，将当前时刻粒子的位置和半径分别标记为r0和碰撞时刻粒子的位置和半径分别标记为r1和由于hs模型不允许粒子之间发生重叠，首先需要满足的碰撞条件是两个粒子质心的距离等于两个粒子的是两个粒子质心的距离等于两个粒子的之和。即，r1和满足关系式：
[0107][0108]
其中，下标1和2表示碰撞的两个粒子。由粒子的速度v和r可以计算得到r1，即：
[0109]
r1＝r0+v
×
δt
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(5)
[0110]
其中，δt表示从当前运动时刻到预测碰撞时刻的时间间隔。由粒子的v
r
和σ
r
可以计算得到即：
[0111][0112]
其次，需要满足的碰撞条件是两个粒子必须有相互碰撞的趋势，即要求满足公式(4)有解的条件。将公式(5)和(6)代入公式(4)后，根据求根公式，需要满足：
[0113][0114]
的条件。其中，v
r1
和v
r2
分别表示粒子1和粒子2的v
r
。由此得到方程的两个解中较小的δt解为：
[0115][0116]
其中，σ
r1
和σ
r2
分别表示粒子1和粒子2的σ
r
。当预测体系中的任何粒子之间的碰撞时，都可以通过公式(7)判断两个粒子是否发生碰撞。若粒子对的位置和速度满足公式(7)，再通过公式(8)求解δt，获得预测的粒子碰撞时刻。
[0117]
在预测获得粒子的碰撞时刻后，可据此建立未来事件列表或在后续步骤中更新未来事件列表。
[0118]
在一具体实施例中，hs模型将体系内粒子参与的所有运动标记为不同的“事件”进行处理，标记的事件类型可以是真实发生的碰撞事件，也可以是为了便于粒子搜索而引入的虚拟事件。hs模型以网格为基本单位，存放在该网格单元内的所有粒子有可能发生的所有事件组成未来事件列表。首先需要对模拟体系内所有可能发生的事件类型(type)进行编号，包括普通粒子间的碰撞事件(type＝0)、粒子穿越虚拟网格边界事件(其中，穿越+x、+y、+z、
‑
x、
‑
y及
‑
z方向分别标记为type＝1、type＝2、type＝3、type＝4、type＝5、type＝6的事件类型编号)、粒子半径变化的粒子间的碰撞事件(type＝7)、反应发生事件(type＝8)、粒子半径发生变化的事件(type＝9)、粒子半径恰好变化到设定半径的事件(type＝10)。然后，预测每一个事件将要发生的时间。每一个事件作为完全二叉树(complete binary tree,cbt)的叶节点，并以事件发生的时间为权值进行排序。基于堆排序，cbt的叶节点以递
归的形式向上逐步更新，最后得到最早发生的事件。hs子区每次只处理这个最早发生的有效事件，并将演化时间相应地更新到当前时刻。待处理完该事件，可能产生一些新事件，也有可能使得未发生的事件失效。因此，需要重新预测粒子所在的网格及周围网格单元内事件发生的时间，重新构建网格的未来事件列表，并且更新相应的cbt叶节点，进而通过更新cbt得到演化系统需要处理的下一个最早发生的有效事件。
[0119]
在步骤s500中，开始模拟计算过程，检测粒子a、粒子b、粒子c和粒子c’各自的反应状态、粒子半径和粒子位置的变化过程，并统计各个粒子的数量信息，当粒子a或粒子b的数量为0时，输出计算结果。
[0120]
在该步骤中，开启模拟计算过程，随着时间的推移，粒子a和粒子b会发生反应，进而生成粒子c’，随后粒子c’逐渐变化至粒子c。在该过程中，各个粒子的反应状态、粒子位置等会不断变化，粒子a和粒子b的粒子半径不发生改变，粒子c’的粒子半径不断变化。各个粒子数量也会在该过程中不断改变，直到粒子a或粒子b的数量为0，反应结束，模拟计算过程终止。
[0121]
在一些实施例中，步骤s500包括：
[0122]
步骤s510，检测是否有粒子发生碰撞，当检测到有两个粒子发生碰撞时，判断发生碰撞的粒子的种类；
[0123]
当发生碰撞的粒子为粒子a和粒子b时，根据预设的反应原则，判断发生碰撞的粒子a和粒子b是否发生反应，并根据判断结果做出对应的处理，反应原则根据反应参数进行设定；当发生碰撞的粒子不是粒子a和粒子b时，则更新粒子碰撞后各自对应的速度、粒子半径和粒子位置；
[0124]
步骤s520，统计各类粒子的数量信息，当粒子a或粒子b的数量为0时，输出计算结果。
[0125]
在步骤s510中，检测是否有粒子发生碰撞，当检测到有两个粒子发生碰撞时，判断发生碰撞的粒子的种类。
[0126]
在a+b
→
c整个化学反应体系中，可能存在粒子a、粒子b、粒子c’和粒子c多种类型的粒子，各个粒子间均有可能发生碰撞。本公开中，针对不同的粒子发生碰撞后，相应地有不同的处理方式。具体地，当发生碰撞的粒子为粒子a和粒子b时，根据预设的反应原则，判断发生碰撞的粒子a和粒子b是否发生反应，并根据判断结果做出对应的处理，反应原则根据反应参数进行设定；当发生碰撞的粒子不是粒子a和粒子b时，则更新粒子碰撞后各自对应的速度、粒子半径和粒子位置。
[0127]
判断发生碰撞的粒子a和粒子b是否发生反应，可根据预设的反应原则进行判断。预设的反应原则可根据反应参数进行设定。在一具体实施例中，反应参数包括反应临界能，基于反应临界能设定反应原则。例如，只有反应物粒子a和粒子b在质心连线方向的相对平动能之和大于设定的反应临界能，才认为粒子a和粒子b要发生化合为粒子c
′
的反应。
[0128]
本公开针对粒子a和粒子b是否发生反应，对应有不同的处理方式。
[0129]
在一具体实施例中，判断发生碰撞的粒子a和粒子b是否发生反应，并根据判断结果做出对应的处理包括：
[0130]
步骤s511，若判断粒子a和粒子b发生反应，则删除发生碰撞的粒子a和粒子b中的其中一者，并将另一者转存为粒子c’，且更新粒子a、粒子b和/或粒子c’在碰撞后各自对应
的速度、粒子半径和粒子位置，以及更新硬球模拟子区的未来事件列表；
[0131]
步骤s512，若判断粒子a和粒子b未发生反应，则更新粒子a和粒子b碰撞后各自对应的速度、粒子半径和粒子位置。
[0132]
在步骤s511中，若判断粒子a和粒子b发生反应，则删除发生碰撞的粒子a和粒子b中的其中一者，并将另一者转存为粒子c’，且更新粒子a、粒子b和/或粒子c’在碰撞后各自对应的速度、粒子半径和粒子位置，以及更新硬球模拟子区的未来事件列表。
[0133]
针对粒子a和粒子b发生反应生成粒子c’，粒子a、粒子b和粒子c’满足第一关系式
[0134][0135]
以及第二关系式
[0136][0137]
其中，e为结余能量，e＞0，m
a
和m
b
分别为粒子a和粒子b在碰撞时的粒子质量，v
a
、v
b
分别为粒子a和粒子b在碰撞时的速度，m
c’为发生反应时粒子c’的粒子质量，v
c
为发生反应时粒子c’的速度。
[0138]
根据上述式(9)可以确定粒子c’的初始速度，式(10)可用于确定粒子c’在后续碰撞过程中的发生碰撞后的速度。
[0139]
具体在一实施例中，将e做为粒子的一个物理属性引入粒子c
′
的存储信息的内存块中。即，在每一个粒子的物理属性中添加能量项energy。一般情况下，粒子的energy＝0，但是对于发生化合反应新生成的粒子c
′
，其energy＝e。储存在粒子c
′
中的结余能量e，可以在粒子c
′
与其他粒子碰撞时，一次释放掉e或者按一定的比例在多次碰撞过程中释放掉e。处理粒子碰撞事件时，需要首先检测粒子的energy。若energy>0，需要在能量守恒方程左边加入粒子现有的energy值，同时在方程的右边加入粒子剩余的energy值。然后，通过求解能量守恒和动量守恒方程，可以得到粒子碰撞后的速度。以一次碰撞全部释放e为例，动量守恒方程为：
[0140]
m1v1+m2v2＝m1v'1+m2v'2ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(11)
[0141]
反应前后的能量守恒方程可描述为：
[0142][0143]
求解方程(11)和(12)，得到碰撞后粒子的速度为：
[0144][0145][0146]
其中，j表示粒子碰撞时的冲量的单位矢量，可表示为：
[0147][0148]
j表示粒子碰撞时的冲量，通过求解方程(11)和(12)，得到表达式为：
[0149][0150]
根据求根公式，得到δ的表达式为：
[0151][0152]
由此可见，守恒方程中包含了结余能量e项，得不到显示的v1'和v2'表达式。
[0153]
进一步地，在步骤s511中，粒子c
′
的粒子半径变化满足第三关系式，
[0154][0155]
其中，t0和t1分别表示粒子运动的上一个时刻及当前时刻，和分别对应于粒子在t0和t1时刻的粒子半径，v
r
为粒子c’的粒子半径增长速率，v
r
>0和v
r
<0分别表示粒子的粒子半径变化过程是一个逐渐增大和逐渐减小的过程。
[0156]
在该步骤中，为避免新生成的粒子与周围的粒子发生位置重叠，在粒子属性中引入了粒子半径增长速率参数(v
r
)。首先比较反应物粒子a和b的粒子半径σ
r
，并删除σ
r
较大的反应物粒子。例如，若满足σ
ra
<σ
rb
，则将粒子a的位置坐标和粒子半径设置为生成物粒子c
′
的位置坐标和粒子半径，以避免粒子的σ
r
突然变化的情况发生。这样处理的优点是不需要检测生成物是否与周围现有的粒子发生重叠，大大节省了因检索合适的生成物位置而产生的计算量。
[0157]
在步骤s511中，当粒子a和粒子b发生反应，可比较粒子a和粒子b反应前的粒子半径的大小，删除其中粒子半径较大的一者，并将其中粒子半径较小的一者转存为粒子c’。同时，该粒子的粒子属性参数也转存为粒子c’的粒子属性参数。具体将存储粒子b(或a)的物理属性的内存块转换为存储粒子c
′
的物理属性并进行相应的修改。
[0158]
针对pp子区，基于反应物的位置坐标，该区域对删除粒子a或粒子b，转存为粒子c’，分两种情况进行处理。一种情况是反应物粒子a和粒子b处在同一个进程中，可以采用类似于hs模型的处理方式，在当前时间步完成粒子的删除操作。这是因为，化合反应中粒子数量的变化是删除粒子，而不是添加粒子，不会出现新添加的粒子不在当前进程中的情况。另外一种情况是，反应物粒子a和粒子b不在同一个进程中，针对这种情况的处理方式是，在当前时间步内比较粒子a和粒子b反应前的粒子半径的大小，并对其中粒子半径较大的粒子予以标记，在下一个时间步同步演化拟颗粒子区所有粒子的状态之前，检测所有粒子的标记信息，删除被予以标记的粒子。具体地，在存储每一个粒子的物理属性中引入一个标记comflag，并且一般情况下粒子的comflag＝0，只有需要被删除的粒子的comflag＝1。当前时间步处理化合反应时，首先比较反应物的粒子半径σ
r
，并将σ
r
较大的反应物粒子标记为将要被删除的粒子。例如，若满足σ
ra
<σ
rb
，则设置粒子b的comflag＝1。同时，将内存块中存储的反应物粒子a的物理属性修改为存储粒子c
′
的物理属性。然后，根据公式(9)得到粒子c
′
的速度，并将多余的能量e储存在粒子c
′
中。这样就完成了在当前时间步内对化合反应的所有相关操作。下一个时间步进入pp子区后，在同步演化pp子区所有粒子的状态之前，首先检测粒子的comflag。若粒子的comflag＝1，则将该粒子从其所在的网格单元的粒子链表中彻底地删除。这种处理方式可避免在相邻进程间传递粒子信息时，出现发送和接收的粒子数量不一致的情况。此外，该处理方式简化了计算过程，提高了计算效率。
[0159]
在此需说明的是，若不考虑计算成本和计算时间，pp子区也可以完全采用与hs子区一样的处理方式，在当前时间步完成粒子的删除操作。此外，本公开为了节省计算开销，当粒子a和粒子b发生化合反应时，只需要删除粒子a(或b)，将存储粒子b(或a)的物理属性的内存块转换为存储粒子c
′
的物理属性并进行相应的修改，而不必从粒子链表中彻底删除反应物粒子a和b。实际上，若不计额外的计算开销，也可以从粒子链表中彻底地删除存储粒子a和b的物理属性的内存块，同时增加存储粒子c
′
的物理属性的内存块。
[0160]
在步骤s511中，更新硬球模拟子区的未来事件列表包括：
[0161]
删除与粒子a和粒子b有关的无效事件，增加与粒子c’有关的事件。例如，增加与粒子c’有关的碰撞事件或粒子半径变化事件。
[0162]
在一具体实施例中，步骤s510中，判断发生碰撞的粒子a和粒子b是否发生反应，并根据判断结果做出对应的处理还包括：
[0163]
步骤s513，监测粒子c’的运动状态和粒子半径，当粒子c’的粒子半径等于粒子c的硬球半径时，将粒子c’的粒子类型编号更改为粒子c对应的粒子类型编号。
[0164]
随着模拟过程的进行，粒子c’可能与其他粒子发生碰撞，在碰撞过程中，粒子c’的半径会不断发生改变。
[0165]
在步骤s513中，hs子区和pp子对粒子c’的粒子半径的处理方式不同。
[0166]
对于hs子区，在化合反应体系后续的模拟演化过程中，hs模型将根据粒子现有的粒子半径σ
r
与粒子c的粒子半径的大小关系，相应地增大或者减小该粒子的σ
r
。即，该粒子在之后与其他粒子发生碰撞时将逐渐变化其σ
r
，直到满足关系式并将生成物前驱体的粒子类型更新为相应生成物的粒子类型。因此，在hs模型的事件类型列表中引入了粒子的σ
r
变化到正好等于的事件类型。处理该类型事件时，不仅需要更新粒子类型及σ
r
，还需要设置v
r
＝0。因为hs模型受事件驱动，体系模拟演化过程中不可能发生或者的情况。这是因为针对某个粒子，发生事件的时间一定早于发生和事件的时间，hs模型将优先处理发生时间更早的的事件。
[0167]
对于pp子区，由于ppm允许粒子之间发生重叠，在同步更新所有粒子的状态之前，首先根据v
r
更新粒子的σ
r
。通过检测粒子的v
r
是否为0，判断是否需要改变粒子的σ
r
。只有当粒子的v
r
≠0时，才需要相应地更新粒子的σ
r
到时刻t1。由于ppm受时间驱动，有可能出现或者的情况。无论出现任何一种情况都需要在当前时间步设置并修改粒子类型参数和设置粒子的v
r
＝0。
[0168]
在此需说明的是，本公开为了简化计算过程，进行化合反应时从粒子链表中删除一个反应物粒子，并将另外一个反应物粒子标记为生成物粒子c
′
，粒子c
′
的粒子半径在随后与体系内其他粒子碰撞的过程中基于粒子半径变化速率而逐渐发生变化直到粒子半径c
′
＝c。实际上，若不考虑计算开销，也可以直接生成粒子c而不是粒子c
′
。
[0169]
在步骤s510中，当发生碰撞的粒子不是粒子a和粒子b时，则更新粒子碰撞后各自对应的速度、粒子半径和粒子位置。
[0170]
具体地，当粒子间发生碰撞时，基于碰撞前后动量和能量的守恒原则，更新粒子碰撞时的速度：
[0171][0172][0173]
其中，e、m、r、v分别表示碰撞的恢复系数、粒子质量、粒子位置及速度，一般在刚性碰撞中e＝1.0，发生碰撞的粒子由下标1和2进行区分。
[0174]
本公开针对简单典型反应类型中普遍存在的化合反应类型a+b
→
c，开发了一种适用于应用在hs
‑
ppm耦合算法框架中的计算方法，适用于对包含化合反应类型的纳微尺度上的化学反应体系进行高效和大规模的并行模拟计算，具有可并行计算、计算效率高、计算规模大、适用范围广、能够获得全系统演化的信息等优点。
[0175]
下面将结合具体实施例说明本公开提供的模拟方法。
[0176]
以下实施例均针对沸石分子筛sapo
‑
34体系进行模拟研究。sapo
‑
34是具有类菱沸石型结构的分子筛，其骨架主要由硅铝磷酸盐组成，具有双六元环、八元环与四元环构成的椭球形cha笼以及三维交叉孔道结构。sapo
‑
34分子筛既具有特殊的吸水性能和质子酸性，同时又有良好的水热稳定性，使得它可用作吸附剂和催化剂载体。可用于小分子的择形催化反应，比如作为甲醇制烯烃的催化剂，还可以作为汽车尾气净化催化剂载体。当发生碰撞的两个反应物粒子在质心连线方向的相对平动能之和大于设定的反应的临界能量时，反应物即按照：
[0177]
a+b
→
c
ꢀꢀꢀꢀ
(反应r1)
[0178]
a+c
→
d
ꢀꢀꢀꢀ
(反应r2)的化合反应类型发生反应。反应r1和r2的临界能量分别标记为e
a1
和e
a2
。初始反应体系中存在大量的反应物粒子a和b，可认为反应r1为主反应，而反应r2可被认为是副反应。周期性边界条件应用在反应体系的x、y和z三个方向。当某粒子从某个方向离开反应体系时，将以相同的粒子类型从对称的位置进入反应体系。为便于讨论分析，在实施例中使用粒子a的质量m
a
(2.66
×
10
‑
26
kg)、(0.341nm)及特征时间(0.322ps)对质量、长度及时间的单位分别进行无量纲化处理，并且将活化能e
a1
无量纲化处理为1.0。
[0179]
实施例1：
[0180]
反应r1和r2的临界能量分别为e
a1
＝1.00和e
a2
＝1.20。体系中初始的反应物为粒子a和b，化合反应后的生成物为粒子c和d。模拟体系演化之前，反应物粒子a和b随机地分布在沸石分子筛的孔道中。设置模拟体系的气体负载量为平均每个沸石分子筛晶胞单元中含有4个气体分子。构建的反应模型中，沸石晶胞数量为8
×8×
8晶胞。模拟体系的反应温度设置为623k。不同粒子的质量m和有效硬球半径σ
hs
如表1所示。基于反应前后质量守恒设定粒子的质量m。粒子的σ
hs
根据粒子的m的变化而等比例发生变化。即，粒子的m越大，其σ
hs
也越大。模拟体系初始化时，各粒子的数量如表1中的n0所示。待完成模拟体系初始化后，发起程序开始计算，当检测到体系中粒子a或者粒子b的分子数量为0时，即认为完成了所有反应，体系演化结束，反应后的粒子数量如表1中的n1所示。模拟结果表明，反应前后粒子数量的变化符合预设的反应关系，验证了算法实现的化合反应的正确性。
[0181]
表1模拟体系中粒子的质量、有效硬球半径、反应前数量及反应后数量
[0182][0183]
注：m、σ
hs
、n0、n1分别是粒子的质量、有效硬球半径、反应前粒子数量及反应后粒子数量。
[0184]
实施例2：
[0185]
反应r1和r2的临界能量分别为e
a1
＝1.00和e
a2
＝1.50。体系中初始的反应物为粒子a和b及少量的粒子c，化合反应后的生成物为粒子c和d。模拟体系演化之前，反应物粒子a、b及c随机地分布在沸石分子筛的孔道中。设置模拟体系的气体负载量为平均每个沸石分子筛晶胞单元中含有5个气体分子。构建的反应模型中，沸石晶胞数量为9
×9×
9晶胞。模拟体系的反应温度设置为723k。不同粒子的质量m和有效硬球半径σ
hs
如表2所示。基于反应前后质量守恒设定粒子的质量m。粒子的σ
hs
根据粒子的m的变化而等比例发生变化。即，粒子的m越大，其σ
hs
也越大。模拟体系初始化时，各粒子的数量如表2中的n0所示。待完成模拟体系初始化后，发起程序开始计算，当检测到体系中粒子a或者粒子b的分子数量为0时，即认为完成了所有反应，体系演化结束，反应后的粒子数量如表2中的n1所示。模拟结果表明，反应前后粒子数量的变化符合预设的反应关系，验证了算法实现的化合反应的正确性。
[0186]
表2模拟体系中粒子的质量、有效硬球半径、反应前数量及反应后数量
[0187][0188]
注：m、σ
hs
、n0、n1分别是粒子的质量、有效硬球半径、反应前粒子数量及反应后粒子数量。
[0189]
实施例3：
[0190]
反应r1和r2的临界能量分别为e
a1
＝1.00和e
a2
＝1.40。体系中初始的反应物为粒子a和b，化合反应后的生成物为粒子c和d。模拟体系演化之前，反应物粒子a和b随机地分布在沸石分子筛的孔道中。设置模拟体系的气体负载量为平均每个沸石分子筛晶胞单元中含有6个气体分子。构建的反应模型中，沸石晶胞数量为10
×
10
×
10晶胞。模拟体系的反应温度设置为673k。不同粒子的质量m和有效硬球半径σ
hs
如表3所示。基于反应前后质量守恒设定粒子的质量m。粒子的σ
hs
根据粒子的m的变化而等比例发生变化。即，粒子的m越大，其σ
hs
也越大。模拟体系初始化时，各粒子的数量如表3中的n0所示。待完成模拟体系初始化后，发起程序开始计算，当检测到体系中粒子a或者粒子b的分子数量为0时，即认为完成了所有反应，体系演化结束，反应后的粒子数量如表3中的n1所示。模拟结果表明，反应前后粒子数量的变化符合预设的反应关系，验证了算法实现的化合反应的正确性。
[0191]
表3模拟体系中粒子的质量、有效硬球半径、反应前数量及反应后数量
[0192][0193]
注：m、σ
hs
、n0、n1分别是粒子的质量、有效硬球半径、反应前粒子数量及反应后粒子数量。
[0194]
实施例4：
[0195]
反应r1和r2的临界能量分别为e
a1
＝1.00和e
a2
＝1.60。体系中初始的反应物为粒子a、b及少量的d，化合反应后的生成物为粒子c和d。模拟体系演化之前，反应物粒子a、b及c随机地分布在沸石分子筛的孔道中。设置模拟体系的气体负载量为平均每个沸石分子筛晶胞单元中含有6个气体分子。构建的反应模型中，沸石晶胞数量为12
×
12
×
12晶胞。模拟体系的反应温度设置为773k。不同粒子的质量m和有效硬球半径σ
hs
如表4所示。基于反应前后质量守恒设定粒子的质量m。粒子的σ
hs
根据粒子的m的变化而等比例发生变化。即，粒子的m越大，其σ
hs
也越大。模拟体系初始化时，各粒子的数量如表4中的n0所示。待完成模拟体系初始化后，发起程序开始计算，当检测到体系中粒子a或者粒子b的分子数量为0时，即认为完成了所有反应，体系演化结束，反应后的粒子数量如表4中的n1所示。模拟结果表明，反应前后粒子数量的变化符合预设的反应关系，验证了算法实现的化合反应的正确性。
[0196]
表4模拟体系中粒子的质量、有效硬球半径、反应前数量及反应后数量
[0197][0198]
注：m、σ
hs
、n0、n1分别是粒子的质量、有效硬球半径、反应前粒子数量及反应后粒子数量。
[0199]
需要说明的是，尽管在附图中以特定顺序描述了本公开中方法的各个步骤，但是，这并非要求或者暗示必须按照该特定顺序来执行这些步骤，或是必须执行全部所示的步骤才能实现期望的结果。附加的或备选的，可以省略某些步骤，将多个步骤合并为一个步骤执行，以及/或者将一个步骤分解为多个步骤执行等，均应视为本公开的一部分。
[0200]
应可理解的是，本公开不将其应用限制到本说明书提出的部件的详细结构和布置方式。本公开能够具有其他实施方式，并且能够以多种方式实现并且执行。前述变形形式和修改形式落在本公开的范围内。应可理解的是，本说明书公开和限定的本公开延伸到文中和/或附图中提到或明显的两个或两个以上单独特征的所有可替代组合。所有这些不同的组合构成本公开的多个可替代方面。本说明书的实施方式说明了已知用于实现本公开的最佳方式，并且将使本领域技术人员能够利用本公开。