粒子模拟装置、粒子模拟方法以及粒子模拟程序与流程

本发明涉及模拟粒子和构造物的动态的粒子模拟装置、粒子模拟方法以及粒子模拟程序。

背景技术：

以往，一直采用DEM(Discrete Element Method：离散元方法)，SPH(Smoothed Particle Hydrodynamics：平滑粒子流体力学)以及MPS(Moving Particle Semi-implicit method：移动粒子半隐式法)等粒子法模拟来进行对粒子和复杂形状的构造物的动态的解析。在这样的模拟中，需要进行粒子与构造物的接触判定。最近，考虑到计算负荷、模拟的正确性等，提出使用符号距离函数(Signed Distance Function)的接触判定(例如，参见非专利文献1)。符号距离函数是使用任意尺寸的正交格子对成为模拟对象的空间内进行分割，将从格子点上到构造物表面为止的最接近距离存储于格子点(例如，参见非专利文献2)。从粒子到构造物为止的距离使用粒子所在的格子和相邻格子的距离函数数据实施内插插值而求取，并基于所计算出的距离来实施接触判定。

现有技术文献

非专利文献

非专利文献1：Y.Shigeto,M.Sakai,Arbitrary-shaped wall boundary modeling based on signed distance functions for granular flow simulations,Chem.Eng.J.231(2013)464-476。

非专利文献2：Y.Wu,J.Man,Z.Xie,A double layer method for constructing signed distance fields from triangle meshes,Graphical Models76(2014)214-223。

技术实现要素：

发明要解决的问题

针对与构造物的时间变化相伴的移动，例如如以下任一项那样处理距离函数。(1)将距离函数固定，使粒子的坐标系与构造物的移动和旋转对应地平行移动和旋转。但是，在该方法中，在存在多个构造物的情况下，即，在存在多个移动轴、旋转轴的情况下，无法使用。(2)每次构造物移动时，都重新创建距离函数。但是，在该方法中，用于重新创建距离函数的计算成本庞大，不易在实际计算中使用。(3)进行使用微分方程式的距离函数的移流计算。但是，在该方法中，由于是格子法所以会产生数值扩散而导致距离函数钝化，因此通常需要以一定时间间隔为单位重新计算距离函数的重新初始化。另外，此时，在此之前在格子点上顺利地变化着的距离函数会不连续地变化，产生不连续地变化的距离振动。因上述数值扩散、数值振动，模拟的准确性降低，在工业、学术领域中，出于计算精度方面的考虑而无法投入实际应用。(4)根据构造物的移动，而事先准备多个距离函数的时序数据。但是，该方法需要格子数×时序数的数据量，预先存储的数据量庞大。

本发明是鉴于上述问题点而完成的，目的在于提供能够以较少的计算负荷准确进行粒子与构造物之间的距离计算，并由此适当地实施粒子与构造物的模拟的粒子模拟装置、粒子模拟方法以及粒子模拟程序。

用于解决问题的手段

为了达成上述目的，本发明的一个实施方式的粒子模拟装置是基于作业空间内的粒子与构造物之间的相互作用力计算该粒子的位置和速度来模拟该粒子和该构造物的动态的粒子模拟装置，其中，具备：构造物初始信息取得单元，取得表示构造物的初始位置和形状的构造物初始信息；假想区域设定单元，基于通过由构造物初始信息取得单元取得的构造物初始信息表示的构造物的初始位置和形状，来在所述构造物的表面的附近设定多个假想区域，并关于所设定的多个假想区域中的各假想区域，取得到构造物的距离；位置信息取得单元，取得表示粒子的位置的位置信息；假想区域确定单元，根据通过由位置信息取得单元取得的位置信息表示的粒子的位置，确定预先设定的距离内的假想区域；粒子距离计算单元，关于由假想区域确定单元确定的假想区域，基于由假想区域设定单元计算出的距离，来计算从粒子到构造物的距离；相互作用力计算单元，基于由粒子距离计算单元计算出的距离，来计算粒子与构造物之间的相互作用力；以及时间迁移单元，基于由相互作用力计算单元计算出的相互作用力，计算下一时间步骤的粒子的位置和速度，并且在维持与所述构造物之间的位置关系的状态下，使假想区域移动。

在本发明的一个实施方式的粒子模拟装置中，基于关于在构造物的表面的附近设定的多个假想区域计算出的距离，进行粒子与构造物之间的距离计算。在模拟中，该假想区域在维持与构造物之间的位置关系的状态下移动。因此，无需重新计算关于假想区域的距离(距离函数)。另外，与以往的距离函数不同，只要针对在构造物的附近设定的假想区域计算距离即可。另外，假想区域针对每个构造物进行设定，因此即使在存在多个构造物的情况下、即存在多个相互不同的移动轴、旋转轴的情况下，也能够使用。另外，由于与距离建立了关联的假想区域移动，因此也不会产生因移流计算所致的数值扩散。即，根据本发明的一个实施方式的粒子模拟装置，能够以较少的计算负荷准确地进行粒子与构造物之间的距离计算，由此能够适当地进行粒子与构造物的模拟。

可以是，假想区域设定单元在与预先设定的距离相应的范围的构造物的表面的附近设定多个假想区域。根据该结构，能够仅设定必要量的假想区域，能够更高效地进行粒子与构造物的模拟。

可以是，时间迁移单元基于由相互作用力计算单元计算出的相互作用力，计算下一时间步骤的构造物的位置和速度。根据该结构，也能够模拟粒子对构造物的影响，因此能够进行更加准确的粒子与构造物的模拟。

可以是，假想区域为粒子状的区域。根据该结构，能够适当且容易地处理假想区域，能够适当且容易地进行粒子与构造物的模拟。

然而，本发明除了如上所述能够记述为粒子模拟装置的发明之外，也能如下所示记述为粒子模拟方法和粒子模拟程序的发明。这只是范畴不同，实质上是相同的发明，起到相同的作用和效果。

即，本发明的一个实施方式的粒子模拟方法是基于作业空间内的粒子与构造物之间的相互作用力计算所述粒子的位置和速度并模拟所述粒子和所述构造物的动态的粒子模拟装置的动作方法，其中，包括：构造物初始信息取得步骤，取得表示构造物的初始位置和形状的构造物初始信息；假想区域设定步骤，基于由在构造物初始信息取得步骤中取得的构造物初始信息表示的构造物的初始位置和形状，来在所述构造物的表面的附近设定多个假想区域，关于所设定的多个假想区域的各个假想区域，计算到构造物的距离；位置信息取得步骤，取得表示粒子的位置的位置信息；假想区域确定步骤，根据由在位置信息取得步骤中取得的位置信息表示的粒子的位置，来确定预先设定的距离内的假想区域；粒子距离计算步骤，关于在假想区域确定步骤中确定出的假想区域，基于在假想区域设定步骤中计算出的距离，来计算从粒子到构造物的距离；相互作用力计算步骤，基于在粒子距离计算单元中计算出的距离，计算粒子与构造物之间的相互作用力；以及时间迁移步骤，基于在相互作用力计算步骤中计算出的相互作用力，计算下一时间步骤的粒子的位置和速度，并且在维持与该构造物之间的位置关系的状态下，使假想区域移动。

另外，本发明的一个实施方式的粒子模拟程序是使计算机作为粒子模拟装置发挥功能的粒子模拟程序，所述粒子模拟装置基于作业空间内的粒子与构造物之间的相互作用力，计算所述粒子的位置和速度，模拟所述粒子和所述构造物的动态，其中，所述粒子模拟程序使计算机作为构造物初始信息取得单元、假想区域设定单元、位置信息取得单元、假想区域确定单元、粒子距离计算单元、相互作用力计算单元以及时间迁移单元发挥功能，所述构造物初始信息取得单元取得表示构造物的初始位置和形状的构造物初始信息，所述假想区域设定单元基于由所述构造物初始信息取得单元取得的构造物初始信息所表示的构造物的初始位置和形状，来在所述构造物的表面的附近设定多个假想区域，关于所设定的多个假想区域的各个假想区域，计算到构造物的距离，所述位置信息取得单元取得表示粒子的位置的位置信息，所述假想区域确定单元根据由所述位置信息取得单元取得的位置信息所表示的粒子的位置，来确定预先设定的距离内的假想区域，所述粒子距离计算单元关于由所述假想区域确定单元确定出的假想区域，基于由所述假想区域设定单元计算出的距离，来计算从粒子到构造物的距离，所述相互作用力计算单元基于由所述粒子距离计算单元计算出的距离，计算粒子与构造物之间的相互作用力，所述时间迁移单元基于由所述相互作用力计算单元计算出的相互作用力，计算下一时间步骤的粒子的位置和速度，并在维持与所述构造物之间的位置关系的状态下，使所述假想区域移动。

发明效果

在本发明的一实施方式中，无需重新计算假想区域的距离(距离函数)。另外，与现有的距离函数不同，仅对在构造物的附近设定的假想区域计算距离即可。另外，假想区域针对每个构造物设定，因此在存在多个构造物的情况下、即存在多个互不相同的移动轴、旋转轴的情况下，也能够使用。另外，由于与距离连接了关联的假想区域移动，因此也不会产生因移流计算所致的数值扩散。即，根据本发明的一个实施方式，能够以较少的计算负荷准确地进行粒子与构造物之间的距离计算，由此，能够适当地实施粒子与构造物的模拟。

附图说明

图1是示出本发明的实施方式的粒子模拟装置的功能结构的图。

图2是示意性地示出与构造物相对的假想粒子的图。

图3是示出供与构造物相对的假想粒子配置的范围和配置了与构造物相对的假想粒子的状态的图。

图4是示意性地示出距模拟对象的粒子在所设定的距离内的假想粒子的图。

图5是示出在本发明的实施方式的粒子模拟装置中执行的处理(粒子模拟方法)的流程图。

图6是示出本发明的实施方式中的模拟的执行中的状态(快照)的图。

图7是将本发明的实施方式的粒子模拟程序的结构与记录介质一起示出的图。

具体实施方式

下面，结合附图详细说明本发明的粒子模拟装置、粒子模拟方法以及粒子模拟程序的实施方式。此外，在附图说明中，对相同要素标注相同的附图标记，并省略重复说明。

图1示出本实施方式的粒子模拟装置10。粒子模拟装置10是模拟(解析)作业空间(计算区域)内的多个球形(从中心到表面的距离一定的形状)的粒子和构造物的动态的装置。具体而言，粒子模拟装置在每个模拟上的时刻即时间步骤，基于各粒子和构造物的位置和速度，计算对各粒子作用的力。对各粒子和构造物作用的力包括例如，因粒子与构造物之间以及粒子间的相互作用即接触(碰撞)所产生的相互作用力亦即接触力。此外，粒子与构造物之间以及粒子间的相互作用并非一定限于接触，只要影响位置和速度，可以设定为任意的相互作用。粒子模拟装置10基于所计算出的力，计算下一时间步骤中的各粒子和构造物的位置和速度。此外，本实施方式中的粒子模拟基于DEM实施，但也可以基于SPH、MPS等其他粒子法实施。另外，根据粒子法的种类不同，也可以不将粒径、旋转速度作为信息保存。

成为本实施方式的粒子模拟装置10的模拟的对象的粒子包括成为现有粒子模拟的对象的任意的粒子。例如，能够以沙土、粉体为对象。或者，也可以将流体、固体假定为由多个粒子构成的物质，而作为对象。成为本实施方式的粒子模拟装置10的模拟的对象的构造物包括任意的物体。例如，以搅拌、混合粉体的装置(的搅拌、混合粉体的部分(搅拌叶片))为对象。成为模拟的对象的构造物也可以是多个。

通过本实施方式的粒子模拟装置10的模拟，能够模拟物理性问题。例如，能够实施有关粉体的搅拌、混合的模拟。或者，也能考虑用于与钻头的形状相应的钻削能力的验证、推土机等建设机械、重型机械的设计。这样，本实施方式的粒子模拟装置10的模拟也能考虑用于土木、粉体等工业领域。另外，也能够用于对堤防、建筑物、地形与海啸、沙土的相互作用、泥石流、雪崩等自然现象的模拟。这样，本实施方式的粒子模拟装置10的模拟也能考虑用于防灾领域。

粒子模拟装置10例如构成为具备CPU(Central Processing Unit：中央处理器)、GPU(Graphics Processing Unit：图形处理器)、存储器、硬盘、显示器等硬件的计算机。这些构成要素通过程序等而动作，由此发挥作为后述的粒子模拟装置10的功能。粒子模拟装置10在能够进行并行运算的装置中特别有效地进行动作。此外，粒子模拟装置10未必需要具备CPU作为运算装置，也可以是仅具备CPU的结构。

如图1所示，粒子模拟装置10具备模拟信息保存部11、构造物初始信息取得部12、假想区域设定部13、位置信息取得部14、假想区域确定部15、粒子距离计算部16、相互作用力计算部17以及时间迁移部18而构成。

在本实施方式中，粒子移动的区域即作业区域例如为三维的空间，被分割(区分)为预先设定了一边的大小的立方体的单元(网格)。此外，单元也可以是立方体以外的形状、例如长方体。粒子模拟装置10在实施模拟的处理之前将作业空间预先分割为单元，预先掌握作业空间是如何被分割为单元的。上述一边的大小例如根据成为模拟的对象的多个粒子的粒径等预先设定。另外，对作业区域内的各单元标注有确定单元的单元号。单元号例如与作业空间内的单元的位置相对应地依次标注。

模拟信息保存部11是保存模拟所使用的信息的单元。具体而言，模拟信息保存部11保存作业区域内的多个粒子各自的粒子信息。粒子信息包括表示粒子的坐标、粒子的速度以及粒子半径的信息。粒子的坐标是表示作业空间中的粒子的位置的三维坐标。粒子的速度包括平移速度和旋转速度。关于粒子的坐标和粒子的速度，模拟的开始时的信息(初始信息)由粒子模拟装置10的用户等预先输入模拟信息保存部11，另外，正在模拟的信息通过后述的时间迁移部18更新。粒子半径由粒子模拟装置10的用户等预先输入模拟信息保存部11。

模拟信息保存部11保存作业区域内的各构造物的构造物信息。构造物信息包括表示构造物的坐标、构造物的速度以及构造物的形状的信息。构造物的坐标是表示作业空间中的构造物的位置(例如，构造物的特定部位的位置)的三维坐标。构造物的速度包括平移速度和旋转速度。构造物的坐标和构造物的速度与粒子的坐标和粒子的速度相同，模拟的开始时的信息(初始信息)由粒子模拟装置10的用户等预先输入模拟信息保存部11，另外，正在模拟的信息通过后述的时间迁移部18更新。表示构造物的形状的信息例如是以多边形表示构造物的形状的CAD(computer-aided design：计算机辅助设计)数据。表示构造物的形状的信息由粒子模拟装置10的用户等预先输入模拟信息保存部11。构造物信息也可以包括表示构造物如何移动(例如，移动轴、旋转轴等)的信息。该信息也由粒子模拟装置10的用户等预先输入模拟信息保存部11。

另外，模拟信息保存部11如后所述，保存假想粒子的信息。另外，模拟信息保存部11也可以预先输入并保存上述粒子信息和构造物信息以外的用于模拟的信息。作为这样的信息，是摩擦系数、弹性系数、粘性阻尼系数、恢复系数等。

构造物初始信息取得部12是取得表示构造物的初始位置和形状的构造物初始信息的构造物初始信息取得单元。例如，构造物初始信息取得部12取得保存于模拟信息保存部11中的构造物的初始信息和表示构造物的形状的信息作为构造物初始信息。构造物初始信息取得部12将所取得的构造物初始信息向假想区域设定部13输出。

假想区域设定部13是基于通过由构造物初始信息取得部12取得的构造物初始信息表示的构造物的初始位置和形状，在该构造物的表面的附近设定多个假想区域，并关于所设定的多个假想区域分别计算到构造物的距离(最接近距离)的假想区域设定单元。如后所述，假想区域是用于计算从模拟对象的粒子到构造物的距离的区域。在本实施方式中，假想区域是球形(从中心到表面的距离一定的形状)的粒子状的区域，叫做假想粒子。多个假想粒子的直径的长度相等，并预先设定。假想粒子可以考虑计算的负荷等，设定为比模拟对象的粒子大。

如上所述，以往以来，从模拟对象的粒子到构造物的距离如非专利文献1所示，使用符号距离函数来求出。例如，符号距离函数是在作业区域被分割为单元(正交格子)而产生的格子点处存储粒子到构造物的距离。图2(a)示出用线连接基于格子点状的距离函数值生成的到构造物100等距离的位置所得到的等值线图。在本实施方式中，如图2(b)所示，在构造物100的表面(构造物100与除此以外的部分的分界)的附近设定多个假想粒子110。

具体而言，假想区域设定部13首先计算由构造物初始信息表示的初始位置和形状的构造物100的符号距离函数。该计算只要使用如非专利文献2中记载的那样的现有方法实施即可。假想区域设定部13基于计算出的符号距离函数，检测作业区域内的构造物100的附近的位置。如图2(b)所示，构造物的附近包括构造物的外侧、表面以及内侧的区域。假想区域设定部13在该附近的区域，以相互接触的方式设定多个假想粒子110。例如，在作业区域为三维的情况下，假想区域设定部13在该区域内设定单纯立方格子或者面心立方格子的格子点，在该格子点设定(配置)假想粒子110。此外，设定假想粒子只要知晓与构造物100的距离即可，未必需要使用符号距离函数，可以以任意方法实施。

设定假想粒子110的构造物100的表面的部分是模拟对象的粒子有可能接触(相互作用)的部分。如图2(b)所示，在粒子有可能遍布构造物100的整个表面接触的情况下，遍及构造物100的整个表面，在表面的附近设定假想粒子110。在预先知晓不可能接触的位置的情况下，可以不在该部分的附近设定假想粒子110。

假想区域设定部13在构造物100的表面的法线方向上，将规定的距离(规定的厚度)的范围设定为构造物的附近。例如，如图3所示，将从构造物100的表面向外侧离开了距离d的位置作为该距离的中心，将从此处向外侧和内侧(图3中的上下方向)分别离开了距离H的区域作为构造物100的附近的区域。距离d是构造物100与模拟对象的粒子101开始相互作用的构造物100的表面与粒子101的中心之间的距离。此外，模拟对象的粒子101的大小并不均一等，在开始相互作用的距离按照每个粒子101不均一的情况下，例如，关于图3中的d-H可以使用它们的最小值，关于d+H可以使用它们的最大值。

在相互作用为接触的情况下，距离d是粒子101的半径。距离H是用于确定为了计算从粒子101到构造物100的距离而使用的假想粒子110的距离(详情后述)。距离H预先设定。如上所述，假想区域设定部13也可以在与距离H相应的范围的构造物100的表面的附近设定多个假想粒子110。此外，在距离d比距离H大的情况下，存在在构造物100的表面和内侧不设定假想粒子110的情况。另外，供设定假想粒子110的区域(法线方向上的距离、厚度)无需如上所述进行设定(计算)，也可以预先设定。

如后所述，构造物100与所设定的假想粒子110在模拟期间，位置关系被固定。

假想区域设定部13关于所设定的各假想粒子110，计算从该假想粒子110到(所配置的表面的)构造物100为止的距离，例如，从假想粒子110的中心的位置到构造物100的表面的距离。该计算能够使用基于运用了以往的符号距离函数的内插插值的方法等任意的方法。此外，从假想粒子110到构造物100的距离的计算也可以结合假想粒子110的配置来进行。所设定的假想粒子110预先设为对每个假想粒子110标注唯一性的编号等，以便能够唯一地确定。假想区域设定部13将表示所设定的假想粒子110的位置和上述距离的信息储存于模拟信息保存部11中进行保存。

假想区域设定部13可以在计算假想粒子110的距离时，预先计算与构造物100的表面相对的法线矢量，并与距离对应地建立起关联。法线矢量的信息能够在后述的内插插值的处理中使用。另外，在假想粒子110中，还能够预先将杨氏模量等机械物性建立关联，而使用于后述的处理。

上述构造物初始信息取得部12和假想区域设定部13的处理在模拟的各时间步骤的运算之前实施。以下的功能部是用于进行模拟的各时间步骤的运算的结构。

位置信息取得部14是取得表示成为模拟对象的粒子101的位置的位置信息的位置信息取得单元。具体而言，位置信息取得部14取得保存于模拟信息保存部11中的各粒子101的当前时间步骤的粒子信息作为粒子101的位置信息。位置信息取得部14将所取得的粒子101的粒子信息输入假想区域确定部15。

位置信息取得部14也取得表示假想粒子110的位置的位置信息。位置信息取得部14取得表示保存于模拟信息保存部11中的各粒子101的当前时间步骤的位置的信息。位置信息取得部14将所取得的假想粒子110的位置信息向假想区域确定部15输出。

假想区域确定部15是根据通过由位置信息取得部14取得的位置信息表示的粒子101的位置来确定预先设定的距离H内的假想粒子110的假想区域确定单元。这里所确定出的假想粒子110用于利用内插插值计算从粒子101到构造物100的距离。即，距离H是表示进行内插插值的范围(插值范围)的距离，预先设定能够利用内插插值适当地计算从粒子101到构造物100的距离那样的值。具体而言，距离H优选为，在假想粒子110设定为单纯立方格子状的情况下，满足H＞D√n，在将假想粒子110设定为面心立方格子状的情况下，满足H＞D。此外，在上述关系式中，如图3所示，D是假想粒子110的直径，n是作业区域的维数。

如图4所示，假想区域确定部15针对每个粒子101，计算由从位置信息取得部14输入的信息表示的粒子101(粒子i)的位置(粒子101的中心的位置)与各假想粒子110(粒子j)的位置(假想粒子110的中心的位置)之间的距离Pij。假想区域确定部15判断所计算出的距离Pij是否在距离H以内。假想区域确定部15基于该判断，针对每个粒子101，确定与粒子101的距离在距离H以内的假想粒子110。在图4所示的例子中，确定4个假想粒子110a。假想区域确定部15将表示针对每个粒子101确定出的假想粒子110的信息向粒子距离计算部16输出。

此外，距离H内的假想粒子110的确定可以利用以往的粒子模拟中的粒子间的接触判定(判定粒子间距离是否在阈值以内)的方法，高效地进行。例如，可以采用，日本特开2015-115567号公报(专利文献1)和日本特开2010-238030号公报(专利文献2)中记载的那样的，基于粒子101和假想粒子110所在的单元来设定粒子101与假想粒子110的对，并关于该对进行距离的判断的方法。

粒子距离计算部16是关于由假想区域确定部15确定的假想粒子110，基于由假想区域设定部13计算出的距离来计算从粒子101到构造物100的距离(最接近距离)的粒子距离计算单元。粒子距离计算部16取得保存于模拟信息保存部11中的、关于针对每个粒子101确定出的假想粒子110的表示到构造物100的距离的信息。粒子距离计算部16基于由所取得的信息表示的距离，计算从粒子101到构造物100的距离(计算图4所示的从粒子101(粒子i)的中心的位置到构造物100的表面的距离disi。该计算例如可以通过非专利文献1中示出的以往的内插插值来进行，也可以使用其他任意方法实施。例如，也能使用在SPH等中使用的核估计。

在没有确定出的假想粒子110或者确定出的假想粒子110的数量不足以实施内插插值(例如，为预先设定的阈值以下)的情况下，关于该粒子101，设为到构造物100的距离未达到至少产生相互作用的距离，不计算到构造物100的距离。在该情况下，也不进行使用了该距离的处理。粒子距离计算部16将针对每个粒子101计算出的表示到构造物100的距离的信息向相互作用力计算部17输出。

另外，粒子距离计算部16通常在计算上述距离时，计算从粒子101相对于构造物100的表面的法线矢量。通过法线矢量，能够利用矢量来表示由相互作用力计算部17计算的相互作用力。该矢量所表示的方向表示相互作用力所作用的方向。法线矢量能够通过对关于对每个粒子101确定出的假想粒子110的到构造物100的距离(距离函数数据)计算空间梯度来求出。此外，根据空间梯度的法线矢量的计算例如可以采用非专利文献1中示出的以往的方法来进行。或者，在对于假想粒子110预先计算了法线矢量的情况下，可以对该信息进行内插插值来计算。也可以采用其他任意方法实施。例如，也能使用在SPH等中使用的核估计的梯度计算法。粒子距离计算部16将计算出的表示法线矢量的信息与表示距离的信息一并向相互作用力计算部17输出。

相互作用力计算部17是基于由粒子距离计算部16计算出的距离，计算粒子101与构造物100之间的相互作用力的相互作用力计算单元。相互作用力计算部17通过判断从粒子101到构造物100的距离是否小于预先设定的阈值，来进行在粒子101与构造物100之间是否会产生相互作用的判定(相互作用判定)。相互作用力计算部17在判断为该距离小于阈值的情况下，判定为在粒子101与构造物100之间产生相互作用。相互作用力计算部17在判断为该距离不小于阈值的情况下，判定为在粒子101与构造物100之间不产生相互作用。

例如，在相互作用为接触的情况下，相互作用力计算部17判断该距离是否小于作为阈值的粒子101的半径d，在判断为该距离小于阈值的情况下，判定为粒子101与构造物100接触，在它们之间产生接触力。相互作用力的计算能够使用与以往的粒子模拟相同的方法、其他任意方法来进行。相互作用力计算部17通常基于由粒子距离计算部16计算出的从粒子101相对于构造物100的表面的法线矢量，计算相互作用力而作为矢量。

相互作用力计算部17可以不仅计算粒子101与构造物100之间的相互作用力，还计算粒子101彼此的相互作用力以及构造物100彼此的相互作用力。粒子101彼此的相互作用力的计算能够使用与在专利文献1、2中记载的以往的粒子模拟相同的方法、其他任意方法来进行。构造物100彼此的相互作用力的计算也能够使用与以往的模拟相同的方法、其他任意的方法来进行。此外，在相互作用力的计算中，也可以取得保存于模拟信息保存部11中的相互作用力的计算所需的参数，并用于相互作用力的计算。

相互作用力计算部17根据如上所述计算出的各个粒子101与构造物100之间的相互作用力、粒子101彼此的相互作用力以及构造物100彼此的相互作用力，计算关于各粒子101和各构造物100的相互作用力的总和。相互作用力计算部17将表示计算出的各粒子101和各构造物100的相互作用力的总和的信息向时间迁移部18输出。

时间迁移部18是基于由相互作用力计算部17计算出的相互作用力，计算下一时间步骤的粒子101和构造物100的位置和速度并且在维持与该构造物100之间的位置关系的状态下使假想粒子110移动的时间迁移单元。具体而言，时间迁移部18针对每个粒子101和构造物100取得保存于模拟信息保存部11中的粒子信息和构造物信息，根据粒子信息和构造物信息中所示的当前时间步骤的坐标和速度以及相互作用力的总和，计算下一时间步骤的坐标和速度。该计算例如能够使用与以往的模拟相同的方法、其他任意的方法来进行。时间迁移部18以计算出的下一时间步骤中的粒子101和构造物100的位置和速度，来针对每个粒子101和构造物100更新保存于模拟信息保存部11中的粒子信息和构造物信息。

时间迁移部18在关于构造物100计算下一时间步骤的位置后，在维持与该构造物100之间的位置关系的状态下(与该构造物100的移动同样地)，使针对该构造物100设定的假想粒子110移动。时间迁移部18以表示移动后的假想粒子110的位置的信息，更新保存于模拟信息保存部11中的表示假想粒子110的信息(将该信息作为表示下一时间步骤的假想粒子110的位置的信息，保存于模拟信息保存部11)。

此外，关于构造物100，也可以不进行基于相互作用力的位置和速度的计算。例如，在粒子101对构造物100的影响极小并且不产生构造物100彼此的相互作用的情况下，也可以不进行该计算。在该情况下，例如，可以设为，表示构造物100预先将如何移动的信息存储于模拟信息保存部11，时间迁移部18据此计算构造物100的移动(下一时间步骤的位置和速度)。另外，时间迁移部18根据构造物100的该移动(根据预先决定的构造物100的表面的轨迹)，使假想粒子110移动。

在通过时间迁移部18进行关于所有粒子101、构造物100以及假想粒子110的更新后，实施下一时间步骤的处理。

在粒子模拟装置10中，每当一个时间步骤中的运算结束时，判断是否满足模拟的结束条件。例如，在预先设定的次数(时间步骤)的运算结束了的情况下，判断为满足了结束条件。在判断为满足了结束条件的情况下，在粒子模拟装置10中，模拟结束。在该情况下，例如，进行将运算结果向显示装置、其他装置的输出等。在判断为没有满足结束条件的情况下，反复实施下一时间步骤的运算。以上为粒子模拟装置10的结构。

接着，使用图5的流程图，说明本实施方式的粒子模拟装置10的动作方法即由粒子模拟装置10执行的处理(粒子模拟方法)。本处理例如通过粒子模拟装置10的用户对粒子模拟装置10进行开始模拟的操作而开始。

在粒子模拟装置10中，通过构造物初始信息取得部12取得保存于模拟信息保存部11中的表示构造物100的初始位置和形状的构造物初始信息(S01，构造物初始信息取得步骤)。图6(a)示出构造物100的形状的例子。所取得的构造物初始信息从构造物初始信息取得部12向假想区域设定部13输出。

接着，通过假想区域设定部13计算构造物初始信息的构造物100的符号距离函数(S02，假想区域设定步骤)。接着，通过假想区域设定部13，使用该符号距离函数来在构造物100的表面的附近配置(设定)多个假想粒子110(S03，假想区域设定步骤)。图6(b)示出在构造物100的表面的附近配置的假想粒子110的例子。接着，通过假想区域设定部13，关于假想粒子110计算从该假想粒子110到构造物100的距离(S04，假想区域设定步骤)。表示假想粒子110的位置和距离的信息保存于模拟信息保存部11中。以上是在模拟的各时间步骤的运算之前进行的处理。

接着，通过位置信息取得部14取得保存于模拟信息保存部11中的表示成为模拟对象的粒子101和假想粒子110的位置的位置信息(S05，位置信息取得步骤)。所取得的位置信息从位置信息取得部14向假想区域确定部15输出。

接着，通过假想区域确定部15，来确定由位置信息表示的距离粒子101的位置在距离H内的假想粒子110(S06，假想区域确定步骤)。针对每个粒子101确定出的表示假想粒子110的信息从假想区域确定部15向粒子距离计算部16输出。

接着，通过粒子距离计算部16，基于关于针对每个粒子101确定出的假想粒子110的到构造物100的距离，来计算从粒子101到构造物100的距离(S07，粒子距离计算步骤)。该计算如上所述例如通过内插插值来进行。针对每个粒子101计算出的表示到构造物100的距离的信息从粒子距离计算部16向相互作用力计算部17输出。

接着，通过相互作用力计算部17，基于针对每个粒子101计算出的到构造物100的距离，计算粒子101与构造物100之间的相互作用力(S08，相互作用力计算步骤)。此时，也可以计算粒子101彼此的相互作用力和构造物100彼此的相互作用力。接着，通过相互作用力计算部17，根据计算出的相互作用力，计算关于各粒子101和各构造物100的相互作用力的总和。表示计算出的各粒子101和各构造物100的相互作用力的总和的信息从相互作用力计算部17向时间迁移部18输出。

接着，通过时间迁移部18，基于计算出的相互作用力计算下一时间步骤的粒子101和构造物100的位置和速度并且在维持与该构造物100之间的位置关系的状态下，使假想粒子110移动(S09，时间迁移步骤)。基于上述运算，将保存于模拟信息保存部11中的粒子信息、构造物信息、表示假想粒子110的位置的信息更新为下一时间步骤的信息。

接着，在粒子模拟装置10中，判断是否满足了模拟的结束条件(S10)。在判断为满足了结束条件的情况下(S10的是)，处理(模拟)结束。在判断为不满足结束条件的情况下(S10的否)，时间步骤前进一步，实施下一时间步骤中的上述处理(S05～S10)。

图6(c)～图6(f)示出各时间步骤中的状态的例子。随着从图6(c)进入图6(f)，时间步骤推进。在该模拟中，在容器的内侧设置有2个用于搅拌粉体的旋转的构造物(搅拌叶片)，粉体从上部掉落到容器中而由构造物搅拌。此外，在各构造物中，能够分别施加其他搅拌速度。以上是在本实施方式的粒子模拟装置10中执行的处理。

如上所述，根据本实施方式，基于关于在构造物100的表面的附近设定的多个假想粒子110而计算出的距离，进行粒子101与构造物100之间的距离计算。在模拟中，该假想粒子110在维持与构造物100之间的位置关系的状态下移动。因此，无需重新计算关于假想粒子110的距离(距离函数)。

另外，与遍布作业空间内全域的格子点处的以往的距离函数不同，只要仅对在构造物100的附近设定的假想粒子110计算距离即可。另外，各个假想粒子110能够针对每个构造物100设定，随着每个构造物100的动作独立地移动。即，在存在多个构造物100的情况下，即，在存在多个互不相同的移动轴、旋转轴的情况下，也能够使用。

另外，由于与距离建立了关联的假想粒子110移动，因此不会产生因移流计算所致的数值扩散，无需距离函数的初始化处理。另外，由于只是移动假想粒子110，因此与移流计算相比，非常简单，计算负荷也较少。即，根据本实施方式，能够以较少的计算负荷准确地实施粒子101与构造物100之间的距离计算，由此能够适当地进行粒子101与构造物100的模拟。

另外，如上所述，由于仅针对在构造物100的附近设定的假想粒子110计算距离并存储即可，因此与以往的距离函数相比，能够压倒性地减少要存储的数据量(存储器使用量)。

以下示出存储器使用量的例子。示出将作业空间的尺寸设为160D×200D×100D，将插值范围H设为3.5d，将假想粒子110的直径(格子尺寸)设为D，将相互作用距离d设为0.5D的情况的例子。在本实施方式的方法中，距离计算所需的粒子数为85，336，总存储器使用量为2，730，752(Byte)(85，336×4×8：每个粒子需要4个(x，y，z坐标和距离函数值)8byte的变量)。在以往的方法(在移动粒子的坐标系的情况下)中，距离计算所需的格子数为3，200，000，总存储器使用量为102，400，000(Byte)(3，200，000×4×8：每个格子需要4个(x，y，z坐标和距离函数值)8byte的变量)。而且，在进行距离函数的移流计算的情况下，例如，若使用CIP法(Constrained Interpolation Profile Scheme：紧致插值曲线方法)，则存储器使用量为102，400，000+3，200，000×7×8＝281，600，000(Byte)(每个格子需要7个8byte的变量)。另外，在使用距离函数的时序数据的情况下，存储器使用量为102，400，000+3，200，000×8＝128，000，000(Byte)(为了实施时间插值，需要2个时刻的距离函数数据，距离函数值为每个格子8Byte，全部需要追加3，200，000×8Byte)。

如上述例子所示，在本实施方式中，能够将存储器使用量减少到以往的约1/40～1/100。

另外，如本实施方式所示，可以将设定(配置)假想粒子110的范围设定为与确定用于计算从粒子101至构造物100的距离的假想粒子110的范围(插值范围)相应。根据该结构，能够仅设定必要量的假想粒子110，能够更高效地实施粒子101与构造物100的模拟。

另外，如本实施方式那样，构造物100可以受到来自粒子101的相互作用力。根据该结构，也能够模拟粒子101对构造物100的影响，因此能够进行更准确的粒子101与构造物100的模拟。

另外，如本实施方式那样，也可以将假想区域设为假想粒子110。根据该结构，能够适当且容易地处理假想区域。例如，如上所述，能够容易地使用以往的粒子基体的模拟的方法，能够使本实施方式的方法高速化。即，能够适当且容易地进行粒子101与构造物100的模拟。其中，未必需要将假想区域设为粒子状的区域，可以设为任意形状的区域。

此外，在上述实施方式中，将作业区域设为了三维的空间，但也可以是三维以外的空间，例如二维的空间。

接着，对用于使计算机执行上述一系列的由粒子模拟装置10实施的处理的粒子模拟程序进行说明。如图7所示，粒子模拟程序30插入计算机而被访问，或者储存于在计算机所具备的记录介质20中形成的程序储存区域21内。

粒子模拟程序30具备模拟信息保存模块31、构造物初始信息取得模块32、假想区域设定模块33、位置信息取得模块34、假想区域特定模块35、粒子距离计算模块36、相互作用力计算模块37以及时间迁移模块38而构成。通过使模拟信息保存模块31、构造物初始信息取得模块32、假想区域设定模块33、位置信息取得模块34、假想区域特定模块35、粒子距离计算模块36、相互作用力计算模块37以及时间迁移模块38执行而实现的功能分别与上述的粒子模拟装置10的模拟信息保存部11、构造物初始信息取得部12、假想区域设定部13、位置信息取得部14、假想区域确定部15、粒子距离计算部16、相互作用力计算部17以及时间迁移部18是同样的。

此外，粒子模拟程序30还可以构成为，其一部分或者全部经由通信线路等传送介质传送，由其他设备接收并记录(包含安装)。另外，粒子模拟程序30的各模块可以不是安装于一台计算机，而是安装于多台计算机中的任一台。在该情况下，通过该多台计算机的计算机系统来进行上述一系列的粒子模拟程序30的处理。

标号说明

10…粒子模拟装置；11…模拟信息保存部；12…构造物初始信息取得部；13…假想区域设定部；14…位置信息取得部；15…假想区域确定部；16…粒子距离计算部；17…互相作用力计算部；18…时间迁移部；20…记录介质；21…程序储存区域；30…粒子模拟程序；31…模拟信息保存模块；32…构造物初始信息取得模块；33…假想区域设定模块；34…位置信息取得模块；35…假想区域确定模块；36…粒子距离计算模块；37…相互作用力计算模块；38…时间迁移模块。