信息处理设备以及用于粒子模拟的方法和系统

技术领域

本文讨论的实施方式涉及信息处理设备、粒子模拟方法和粒子模拟系统。

背景技术

离散元素方法(Discrete element method，DEM)，也称为不同元素方法，是一种用于分析对象的行为的模拟方案。DEM方案将要分析的对象建模为呈例如圆盘或球形形状的粒子的集合，并且考虑到粒子之间的接触和滑动来跟踪粒子的运动。DEM方案有时用于分析粉末例如沉积物颗粒、化学物质和打印机墨粉粒子的行为。

DEM方案检测给定粒子与之接触的粒子(particle)，然后计算与该粒子接触的所有相邻粒子所施加的接触力，并基于接触力找到后来的预定时间的粒子的位置。在模拟过程中采用具有预定时间段的时间步(time step)，并且计算循环在接触的相邻粒子的检测、接触力的计算以及粒子的位置的更新之间交替。因此，DEM方案模拟了每个粒子的位置沿时间轴的变化。

提出了一种用于实现节省与DEM模拟相关联的计算工作量的分析设备。为了降低计算复杂度，分析设备将要分析的区域划分为多个单元，然后如果要检查的粒子与目标粒子接触(该过程被称为“接触检测”)，则将该要检查的粒子限制为与包含目标粒子的单元邻接的相邻单元中的粒子。模拟系统是另一种提出的可以节省DEM模拟涉及的计算工作量的技术。为此，模拟系统生成指示在特定时间段期间可能与目标粒子接触的粒子的相邻列表，并在每个时间步仅针对相邻列表中包括的粒子执行接触检测。还提出另一种模拟系统，该系统将粒子分配给多个节点并且并行运行DEM计算。

参见，例如，日本公开特许公报第2010-72379号；国际公布小册子第WO 2012/034176号；以及国际公布小册子第WO 2014/094410号。

在DEM方案中，为了计算由另一接触的粒子施加在给定粒子上的接触力，有时使用两个粒子接触期间的累积位移。例如，在计算分解为法向分量和切向分量的接触力时，切向接触力可以取决于在形成接触之后直到当前时间的沿切线方向的累积相对行进距离。因此，当模拟时钟上的时间前进时，将发生复制处理，其中针对每对连续接触的粒子，将先前时间的位移复制到当前时间。

涉及识别给定粒子在先前时间点和当前时间点都与之接触的粒子并且复制所识别的接触的粒子中的每一个的位移的复制处理是一项计算量很大的任务。例如，最简单的复制处理方法之一是基于各个粒子在当前时间点的位置提取在当前时间点与给定粒子接触的粒子，然后从指示在先前时间点处于接触的粒子的接触数据中搜索这些提取的粒子。该方法的缺点在于重复扫描先前时间点的接触数据，并且这样的比较任务变得在计算上是密集的。

鉴于此，改变复制处理的过程从而降低计算复杂度可能被认为是合理的。然而，改变并行粒子模拟中的复制处理过程并不容易，在所述并行粒子模拟中，粒子可能分布的整个区域被划分为子区域，并且多个节点被分配不同的子区域以进行并行执行。

在并行粒子模拟中，属于子区域的粒子可能会在两个子区域之间的边界处与属于相邻子区域的其他粒子接触。因此，多个节点需要交换关于在其边界附近的粒子的数据。然而，应该注意，每个节点都能够有效地跟踪粒子在其分配的子区域中而不是在相邻子区域中的运动。出于该原因，如果以统一的方式执行复制处理而不区分每个接触粒子是属于其分配的子区域还是相邻的子区域，则很难降低计算复杂度。

发明内容

实施方式的一方面在于提供实现粒子模拟的计算复杂度的降低的信息处理设备、粒子模拟系统和粒子模拟方法。

根据一个方面，提供了一种信息处理设备，包括：存储单元，其用于存储与属于第一区域的第一粒子相关联的第一接触数据，第一接触数据指示在第一时间处与第一粒子接触的一个或更多个接触粒子以及在形成接触之后直到第一时间的接触时段上发生的第一粒子相对于接触粒子中的每个接触粒子的累积位移；以及处理单元，其用于执行包括以下的处理：计算指示在第一时间之后的第二时间处属于第一区域的粒子的位置的第一位置数据，以及从不同的信息处理设备获取指示在第二时间处属于第二区域的粒子的位置的第二位置数据，第二区域与第一区域接界；基于从第一接触数据和第一位置数据获得的计算结果来检测第二粒子，第二粒子在第一时间处与第一粒子接触并属于第一区域并且在第二时间处属于第一区域；基于第一位置数据检查第一粒子和第二粒子在第二时间处是否彼此接触，并且在第一粒子和第二粒子在第二时间处接触的情况下，将第二粒子的累积位移从第一接触数据复制到与第二时间相对应的第二接触数据；基于第一位置数据和第二位置数据检测第三粒子，第三粒子是除了第二粒子以外的各自属于第一区域或第二区域并且在第二时间处与第一粒子接触的粒子；以及在第一接触数据中搜索第三粒子，并且当每个第三粒子被登记在第一接触数据中时，将所登记的第三粒子的累积位移从第一接触数据复制到第二接触数据。

附图说明

图1示出了根据第一实施方式的示例性信息处理系统；

图2示出了根据第二实施方式的示例性信息处理系统；

图3是示出计算节点的示例性硬件配置的框图；

图4示出了DEM方案中的示例性粒子模型；

图5示出了模拟区域的示例性划分；

图6是示出相邻粒子的示例性运动的第一图；

图7是示出相邻粒子的示例性运动的第二图；

图8示出了沿切线方向的位移的示例性复制；

图9是示出计算节点的示例性功能的框图；

图10示出了示例性结果表；

图11示出了示例性粒子表和粒子编号映射表；

图12示出了示例性接触粒子列表；

图13是示出DEM方案的示例性过程的流程图；

图14是示出接触粒子列表的生成的示例性过程的流程图的第一部分；

图15从图14继续，是示出接触粒子列表的生成的示例性过程的流程图的第二部分；以及

图16从图15继续，是示出接触粒子列表的生成的示例性过程的流程图的第三部分。

具体实施方式

下面将参照附图来描述若干实施方式。

(a)第一实施方式

在下文中描述第一实施方式。

图1示出了根据第一实施方式的示例性信息处理系统。

根据第一实施方式的信息处理系统使用DEM方案执行用于模拟粉末粒子的运动的粒子模拟。第一实施方式的信息处理系统包括多个信息处理器，信息处理器包括信息处理器10和10-1。第一实施方式的信息处理系统也可以被称为并行处理系统、模拟系统或模拟器。在粒子模拟的并行处理中，粒子可能分布的模拟空间被划分为多个区域，并且不同的信息处理器被分配不同的区域。包括信息处理器10和10-1的信息处理器中的每一个模拟属于其分配的区域(目标区域)的粒子的运动。每个信息处理器也可以称为计算机、节点或模拟单元。多个信息处理器连接至网络，从而能够交换数据。如下所描述的，信息处理器10负责区域13，信息处理器10-1负责区域13-1。

信息处理器10包括存储单元11和处理单元12。存储单元11可以是易失性半导体存储器诸如随机存取存储器(RAM)，或者可以是非易失性存储装置诸如硬盘驱动器(HDD)或闪存。处理单元12例如是诸如中央处理单元(CPU)、图形处理单元(GPU)或数字信号处理器(DSP)的处理器。然而，应该注意，处理单元12可以包括针对特定用途而设计的电子电路，例如专用集成电路(ASIC)或现场可编程门阵列(FPGA)。处理器执行存储在诸如RAM的存储器中(或存储单元11中)的程序。术语“多处理器”或简单地“处理器”可以用于指代多个处理器的集合。

根据DEM方案，信息处理器10重复以下过程：检查在当前时间点与目标粒子接触的接触粒子；计算接触粒子施加在目标粒子上的接触力；以及基于接触力计算目标粒子在下一时间点的位置。在每个模拟时间点，信息处理器10生成指示接触粒子的接触数据。

接触力的计算使用在目标粒子和每个接触粒子彼此接触期间发生的累积位移。累积位移是例如从形成接触的时间开始在每个时间步上目标粒子相对于每个接触粒子沿切线方向的增量行进距离的总和。因此，信息处理器10在接触数据中记录每个模拟时间点的累积位移。此时，对于与目标粒子保持接触的每个接触粒子，需要将记录在先前时间点的接触数据中的累积位移复制到当前时间点的接触数据中，然后再将当前时间点的增量位移添加到复制的累积位移。

存储单元11在其中存储接触数据14(第一接触数据)。存储单元11还可以存储将在下面描述的位置数据15(第一位置数据)、位置数据15-1(第二位置数据)和接触数据16(第二接触数据)。

接触数据14是与属于被分配给信息处理器10的区域13的粒子J1(第一粒子)相关联的模拟时间点t1(第一时间)的接触数据。接触数据16是与粒子J1相关联的在时间点t1之后的模拟时间点t2(第二时间)的接触数据。时间点t2例如比时间点t1晚预定单位时间出现。接触数据14指示在时间点t1与粒子J1接触的接触粒子，以及在形成接触之后直到时间点t1粒子J1相对于接触粒子的位移。接触数据16指示在时间点t2与粒子J1接触的接触粒子，以及在形成接触之后直到时间点t2粒子J1相对于接触粒子的位移。例如，接触数据14和16中的每一个包括接触粒子的标识信息和与接触粒子相关联的位移。如果存在多个与粒子J1接触的接触粒子，则接触数据14和16列举例如彼此映射的、接触粒子的多个标识信息和多个位移。

位置数据15指示在时间点t2处属于信息处理器10负责的区域13的粒子的位置。位置数据15-1指示在时间点t2处属于信息处理器10-1负责的区域13-1的粒子的位置。区域13-1与区域13接界。位置数据15-1可以仅包括属于区域13-1的粒子之中的在距区域13与13-1之间的边界预定距离内的一些粒子。位置数据15由信息处理器10生成，而位置数据15-1由信息处理器10-1生成。

处理单元12将模拟时钟上的时间从时间点t1前进到时间点t2。那时，处理单元12生成位置数据15。例如，处理单元12基于时间点t1的接触数据例如接触数据14计算作用在属于区域13的每个粒子上的接触力，然后基于计算出的接触力和指示区域13中的每个粒子在时间点t1的位置的位置数据来计算粒子在时间点t2的位置。另外，处理单元12获取从信息处理器10-1接收的位置数据15-1。位置数据15-1由信息处理器10-1以与位置数据15相同的方式生成。

当模拟时钟上的时间从时间点t1前进到时间点t2时，某些粒子可能会越过区域13与13-1之间的边界到达另一区域。在时间点t1属于区域13的粒子可能在时间点t2仍保留在区域13中或已移动到区域13-1。类似地，在时间点t1属于区域13-1的粒子可能在时间点t2仍保留在区域13-1中或已移动到区域13。每个粒子由负责粒子的目的地区域的信息处理器管理。

接下来，处理单元12生成时间点t2的接触数据16。那时，处理单元12基于根据接触数据14和位置数据15获得的计算结果来检测在时间点t1处与粒子J1接触并属于区域13然后在时间点t2处仍然属于区域13的第二粒子。

例如，处理单元12针对时间点t1和t2向每个粒子分配不同标识信息。在这方面，当计算位置数据15时，处理单元12生成映射数据，该映射数据针对属于区域13的每个粒子将时间点t1的标识信息映射到时间点t2的标识信息。因为区域13是信息处理器10跟踪粒子的运动的唯一区域，所以映射数据包括在时间点t1和t2均属于区域13的粒子。即，映射数据不包括已经从区域13移出到区域13-1中的粒子、已经从区域13-1移出到区域13中的粒子以及保留在区域13-1中的粒子。在这种情况下，处理单元12在映射数据中搜索映射到在接触数据14中包括的时间点t1的每个标识信息的时间点t2的标识信息，从而检测满足上述条件的第二粒子。

处理单元12基于位置数据15检查粒子J1和前述第二粒子在时间点t2是否彼此接触。如果它们接触，则处理单元12将与第二粒子相关联的位移从接触数据14复制到接触数据16。例如，处理单元12检测在时间点t1与粒子J1接触并且属于区域13然后在时间点t2属于区域13的粒子J2。在时间点t2，粒子J1和J2彼此接触。因此，处理单元12将与粒子J2相关联的位移D2从接触数据14复制到接触数据16。

接下来，处理单元12基于位置数据15和15-1检测第二粒子以外的各自属于区域13或13-1并且在时间点t2与粒子J1接触的第三粒子。处理单元12在接触数据14中搜索第三粒子。如果每个以上描述的第三粒子被登记在接触数据14中，则处理单元12将与该第三粒子相关联的位移从接触数据14复制到接触数据16。

例如，处理单元12检测在时间点t2与粒子J1接触的粒子J3-1、J3-2和J3-3。粒子J3-1在时间点t1属于区域13，并且在时间点t2已经移动到区域13-1。粒子J3-2在时间点t1属于区域13-1，并且在时间点t2已经移动到区域13。粒子J3-3在时间点t1和t2都属于区域13-1。处理单元12在接触数据14中搜索各个粒子J3-1、J3-2和J3-3。粒子J3-1、J3-2和J3-3被登记在接触数据14中。因此，处理单元12将分别与粒子J3-1、J3-2和J3-3相关联的位移D3-1、D3-2和D3-3从接触数据14复制到接触数据16。注意，接触数据14和16可以区分属于区域13的接触粒子和属于区域13-1的接触粒子。

根据第一实施方式的信息处理系统，在粒子模拟中，生成接触数据16，该接触数据16指示在时间点t2与粒子J1接触的接触粒子以及分别与这些接触粒子相关联的位移。这使得能够有效地计算在时间点t2由接触粒子施加在粒子J1上的接触力，这进而便于有效地计算粒子J1在时间点t2之后的时间点t3的位置。

另外，当生成接触数据16时，利用时间点t1的接触数据14来检测在时间步t1处与粒子J1接触并属于区域13然后在时间点t2保持在区域13中的第二粒子。对于第二粒子，在确认第二粒子在时间点t2与粒子J1接触之后，将与第二粒子相关联的位移从时间点t1的接触数据14复制到时间点t2的接触数据16。这消除了在接触数据14中重复搜索区域13中保留的粒子的需要，从而降低了由比较任务导致的计算复杂度。因此，可以提供将位移从接触数据14复制到接触数据16的更有效的方式。

此外，除了上述第二粒子之外，还检测在时间点t2与粒子J1接触的第三粒子。第三粒子包括已经从区域13移动到区域13-1中的粒子；已经从区域13-1移动到区域13中的粒子；保留在区域13-1中的粒子。对于这些第三粒子，从接触数据14中进行搜索。然后，对于在接触数据14中找到的每个第三粒子，将与第三粒子相关联的位移从接触数据14复制到接触数据16。

如上所描述的，针对保留在区域13中的接触粒子和其他接触粒子采用不同的复制过程。这是因为，对于属于区域13的粒子，由于区域13是信息处理器10的目标区域，因此它们的运动被有效地跟踪；相反，对于属于相邻区域13-1的粒子，信息处理器10难以有效地跟踪它们的运动。根据第一实施方式，将保留在区域13中的接触粒子与其他接触粒子分离会使得复制过程有效，至少对于前者接触粒子会使得复制过程有效。这致使降低了复制处理中的计算复杂度。

(b)第二实施方式

接下来描述第二实施方式。

图2示出了根据第二实施方式的示例性信息处理系统。

根据第二实施方式的信息处理系统是用于使用DEM方案对大量粒子的运动进行模拟的并行处理系统。考虑到粒子之间的接触，DEM方案跟踪每个粒子的位置和速度的变化。DEM方案以每个粒子都放置在其初始位置的第一时间步开始，并且模拟时钟以非常短的持续时间的离散时间步前进。沿时间轴在每个时间步更新每个粒子的位置和速度。在DEM方案中建模的粒子的几何形状可以是二维(2D)圆盘或三维(3D)球体。然而，注意，为了便于说明，第二实施方式假定模拟的粒子是二维圆盘。

根据第二实施方式的信息处理系统包括多个计算节点例如计算节点100、100-1、100-2和100-3以及管理节点200。多个计算节点和管理节点200各自连接至网络30，该网络30是数据通信网络例如局域网(LAN)。注意，计算节点100对应于第一实施方式的信息处理器10。计算节点100-1对应于第一实施方式的信息处理器10-1。

多个计算节点是用于以分散的方式计算每个时间步的大量粒子的位置和速度的服务器计算机。粒子可能分布的整个模拟区域被划分成多个子区域，并且每个子区域被分配给一个处理。DEM方案中使用的每个计算节点执行一个子区域被分配至的一个处理。因此，使用多个计算节点来并行执行与多个子区域相对应的多个处理。为了便于说明，第二实施方式假设整个区域被划分为四个子区域并且使用四个计算节点100、100-1、100-2和100-3并行执行四个处理的情况。然而，注意，可以增加子区域的数量并使用增加数量的计算节点。

计算节点100、100-1、100-2和100-3中的每一个前进一个时间步，并且更新存在于其被分配的子区域中的粒子的位置和速度。作为结果，在先前时间步期间存在于被分配的子区域中的粒子可能会移出被分配的子区域。鉴于此，计算节点100、100-1、100-2和100-3在每个时间步彼此通信并且交换已经跨过边界移动到相邻子区域的粒子的粒子数据。另外，在分配的子区域中的粒子之中，位于与相邻子区域的边界附近的那些粒子可能与相邻子区域中的粒子接触。因此，计算节点100、100-1、100-2和100-3在每个时间步彼此通信，并且复制在其边界附近的粒子的粒子数据并共享复制的粒子数据。注意，第二实施方式着重于计算节点100来描述DEM处理。

管理节点200控制采用多个计算节点的DEM方案的并行处理。具体地，管理节点200从用户接收包括各个粒子的初始位置的输入数据。管理节点200将整个模拟区域划分为多个子区域，然后根据各个粒子的初始位置将粒子数据分配给每个子区域，并针对子区域中的每一个生成输入数据。管理节点200将子区域之一分配给计算节点100、100-1、100-2和100-3中的每一个，并将每个分配的子区域的输入数据发送至对应的计算节点。

在所有时间步完成后，管理节点200从计算节点100、100-1、100-2和100-3收集结果数据。结果数据包括每个时间步处粒子的位置。计算节点100、100-1、100-2和100-3中的每一个将结果数据发送至管理节点200，该结果数据包括在粒子各自存在于分配的子区域中期间粒子的位置。管理节点200可以将结果数据存储在非易失性存储装置诸如HDD中。另外，管理节点200可以将结果数据呈现在显示单元上。例如，管理节点200可以显示每个粒子由点表示并且粒子的位置的变化由点的运动表示的视频。另外，管理节点200可以将结果数据发送至不同的信息处理器。

注意，计算节点100、100-1、100-2和100-3以及管理节点200可以不被实现为独立的计算机，而是被实现为布置在包括大量处理器的并行处理单元中的不同处理器。另外，计算节点100、100-1、100-2和100-3之一可以承担管理节点200的处理。

图3是示出计算节点的示例性硬件配置的框图。

计算节点100包括CPU 101、RAM 102、HDD 103、图像接口104、输入装置接口105、介质读取器106和通信接口107。这些单元分别连接至总线。CPU 101对应于第一实施方式的处理单元12。RAM 102或HDD 103对应于第一实施方式的存储单元11。图3的计算节点100的相同硬件配置类似地应用于其它计算节点和管理节点200。

CPU 101是被配置成执行程序指令的处理器。CPU 101读出存储在HDD 103中的程序和数据的至少一部分，将它们加载到RAM 102中，并且执行加载的程序。注意，CPU 101可以包括两个或更多个处理器核，并且计算节点100可以包括两个或更多个处理器。术语“多处理器”或简单地“处理器”可以用于指代处理器的集合。

RAM 102是易失性半导体存储器，其用于在其中临时存储要由CPU 101执行的程序和要由CPU 101使用以用于其计算的数据。注意，信息处理器100可以设置有除了RAM之外的不同类型的存储器，或者可以设置有两个或更多个存储器装置。

HDD 103是非易失性存储装置，其用于在其中存储软件程序诸如操作系统(OS)、中间件和应用软件以及各种类型的数据。注意，计算节点100可以设置有不同类型的存储装置诸如闪存或固态驱动器(SSD)，或者可以设置有两个或更多个存储装置。

图像接口104根据来自CPU 101的绘制命令产生视频图像，并将它们显示在耦接至计算节点100的显示装置111的屏幕上。显示装置111可以是任何类型的显示器，例如阴极射线电子管(CRT)显示器；液晶显示器(LCD)；有机电致发光(OEL)显示器或投影仪。另外，除了显示装置111之外的输出装置诸如打印机也可以连接至计算节点100。

输入装置接口105从连接至计算节点100的输入装置112接收输入信号。各种类型的输入装置可以用作输入装置112，例如鼠标、触摸面板、触摸板或键盘。多种类型的输入装置可以连接至计算节点100。

介质读取器106是用于读取记录在存储介质113上的程序和数据的装置。各种类型的存储介质可以用作存储介质113，例如磁盘诸如软盘(FD)或HDD、光盘诸如致密盘(CD)或数字通用盘(DVD)和半导体存储器。介质读取器106将从存储介质113读出的程序和数据复制到不同的存储介质例如RAM 102或HDD 103。读取的程序例如由CPU 101执行。注意，存储介质113可以是便携式存储介质并且可以用于分发程序和数据。另外，有时将存储介质113和HDD 103称为计算机可读存储介质。

通信接口107连接至网络30，并且经由网络30与不同的计算节点和管理节点200通信。通信接口107是连接至有线通信装置诸如交换机或路由器的有线通信接口。

接下来描述DEM方案。

图4示出了DEM方案中的示例性粒子模型。

DEM方案涉及以下计算步骤：计算作用在目标粒子上的净力；基于计算出的净力和目标粒子的质量计算目标粒子的加速度；基于计算出的加速度和当前速度计算下一时间步的速度；基于计算出的下一时间步的速度和当前位置计算下一时间步的位置。作用在目标粒子上的净力被计算为外力例如重力和由与目标粒子接触的一个或更多个相邻粒子(接触粒子)所施加的接触力之和。因此，DEM方案使用表示粒子之间的接触力的模型。

在此考虑粒子41与42彼此接触的情况。粒子41的位置由粒子41的中心的坐标给出。类似地，粒子42的位置由粒子42的中心的坐标给出。粒子41的半径可以与粒子42的半径相同或者不同。为了简化模拟，可以假定所有粒子都具有相同的半径。当粒子41和42彼此强烈碰撞时，它们可能在交叠区域侵入(penetrate)彼此。鉴于此，如果粒子41的中心与粒子42的中心之间的距离小于或等于粒子41和42的半径之和，则粒子41和42被确定为彼此接触。

当粒子41和42接触时，由粒子42施加在粒子41上的接触力被分解为法向分量(normal component)和切向分量(tangential component)。法向接触力由并联布置的阻尼器(dashpot)43和弹簧45建模。切向接触力由并联布置的阻尼器44和弹簧46建模。阻尼器43和44中的每一个是用于提供沿正速度的相反方向的与速度成比例的粘滞阻力的元件。速度越高，粘滞阻力越大，速度越低，粘滞阻力越小。弹簧45和46中的每一个是用于提供沿位移的相反方向的与从原始弹簧长度的位移成比例的弹性阻力的元件。位移越大，弹性阻力越大，位移越小，弹性阻力越小。

法向接触力是阻尼器43的粘滞阻力和弹簧45的弹性阻力之和。因此，在时间步t处的法向接触力F

F

位移D

切向接触力是阻尼器44的粘滞阻力和弹簧46的弹性阻力之和。因此，在时间步t处的切向接触力F

F

位移D

然后，将上述法向接触力和切向接触力相结合以获得由粒子42施加在粒子41上的接触力。如果粒子41与又一粒子接触，则还计算由该不同的粒子施加在粒子41上的接触力。由一个或更多个接触粒子施加在粒子41上的接触力与外力(例如重力)结合，从而计算作用在粒子41上的净力。

如上所述，沿法线方向的位移D

另一方面，通过将时间步t-1与t之间的沿切线方向的行进距离与时间步t-1处的位移D

图5示出了模拟区域的示例性划分。

整个模拟区域例如被划分为子区域51、52、53和54。子区域51、52、53和54中的每一个是正方形。子区域51位于整个区域的左上方；子区域52位于整个区域的右上方；子区域53位于整个区域的左下方；子区域54位于整个区域的右下方。因此，子区域51在右侧与子区域52接界；在下侧与子区域53接界；在右下点与子区域54接界。子区域52在左侧与子区域51接界；在左下点与子区域53接界；在下侧与子区域54接界。子区域53在上侧与子区域51接界；在右上点与子区域52接界；在右侧与子区域54接界。子区域54在左上点与子区域51接界；在上侧与子区域52接界；在左侧与子区域53接界。

子区域51被分配给计算节点100。子区域52被分配给计算节点100-1。子区域53被分配给计算节点100-2。子区域54被分配给计算节点100-3。

子区域51包括边界区域51a、51b和51c，边界区域51a、51b和51c中的每一个是在距子区域51与另一子区域或其他子区域的边界预定距离以内的点的集合。边界区域51c是在距子区域51的下边界预定距离以内并且距右边界预定距离以内的点的集合。边界区域51a是除了边界区域51c的点以外、在距子区域51的下边界预定距离以内的点的集合。边界区域51b是除了边界区域51c的点以外、在距子区域51的右边界预定距离以内的点的集合。

子区域52包括边界区域52a、52b和52c，边界区域52a、52b和52c中的每一个是在距子区域51与另一子区域或其他子区域的边界预定距离以内的点的集合。边界区域52c是在距子区域52的下边界预定距离以内并且距左边界预定距离以内的点的集合。边界区域52a是除了边界区域52c的点以外、在距子区域52的下边界预定距离以内的点的集合。边界区域52b是除了边界区域52c的点以外、在距子区域52的左边界预定距离以内的点的集合。

子区域53包括边界区域53a、53b和53c，边界区域53a、53b和53c中的每一个是在距子区域53与另一子区域或其他子区域的边界预定距离以内的点的集合。边界区域53c是在距子区域53的上边界预定距离以内并且距右边界预定距离以内的点的集合。边界区域53a是除了边界区域53c的点以外、在距子区域53的上边界预定距离以内的点的集合。边界区域53b是除了边界区域53c的点以外、在距子区域53的右边界预定距离以内的点的集合。

子区域54包括边界区域54a、54b和54c，边界区域54a、54b和54c中的每一个是在距子区域54与另一子区域或其他子区域的边界预定距离以内的点的集合。边界区域54c是在距子区域54的上边界预定距离以内并且距左边界预定距离以内的点的集合。边界区域54a是除了边界区域54c的点以外、在距子区域54的上边界预定距离以内的点的集合。边界区域54b是除了边界区域54c的点以外、在距子区域54的左边界预定距离以内的点的集合。

如上所描述的，属于一个子区域的粒子可以跨两个子区域的边界与属于另一子区域的不同粒子接触。鉴于此，计算节点100、100-1、100-2和100-3共享属于边界区域的粒子的粒子数据。

具体地，计算节点100将与边界区域51b和51c相关联的粒子数据发送至计算节点100-1；将与边界区域51a和51c相关联的粒子数据发送至计算节点100-2；将与边界区域51c相关联的粒子数据发送至计算节点100-3。计算节点100-1将与边界区域52b和52c相关联的粒子数据发送至计算节点100；将与边界区域52c相关联的粒子数据发送至计算节点100-2；将与边界区域52a和52c相关联的粒子数据发送至计算节点100-3。计算节点100-2将与边界区域53a和53c相关联的粒子数据发送至计算节点100；将与边界区域53c相关联的粒子数据发送至计算节点100-1；将与边界区域53b和53c相关联的粒子数据发送至计算节点100-3。计算节点100-3将与边界区域54c相关联的粒子数据发送至计算节点100；将与边界区域54a和54c相关联的粒子数据发送至计算节点100-1；将与边界区域54b和54c相关联的粒子数据发送至计算节点100-2。

因此，负责子区域51的计算节点100除了属于子区域51的粒子的粒子数据以外，还具有属于边界区域52b、52c、53a、53c和54c的粒子的粒子数据。计算节点100针对属于子区域51的每个粒子确定接触粒子，然后计算位置和速度。从计算节点100的角度来看，以下将子区域51有时称为“主区域”，将边界区域52b、52c、53a、53c和54c一起有时统称为“裙摆区域(skirt region)”。

类似地，负责子区域52的计算节点100-1除了属于子区域52的粒子的粒子数据以外，还具有属于边界区域51b、51c、53c、54a和54c的粒子的粒子数据。计算节点100-1针对属于子区域52的每个粒子确定接触粒子，然后计算位置和速度。从计算节点100-1的角度来看，以下将子区域52有时称为“主区域”，将边界区域51b、51c、53c、54a和54c一起有时统称为“裙摆区域”。

负责子区域53的计算节点100-2除了属于子区域53的粒子的粒子数据以外，还具有属于边界区域51a、51c、52c、54b和54c的粒子的粒子数据。计算节点100-2针对属于子区域53的每个粒子确定接触粒子，然后计算位置和速度。从计算节点100-2的角度来看，以下将子区域53有时称为“主区域”，将边界区域51a、51c、52c、54b和54c一起有时统称为“裙摆区域”。

负责子区域54的计算节点100-3除了属于子区域54的粒子的粒子数据以外，还具有属于边界区域51c、52a、52c、53b和53c的粒子的粒子数据。计算节点100-3针对属于子区域54的每个粒子确定接触粒子，然后计算位置和速度。从计算节点100-3的角度来看，以下将子区域54有时称为“主区域”，将边界区域51c、52a、52c、53b和53c一起有时统称为“裙摆区域”。

另外，如上所描述的，当模拟前进一个时间步时，属于子区域的粒子可能会移动到不同的子区域。注意，每个粒子所属的子区域由粒子中心的坐标确定。计算节点100、100-1、100-2和100-3相互交换已经移出其被分配的子区域的粒子的粒子数据。

具体地，计算节点100将已经从右边界离开子区域51的粒子的粒子数据发送至计算节点100-1；将已经从下边界离开的粒子的粒子数据发送至计算节点100-2；并且将已经从右下点离开的粒子的粒子数据发送至计算节点100-3。计算节点100-1将已经从左边界离开子区域52的粒子的粒子数据发送至计算节点100；将已经从右下点离开的粒子的粒子数据发送至计算节点100-2；并且将已经从下边界离开的粒子的粒子数据发送至计算节点100-3。计算节点100-2将已经从上边界离开子区域53的粒子的粒子数据发送至计算节点100；将已经从右上点离开的粒子的粒子数据发送至计算节点100-1；并且将已经从右边界离开的粒子的粒子数据发送至计算节点100-3。计算节点100-3将已经从左上点离开子区域54的粒子的粒子数据发送至计算节点100；将已经从上边界离开的粒子的粒子数据发送至计算节点100-1；并且将已经从左边界离开的粒子的粒子数据发送至计算节点100-2。

图6是示出相邻粒子的示例性运动的第一图。

主区域61对应于图5的子区域51。裙摆区域62对应于图5的边界区域52b、52c、53a、53c和54c的集合。着重于属于主区域61的一个粒子，在目标粒子周围存在的相邻粒子根据它们中的每一个属于主区域61还是裙摆区域62以及是否与目标粒子接触而被分类为四种状态。然后，根据先前时间步处和当前时间步处的状态的组合将相邻粒子的运动分类为十六种模式。图6描绘了十六种模式中的四种模式。

粒子63是主区域61中的目标粒子。粒子63-1至63-4是在粒子63周围存在的相邻粒子。在先前时间步处，粒子63-1和63-3属于主区域61，并且与粒子63接触。粒子63-2和63-4属于裙摆区域62，并且与粒子63接触。在当前时间步处，粒子63-1属于主区域61，并且与粒子63接触。粒子63-2已经从裙摆区域62移动到主区域61中，并且与粒子63接触。粒子63-3已经从主区域61移动到裙摆区域62中，并且与粒子63接触。粒子63-4属于裙摆区域62，并且与粒子63接触。

以上描述的粒子63-1至63-4在先前时间步处和当前时间步处均与粒子63接触。因此，在计算作用在粒子63上的净力时，将粒子63-1至63-4中的每一个的沿切线方向的位移从先前时间步复制到当前时间步。

图7是示出相邻粒子的示例性运动的第二图。

图7描绘了前述十六个模式中的剩余十二个模式。粒子63-5至63-16是在粒子63周围存在的相邻粒子。

在先前时间步处，粒子63-5和63-11属于主区域61，并且与粒子63接触。粒子63-6、63-7、63-12和63-13属于主区域61，并且不与粒子63接触。粒子63-8和63-14属于裙摆区域62，并且与粒子63接触。粒子63-9、63-10、63-15和63-16属于裙摆区域62，并且不与粒子63接触。

在当前时间步处，粒子63-5属于主区域61，并且移动为不与粒子63接触。粒子63-6属于主区域61，并且与粒子63接触。粒子63-7属于主区域61，并且不与粒子63接触。粒子63-8已经从裙摆区域62移动到主区域61中，并且移动为不与粒子63接触。

粒子63-9已经从裙摆区域62移动到主区域61中，并且与粒子63接触。粒子63-10已经从裙摆区域62移动到主区域61中，并且不与粒子63接触。粒子63-11已经从主区域61移动到裙摆区域62中，并且移动为不与粒子63接触。粒子63-12已经从主区域61移动到裙摆区域62中，并且与粒子63接触。

粒子63-13已经从主区域62移动到裙摆区域61中，并且不与粒子63接触。粒子63-14属于裙摆区域62，并且移动为不与粒子63接触。粒子63-15属于裙摆区域62，并且与粒子63接触。粒子63-16属于裙摆区域62，并且不与粒子63接触。

以上描述的粒子63-5至63-16在先前时间步和当前时间步的至少之一处不与粒子63接触。因此，作用在粒子63上的净力的计算不涉及复制这些粒子63-5至63-16的沿切线方向的位移。

图8示出了沿切线方向的位移的示例性复制。

在这里，让我们看一种复制沿切线方向的位移的相对简单的方法。

对于主区域中的一个粒子(目标粒子)，生成先前时间步的接触粒子列表71和位移列表72。接触粒子列表71是列举在先前时间步处与目标粒子接触的接触粒子的粒子编号的数据。接触粒子各自属于主区域或裙摆区域。位移列表72是列举沿切线方向的位移的数据，其映射到在接触粒子列表71中登记的粒子编号。位移列表72中的每个位移是直到先前时间步的由相对应的粒子编号指示的接触粒子的位移。

当模拟前进一个时间步时，即，模拟时钟上的时间前进一个单位时间步时，针对上述相同的粒子生成当前时间步的接触粒子列表73和位移列表74。接触粒子列表73是列举在当前时间步处与目标粒子接触的接触粒子的粒子编号的数据。位移列表74是列举沿切线方向的位移的数据，其映射到在接触粒子列表73中登记的粒子编号。位移列表74中的每个位移是直到当前时间步的由相对应的粒子编号指示的接触粒子的位移。

在完成当前时间步的接触粒子列表73和位移列表74的生成时，可以删除先前时间步的接触粒子列表71和位移列表72。然而，注意，在存在登记在接触粒子列表71和73两者中的各个接触粒子的情况下，将与这些接触粒子相对应的位移从位移列表72复制到位移列表74。以这种方式，在当前时间步开始时完成直到先前时间步的位移的复制。

然后，向复制的位移中的每一个添加在先前时间步与当前时间步之间沿切线方向的相对行进距离(增量位移)，从而获得相对应的接触粒子的直到当前时间步的位移。例如，通过从目标粒子在当前时间步处的速度的切向分量与单位时间步的乘积中减去接触粒子在当前时间步处的速度的切向分量与单位时间步的乘积来计算沿切线方向的增量位移。

以下是一种复制沿切线方向的位移的方法。首先，与先前时间步的接触粒子列表71无关，通过针对属于主区域或裙摆区域的每个粒子检查其是否与目标粒子接触来检测当前时间步处的接触粒子。由此，生成当前时间步的接触粒子列表73。

随着接触粒子列表73的生成，针对在接触粒子列表73中登记的每个接触粒子扫描接触粒子列表71，以在接触粒子列表71中找到相同的接触粒子。当在接触粒子列表71中检测到相同的接触粒子时，从位移列表72内与在接触粒子列表71中检测到的接触粒子的位置相对应的位置读取位移。然后，将读取的位移复制到位移列表74内与在接触粒子列表73中登记的接触粒子的位置相对应的位置。另一方面，如果在接触粒子列表71中未包括相同的接触粒子，则将位移0写入位移列表74内与在接触粒子列表73中登记的接触粒子的位置相对应的位置。

例如，当接触粒子列表73的第一元素表示与接触粒子列表71的第一元素相同的接触粒子时，将位移列表72的第一元素复制成为位移列表74的第一元素。另外，当接触粒子列表73的第二元素表示与接触粒子列表71的第三元素相同的接触粒子时，将位移列表72的第三元素复制成为位移列表74的第二元素。如果接触粒子列表71不包括表示与接触粒子列表73的第三元素和第四元素相同的接触粒子的元素，则位移列表74的第三元素和第四元素的位移被设置为0。

在以上描述的位移复制方法中，在接触粒子列表71与73之间比较粒子编号的比较任务的计算成本是O(N

注意，包括在接触粒子列表71和73中的粒子编号是一个计算节点上的局部唯一标识编号。例如，计算节点100独立于计算节点100-1、100-2和100-3来将粒子编号分配给属于计算节点100的主区域和裙摆区域的粒子。当粒子在子区域之间移入和移出时，不同的粒子编号可能被分配给在先前时间步和当前时间步中相同的粒子。然而，注意，计算节点100针对保留在主区域中的每个粒子包含先前时间步和当前时间步的粒子编号之间的映射。另外，根据第二实施方式，除了粒子编号之外，全局标识符(ID)也用作粒子的标识编号。全局ID是针对所有计算节点的全局唯一标识编号。

接下来描述计算节点100的功能。

图9是示出计算节点的示例性功能的框图。

计算节点100包括模型存储单元121、中间数据存储单元122、结果数据存储单元123、数据输入和输出单元124、时间步执行单元125以及处理间通信单元126。模型存储单元121、中间数据存储单元122和结果数据存储单元123使用例如固定在RAM 102或HDD 103中的存储区域来实现。数据输入和输出单元124、时间步执行单元125和处理间通信单元126例如使用由CPU 101执行的程序来实现。注意，计算节点100-1、100-2和100-3还各自具有与在计算节点100中发现的模块相同的模块。

模型存储单元121在其中存储表示分配给计算节点100的子区域的模型的输入数据。输入数据包括区域数据诸如子区域的大小和坐标以及壁的位置。另外，输入数据包括最初在子区域中分布的粒子的位置。输入数据还包括物理参数例如重力、粒子半径以及接触模型的弹性模量和粘滞阻尼系数。此外，输入数据包括控制参数例如单位时间步持续时间和最大时间步编号。

中间数据存储单元122在其中存储在模拟期间生成的中间数据。中间数据包括指示属于主区域的粒子的位置和速度的粒子数据，以及指示属于裙摆区域的粒子的位置和速度的粒子数据。另外，中间数据包括先前时间步的粒子编号与当前时间步的粒子编号之间的信息映射。中间数据还包括属于主区域的每个粒子的接触粒子列表和位移列表。

结果数据存储单元123在其中存储指示分配给计算节点100的子区域的模拟结果的结果数据。结果数据包括表示针对每个时间步存在于子区域中的粒子及其在该时间步处的位置的信息。按时间顺序跟踪每个粒子的位置使得能够对粉末粒子的运动进行可视化。

数据输入和输出单元124与控制并行处理的管理节点200通信。具体地，数据输入和输出单元124从管理节点200接收分配给计算节点100的子区域的输入数据，然后将其存储在模型存储单元121中。另外，当模拟时间达到最大时间步编号时，数据输入和输出单元124从结果数据存储单元123读取结果数据，并将其发送至管理节点200。然后，管理节点200将多个子区域的结果数据进行聚合，并对整个区域上的粉末粒子的运动进行可视化。

时间步执行单元125使用存储在模型存储单元121中的输入数据沿时间轴逐个时间步地计算分配的子区域中的粒子的运动。将模拟前进一个时间步对应于将模拟时钟上的时间前进预定的非常短的时间段。每个时间步包括对属于子区域的每个粒子的以下计算过程：检查接触粒子并计算由接触粒子施加的接触力；通过将接触力和重力结合来计算净力；以及基于当前位置和速度以及净力计算下一时间步的位置和速度。时间步执行单元125相应地更新存储在中间数据存储单元122中的中间数据。此外，时间步执行单元125将每个时间步处粒子的位置写入到结果数据存储单元123中。

处理间通信单元126实现由计算节点100的时间步执行单元125执行的处理与由计算节点100-1、100-2和100-3各自执行的处理之间的处理间通信(interprocesscommunication)。具体地，处理间通信单元126从中间数据存储单元122读取一些中间数据，并将其发送至计算节点100-1、100-2和100-3。另外，处理间通信单元126从计算节点100-1、100-2和100-3接收中间数据，并将其存储在中间数据存储单元122中。

要发送至计算节点100-1、100-2和100-3的中间数据包括已经从分配给计算节点100的子区域移出的粒子的粒子数据。此外，要发送的中间数据包括与计算节点100-1、100-2和100-3的裙摆区域相对应的边界区域中的粒子的粒子数据。另一方面，从计算节点100-1、100-2和100-3接收的中间数据包括已经移动到分配给计算节点100的子区域中的粒子的粒子数据。此外，所接收的中间数据包括与计算节点100的裙摆区域相对应的边界区域中的粒子的粒子数据。

图10示出了示例性结果表。

结果表131被存储在结果数据存储单元123中。结果表131包括名称为时间步编号；全局ID；和位置的列。时间步编号下的行包含一个或更多个用于标识时间步的整数。在时间步编号下的条目“1”对应于模拟开始时的初始状态。全局ID下的行包含全局唯一标识编号，全局唯一标识编号中的每一个在模拟区域中的所有粒子中唯一地标识粒子。位置下的行包含表示相对应的粒子的位置的坐标。注意，提供的输入数据包括与时间步编号下的条目“1”相对应的全局ID和位置下的条目，作为关于粒子的初始分布的信息。

图11示出了示例性粒子表和粒子编号映射表。

主粒子表132存储在中间数据存储单元122中。主粒子表132包含属于主区域的粒子的最新粒子数据。主粒子表132包括名称为全局ID；位置；和速度的列。全局ID下的行包含分配给属于主区域的粒子的全局ID。位置下的行包含表示属于主区域的粒子的位置的坐标。速度下的行包含属于主区域的粒子的速度向量。当模拟区域是二维时，坐标和速度向量是二维向量。

指示主粒子表132中的记录的顺序的索引被用作属于主区域的粒子的粒子编号。例如，由主粒子表132的第一记录指示的粒子具有粒子编号1。当粒子移出主区域时，从主粒子表132中删除与该粒子相对应的记录。当粒子移入到主区域中时，将与该粒子相对应的记录插入到主粒子表132中。因此，当模拟前进一个时间步时，分配给现有记录的索引可能会更改，从而导致现有粒子的粒子编号发生变化。

注意，主粒子表132可以被组织为多维阵列。可替选地，主粒子表132可以被组织为包括全局ID阵列、位置阵列和速度阵列的三个一维阵列。索引对应于一个或多个阵列中的下标。当主粒子表132被组织为三个一维阵列时，由相同索引表示的全局ID、位置和速度彼此映射作为关于相同粒子的信息。

裙摆粒子表133存储在中间数据存储单元122中。裙摆粒子表133包含属于裙摆区域的粒子的最新粒子数据。裙摆粒子表133包括名称为全局ID；位置；和速度的列。全局ID下的行包含分配给属于裙摆区域的粒子的全局ID。位置下的行包含表示属于裙摆区域的粒子的位置的坐标。速度下的行包含属于裙摆区域的粒子的速度向量。

向指示裙摆粒子表133中的记录的顺序的索引的每一个中添加预定的前缀，得到的索引用作属于裙摆区域的粒子的粒子编号。给出前缀使得得到的粒子编号与属于主区域的粒子的粒子编号不冲突。例如，由裙摆粒子表133中的第一记录指示的粒子具有粒子编号900001。裙摆粒子表133的粒子数据是从计算节点100-1、100-2和100-3接收的粒子数据的汇编。

注意，裙摆粒子表133可以被组织为多维阵列。可替选地，裙摆粒子表133可以被组织为包括全局ID阵列、位置阵列和速度阵列的三个一维阵列。索引对应于一个或多个阵列中的下标。当裙摆粒子表133被组织为三个一维阵列时，由相同索引表示的全局ID、位置和速度彼此映射作为关于相同粒子的信息。

主粒子表132中的位置和速度由计算节点100进行更新。此外，主粒子表132中登记的记录可以保留用于下一时间步。另一方面，由计算节点100-1、100-2和100-3计算的裙摆粒子表133中的位置和速度不被计算节点100更新。此外，在裙摆粒子表133中登记的记录在每个时间步被重置，并且不保留用于下一时间步。计算节点100可以在每个时间步通过将从计算节点100-1、100-2和100-3接收的裙摆区域的粒子数据放在一起以生成新的裙摆粒子表133。

粒子编号映射表134被存储在中间数据存储单元122中。粒子编号映射表134针对在先前时间步处属于主区域并且在当前时间步处仍保持在主区域中的每个粒子将先前时间步的粒子编号(称为“先前粒子编号”)映射到当前时间步的粒子编号。粒子编号映射表134是先前粒子编号的一维阵列。粒子编号映射表134的索引对应于当前时间步的粒子编号。

一个索引和由该索引表示的先前粒子编号分别表示同一粒子的当前时间步的粒子编号和先前时间步的粒子编号。例如，当粒子编号映射表134的第二记录指示先前粒子编号为3时，先前时间步的粒子编号为3的粒子与当前时间步的粒子编号为2的粒子相同。可以形成粒子编号映射表134例如以使得能够基于先前粒子编号快速检索当前时间步的粒子编号。

注意，当在主粒子表132中反映已经移入和移出主区域的粒子时，计算节点100生成粒子编号映射表134。因为计算节点100处理主粒子表132中的记录的插入和删除，所以它包含保留在主区域中的粒子的粒子编号的变化。另一方面，因为计算节点100在每个时间步都重置裙摆粒子表133，所以它未掌握保留在裙摆区域中的粒子的粒子编号的变化以及已经在主区域与裙摆区域之间移动的粒子的粒子编号的变化。

图12示出了示例性接触粒子列表。

根据第二实施方式，计算节点100针对属于主区域的每个粒子生成接触粒子列表、全局ID列表和位移列表。直到当前时间步的接触粒子列表、全局ID列表和位移列表的生成完成，先前时间步的接触粒子列表、全局ID列表和位移列表暂时共存。此外，根据第二实施方式，计算节点100将与目标粒子接触的接触粒子分类为属于主区域的接触粒子和属于裙摆区域的接触粒子。然后，计算节点100将接触粒子列表、全局ID列表和位移列表中的每一个划分为与主区域相关联的列表和与裙摆区域相关联的列表。

具体地，中间数据存储单元122针对先前时间步存储主接触粒子列表141、全局ID列表142、位移列表143、裙摆接触粒子列表144、全局ID列表145和位移列表146。此外，中间数据存储单元122针对当前时间步存储主接触粒子列表151、全局ID列表152、位移列表153、裙摆接触粒子列表154、全局ID列表155和位移列表156。

主接触粒子列表141是列举在先前时间步处与目标粒子接触的接触粒子中属于主区域的接触粒子的粒子编号的数据。全局ID列表142是列举在主接触粒子列表141中登记的接触粒子的全局ID的数据。位移列表143是列举与目标粒子和在主接触粒子列表141中登记的各个接触粒子相关联的、直到先前时间步的沿切线方向的位移的数据。所有主接触粒子列表141、全局ID列表142和位移列表143具有相同的长度。主接触粒子列表141中的第n个粒子编号(n是大于或等于1的整数)、全局ID列表142中的第n个全局ID和位移列表143中的第n个位移是关于同一接触粒子的信息。

裙摆接触粒子列表144是列举在先前时间步处与目标粒子接触的接触粒子中属于裙摆区域的接触粒子的粒子编号的数据。全局ID列表145是列举在裙摆接触粒子列表144中登记的接触粒子的全局ID的数据。位移列表146是列举与目标粒子和在裙摆接触粒子列表144中登记的各个接触粒子相关联的、直到先前时间步的沿切线方向的位移的数据。所有裙摆接触粒子列表144、全局ID列表145和位移列表146具有相同的长度。裙摆接触粒子列表144中的第n个粒子编号、全局ID列表145中的第n个全局ID和位移列表146中的第n个位移是关于同一接触粒子的信息。

主接触粒子列表151是列举在当前时间步处与目标粒子接触的接触粒子中属于主区域的接触粒子的粒子编号的数据。全局ID列表152是列举在主接触粒子列表151中登记的接触粒子的全局ID的数据。位移列表153是列举与目标粒子和在主接触粒子列表151中登记的各个接触粒子相关联的、直到当前时间步的沿切线方向的位移的数据。所有主接触粒子列表151、全局ID列表152和位移列表153具有相同的长度。主接触粒子列表151中的第n个粒子编号、全局ID列表152中的第n个全局ID和位移列表153中的第n个位移是关于同一接触粒子的信息。

裙摆接触粒子列表154是列举在当前时间步处与目标粒子接触的接触粒子中属于裙摆区域的接触粒子的粒子编号的数据。全局ID列表155是列举在裙摆接触粒子列表154中登记的接触粒子的全局ID的数据。位移列表156是列举与目标粒子和在裙摆接触粒子列表154中登记的各个接触粒子相关联的、直到当前时间步的沿切线方向的位移的数据。所有裙摆接触粒子列表154、全局ID列表155和位移列表156具有相同的长度。裙摆接触粒子列表154中的第n个粒子编号、全局ID列表155中的第n个全局ID和位移列表156中的第n个位移是关于同一接触粒子的信息。

注意，针对在属于主区域的粒子中靠近裙摆区域的粒子生成的裙摆接触粒子列表144和154、全局ID列表145和155以及位移列表146和156可以各自包含一个或更多个元素。另一方面，针对远离裙摆区域的粒子生成的裙摆接触粒子列表144和154、全局ID列表145和155以及位移列表146和156期望为空，不包含任何元素。因此，针对大量粒子生成的裙摆接触粒子列表144和154、全局ID列表145和155以及位移列表146和156为空。

接下来描述计算节点100的处理过程。

图13是示出DEM方案的示例性过程的流程图。

(步骤S10)数据输入输出单元124从管理节点200接收计算节点100负责的子区域的输入数据。时间步执行单元125基于输入数据对主粒子表132进行初始化。经初始化的主粒子表132登记输入数据中包括的每个粒子的全局ID和位置。然后，在主粒子表132中，将每个粒子的速度设置为0。然而，注意，如果输入数据包括速度的初始值，则代替地，将这些值登记在主粒子表132中。

(步骤S11)处理间通信单元126在主粒子表132中搜索属于主区域的粒子中的位于边界区域中的粒子的粒子数据。边界区域是在主区域中的距裙摆区域预定距离内的区域。根据输入数据中包括的子区域的坐标范围和主粒子表132中包括的粒子的位置的坐标来确定粒子是属于一个边界区域还是属于另一边界区域。一个粒子的粒子数据对应于主粒子表132中的一个记录，并且包括粒子的全局ID、位置和速度。处理间通信单元126复制在主粒子表132中找到的边界区域中的粒子的粒子数据，并将该粒子数据发送至计算节点100-1、100-2和100-3。

(步骤S12)处理间通信单元126从计算节点100-1、100-2、100-3接收边界区域中的粒子的粒子数据。这里接收到的粒子数据是属于计算节点100的裙摆区域的粒子的粒子数据。处理间通信单元126将接收到的粒子数据登记在裙摆粒子表133中。一个粒子的粒子数据对应于裙摆粒子表133中的一个记录，并且包括粒子的全局ID、位置和速度。注意，如果先前时间步的粒子数据保留在裙摆粒子表133中，则该粒子数据被丢弃。

(步骤S13)时间步执行单元125针对属于主区域的每个粒子生成接触粒子列表、全局ID列表和位移列表。具体地，时间步执行单元125生成当前时间步的以下列表：主接触粒子列表151；裙摆接触粒子列表154；全局ID列表152和155；以及位移列表153和156。时间步执行单元125在完成复制处理之后丢弃先前时间步的以下列表：主接触粒子列表141；裙摆接触粒子列表144；全局ID列表142和145；以及位移列表143和146。复制处理可以涉及将位移从先前时间步的位移列表143和146复制到当前时间步的位移列表153和156。接触粒子列表的生成将在后面详细描述。

(步骤S14)时间步执行单元125更新属于主区域的每个粒子的位移列表153和156中包括的位移。在位移列表153的更新中，时间步执行单元125从主粒子表132读取每个接触粒子的速度，然后计算目标粒子相对于接触粒子的相对速度，并计算相对速度的切向分量。时间步执行单元125基于相对速度的切向分量计算从先前时间步到当前时间步的沿切线方向的增量位移，并将所计算的增量位移与位移列表153中的相对应的位移相加。类似地，在位移列表156的更新中，时间步执行单元125从裙摆粒子表133读取每个接触粒子的速度，然后计算从先前时间步到当前时间步的沿切线方向的增量位移，并将计算出的增量位移与位移列表156中的相对应的位移相加。注意，位移列表153和156的更新可以在步骤S13中执行。

(步骤S15)时间步执行单元125计算作用在属于主区域的每个粒子上的净力。在计算净力时，时间步执行单元125从主粒子表132中读取与在主接触粒子列表151中登记的每个粒子编号相关联的粒子数据，从裙摆粒子表133中读取与在裙摆接触粒子列表154中登记的每个粒子编号相关联的粒子数据。对于每个接触粒子，时间步执行单元125根据上述接触模型基于目标粒子的位置和速度、接触粒子的位置和速度、接触粒子沿切线方向的位移以及包括在输入数据中的弹性模量和粘滞阻尼系数计算接触力。时间步执行单元125将由一个或更多个接触粒子施加在目标粒子上的接触力与外力诸如重力相结合，从而计算作用在目标粒子上的净力。

(步骤S16)对于属于主区域的每个粒子，时间步执行单元125基于在步骤S15中计算的净力和粒子的质量来计算加速度，然后基于计算的加速度以及当前位置和速度来计算之后的一个单位时间步的位置和速度。时间步执行单元125更新在主粒子表132中登记的每个粒子的位置和速度。此外，时间步执行单元125将针对下一时间步计算出的每个粒子的位置添加到结果表131中。

(步骤S17)处理间通信单元126从主粒子表132中提取已经从分配给计算节点100的子区域(主区域)中移出的粒子的粒子数据。根据输入数据中包括的子区域的坐标范围和主粒子表132中包括的粒子的位置的坐标来确定粒子是否已经移出子区域。处理间通信单元126将从主粒子表132提取的粒子数据发送至计算节点100-1、100-2和100-3。要发送的粒子数据包括全局ID、位置和速度。从主粒子表132中删除提取的粒子数据。注意，当从主粒子表132中删除一个或更多个记录时，剩余记录可以被移位以填充由被删除的记录留下的间隙，以便没有空隙残留。在这种情况下，剩余记录的索引可能会改变。

(步骤S18)处理间通信单元126从计算节点100-1、100-2和100-3接收已经移动到分配给计算节点100的子区域(主区域)中的粒子的粒子数据。所接收的粒子数据包括全局ID、位置和速度。处理间通信单元126将接收到的粒子数据插入到主粒子表132中。例如，新记录被添加到主粒子表132的末尾。

(步骤S19)处理间通信单元126监测在步骤S17和S18中主粒子表132中的记录的删除和插入，并跟踪保留的未被删除的现有记录的索引的改变。现有记录表示在当前时间步处属于主区域并在下一时间步处保留在主区域中的粒子。处理间通信单元126生成使当前时间步的粒子编号与下一时间步的粒子编号相关的新粒子编号映射表134。可以丢弃旧粒子编号映射表134。

(步骤S20)时间步执行单元125确定是否结束时间步的迭代。例如，时间步执行单元125确定时间步的编号是否已经达到输入数据中指定的最大时间步编号。如果时间步执行单元125确定结束时间步的迭代，则处理进入步骤S21。如果不是，则处理返回到步骤S11。

(步骤S21)数据输入和输出单元124从结果数据存储单元123读取结果表131，并将结果表131的结果数据发送至管理节点200。

图14是示出接触粒子列表的生成的示例性过程的流程图的第一部分。

在以上描述的步骤S13中执行接触粒子列表的生成。针对属于主区域的每个粒子执行接触粒子列表的生成。

接下来的步骤S30至S37包括复制在先前时间步处于接触同时属于主区域然后在当前时间步处于接触同时属于主区域的每个接触粒子(下文中有时称为“主-主粒子”)的位移。该接触粒子对应于前述粒子63-1。

(步骤S30)时间步执行单元125在先前时间步的主接触粒子列表141中登记的先前粒子编号中选择一个先前粒子编号。

(步骤S31)时间步执行单元125在粒子编号映射表134中搜索与在步骤S30中选择的先前粒子编号相对应的当前时间步的粒子编号。

(步骤S32)在步骤S31中，时间步执行单元125确定在粒子编号映射表134中是否找到适当的粒子编号。如果找到了适当的粒子编号，则处理进入步骤S33。如果没有，则处理进入步骤S37。

(步骤S33)时间步执行单元125从主粒子表132读取与在步骤S31中找到的粒子编号相对应的位置。时间步执行单元125基于目标粒子的位置和由找到的粒子编号指示的粒子的位置来检查关注的粒子是否与目标粒子接触。如果目标粒子与关注的粒子之间的距离小于或等于两个粒子的半径之和，则确定它们彼此接触。如果它们之间的距离大于两个粒子的半径之和，则确定它们彼此不接触。注意，预先从主粒子表132读取目标粒子的位置。

(步骤S34)时间步执行单元125确定两个粒子是否彼此接触。如果它们彼此接触，则处理进入步骤S35。如果没有，则处理进入步骤S37。

(步骤S35)时间步执行单元125将关注的粒子的粒子编号添加到当前时间步的主接触粒子列表151的末尾。另外，时间步执行单元125将关注的粒子的全局ID添加到当前时间步的全局ID列表152的末尾。可以通过从先前时间步的全局ID列表142中读取与在步骤S30中选择的先前粒子编号相对应的全局ID或者从主粒子表132中读取与当前时间步的粒子编号相对应的全局ID来获得全局ID。

(步骤S36)时间步执行单元125从先前时间步的位移列表143中读取与在步骤S30中选择的先前粒子编号相对应的位移。可使用与主接触粒子列表141的先前粒子编号相同的索引访问位移列表143中的期望位移。时间步执行单元125将读取的位移添加到当前时间步的位移列表153的末尾。

(步骤S37)时间步执行单元125确定在步骤S30中是否已经选择了主接触粒子列表141的所有先前粒子编号。如果已经选择了所有先前粒子编号，则处理进入步骤S38。如果没有，则处理返回到步骤S30。

图15从图14继续，是示出了接触粒子列表的生成的示例性过程的流程图的第二部分。

接下来的步骤S38至S46包括复制在先前时间步处于接触同时属于裙摆区域然后在当前时间步处于接触同时属于主区域的每个接触粒子(下文中有时称为“裙摆-主粒子”)的位移。该接触粒子对应于前述粒子63-2。

(步骤S38)时间步执行单元125从主粒子表132中选择一个不是目标粒子的粒子编号且尚未进行步骤S30至S37的粒子编号。这里可以使用标记来指示各个粒子编号是否已经进行了步骤S30至S37。这里由选择的粒子编号指示的粒子是属于主区域的尚未检查与目标粒子的接触的粒子。

(步骤S39)时间步执行单元125从主粒子表132读取与在步骤S38中选择的粒子编号相对应的位置。时间步执行单元125基于目标粒子的位置和由粒子编号指示的粒子的位置来检查关注的粒子是否与目标粒子接触。

(步骤S40)时间步执行单元125确定两个粒子是否彼此接触。如果它们彼此接触，则处理进入步骤S41。如果没有，则处理进入步骤S46。

(步骤S41)时间步执行单元125将关注的粒子的粒子编号添加到当前时间步的主接触粒子列表151的末尾。另外，时间步执行单元125从主粒子表132中读取与关注的粒子的粒子编号相对应的全局ID，并将关注的粒子的全局ID添加到当前时间步的全局ID列表152的末尾。

(步骤S42)时间步执行单元125从上到下扫描先前时间步的全局ID列表145，以搜索关注的粒子的全局ID。全局ID列表145指示在先前时间步属于裙摆区域的粒子。

(步骤S43)时间步执行单元125确定在步骤S42中是否找到适当的全局ID。如果找到适当的全局ID，则处理进入步骤S45。如果没有，则处理进入步骤S44。

(步骤S44)时间步执行单元125将指示该位移为0的数值添加到当前时间步的位移列表153的末尾。然后，处理进入步骤S46。

(步骤S45)时间步执行单元125从先前时间步的位移列表146中读取与在步骤S42中检测到的全局ID相对应的位移。可使用与全局ID列表145的全局ID相同的索引来访问位移列表146中的期望位移。时间步执行单元125将读取的位移添加到当前时间步的位移列表153的末尾。

(步骤S46)时间步执行单元125确定在步骤S38中是否已经选择了主粒子表132的所有粒子编号。如果已经选择了所有粒子编号，则处理进入步骤S47。如果没有，则处理返回到步骤S38。

图16从图15继续，是示出接触粒子列表的生成的示例性过程的流程图的第三部分。

接下来的步骤S47至S57包括复制在先前时间步处于接触同时属于主区域然后在当前时间步处于接触同时属于裙摆区域的每个接触粒子(下文中有时称为“主-裙摆粒子”)的位移。该接触粒子对应于前述粒子63-3。此外，步骤S47至S57包括复制在先前时间步处于接触同时属于裙摆区域以及在当前时间步处于接触同时属于裙摆区域的每个接触粒子(下文中有时称为“裙摆-裙摆粒子”)的位移。该接触粒子对应于前述粒子63-4。

(步骤S47)时间步执行单元125参考裙摆粒子表133从在当前时间步属于裙摆区域的粒子的粒子编号中选择一个粒子编号。通过将预定的前缀添加到裙摆粒子表133的索引来形成粒子编号。

(步骤S48)时间步执行单元125从裙摆粒子表133读取与在步骤S47中选择的粒子编号相对应的位置。时间步执行单元125基于目标粒子的位置和由所选择的粒子编号指示的粒子的位置来检查关注的粒子是否与目标粒子接触。

(步骤S49)时间步执行单元125确定两个粒子是否彼此接触。如果它们彼此接触，则处理进入步骤S50。如果没有，则处理进入步骤S57。

(步骤S50)时间步执行单元125将关注的粒子的粒子编号添加到当前时间步的裙摆接触粒子列表154的末尾。另外，时间步执行单元125从裙摆粒子表133中读取与关注的粒子的粒子编号相对应的全局ID，并将关注的粒子的全局ID添加到当前时间步的全局ID列表155的末尾。

(步骤S51)时间步执行单元125从上到下扫描先前时间步的全局ID列表142，以搜索关注的粒子的全局ID。全局ID列表142指示在先前时间步属于主区域的粒子。

(步骤S52)时间步执行单元125确定在步骤S51中是否找到适当的全局ID。如果找到适当的全局ID，则处理进入步骤S56。如果没有，则处理进入步骤S53。

(步骤S53)时间步执行单元125从上到下扫描先前时间步的全局ID列表145，以搜索关注的粒子的全局ID。全局ID列表145指示在先前时间步属于裙摆区域的粒子。

(步骤S54)时间步执行单元125确定在步骤S53中是否找到适当的全局ID。如果找到适当的全局ID，则处理进入步骤S56。如果没有，则处理进入步骤S55。

(步骤S55)时间步执行单元125将指示位移为0的数值添加到当前时间步的位移列表156的末尾。然后，处理进入步骤S57。

(步骤S56)时间步执行单元125从先前时间步的位移列表143中读取与在步骤S51中检测到的全局ID相对应的位移。可使用与全局ID列表142的全局ID相同的索引来访问位移列表143中的期望位移。可替选地，时间步执行单元125从先前时间步的位移列表146中读取与在步骤S53中检测到的全局ID相对应的位移。可使用与全局ID列表145的全局ID相同的索引来访问位移列表146中的期望位移。时间步执行单元125将读取的位移添加到当前时间步的位移列表156的末尾。

(步骤S57)时间步执行单元125确定在步骤S47中是否已经选择了裙摆粒子表133的所有粒子编号。如果已经选择了所有粒子编号，则接触粒子列表的生成结束。如果没有，则处理返回到步骤S47。

根据第二实施方式的信息处理系统，在DEM方案的粒子模拟中，针对属于主区域的每个粒子生成列举接触粒子的接触粒子列表和列举接触粒子的沿切线方向的位移的位移列表。因此，可以有效地计算由接触粒子施加在属于主区域的每个粒子上的接触力。另外，将接触粒子列表划分为指示属于主区域的接触粒子的主接触粒子列表和指示属于裙摆区域的接触粒子的裙摆接触粒子列表。这使得能够根据各个接触粒子所属区域来有效地管理关于各个接触粒子的信息。

另外，在生成当前时间步的主接触粒子列表时，优选地针对在先前时间步的主接触粒子列表中登记的粒子检查与目标粒子的接触。随后，针对未在先前时间步的主接触粒子列表中登记的粒子检查与目标粒子的接触。因此，对于主-主粒子，可以有效地将其位移从先前时间步的位移列表复制到当前时间步的位移列表。

例如，从主粒子表中检测在当前时间步属于主区域的每个接触粒子并且从先前时间步的主接触粒子列表中搜索接触粒子涉及对先前时间步的主接触粒子列表的重复扫描。计算复杂度为O(N

另外，对于属于主区域的每个粒子，除了接触粒子列表和位移列表之外，还生成列举接触粒子的全局ID的全局ID列表。这有利于针对裙摆-主粒子、主-裙摆粒子和裙摆-裙摆粒子从先前时间步的主接触粒子列表和裙摆接触粒子列表中有效搜索期望的接触粒子。作为结果，尽管每个粒子的粒子编号在先前时间步与当前时间步之间可能不同，但是对于除了主-主粒子之外的接触粒子也实现了有效的位移复制。

根据一个方面，可以降低粒子模拟的计算复杂度。