一种离散元仿真中多球颗粒的两层网格搜索接触检测方法

摘要

本发明涉及一种离散元仿真中多球颗粒的两层网格搜索接触检测方法，包括以下步骤：全局空间被划分为尺寸均为D的方形或立方体网格，然后进行第一层网格搜索，如果目标颗粒的包络球与某个候选颗粒的包络球相交，则第一层网格搜索完成，反之则进行目标颗粒的包络球与其余候选颗粒的包络球之间的相交检测；相交时，一个尺寸为(D+d)的方形或立方体局部空间被划分为尺寸均为d的方形或立方体网格，并供第二层网格搜索所用，检测目标颗粒的单元球与候选颗粒的单元球之间的接触。本发明在实现接触检测的内存消耗大幅降低的同时，维持接触检测的时间消耗不增加，为多球颗粒离散元仿真在工业领域的大规模应用提供了有效的技术手段。

说明书

技术领域

本发明涉及基于离散元法的非球形颗粒介质动力学仿真领域，尤其是涉及一种离 散元仿真中多球颗粒的两层网格搜索接触检测方法。

背景技术

由于几何形状简单，球形颗粒被广泛地应用于三维离散元仿真之中。然而，真实的 颗粒形状难以由简单的球体来有效地表达。因此，非球形颗粒的形状表达已在相关专利和 文献资料中进行了广泛的讨论。非球形颗粒的形状表达方法基本上可以分为两类：单颗粒 方法和多球颗粒方法。

在单颗粒方法中，一个非球形颗粒是由一个椭圆、椭球或超二次曲面等复杂几何 形状的单颗粒来近似表达。二维椭圆颗粒和三维椭球颗粒的几何形状可由数学方程来描 述，二维椭圆颗粒间或三维椭球颗粒间的接触检测通常需要通过迭代求解多项式方程的根 来实现。超二次曲面颗粒的几何形状也可由数学方程来描述，在三维方面，超二次曲面颗粒 是椭球颗粒的一种延伸，而在二维方面，超二次曲线颗粒是椭圆颗粒的一种扩展；超二次曲 面颗粒为三维或二维颗粒凸凹面的形状表达提供了丰富的手段；与椭圆和椭球颗粒相类 似，超二次曲面颗粒间的接触检测需要通过迭代求解非线性方程组的根来实现。椭圆、椭球 或超二次曲面等复杂几何形状的颗粒间的接触检测在大规模的离散元仿真中会消耗大量 的内存和计算时间。

在多球颗粒方法中，一个非球形颗粒可以通过若干个形状简单的单元球组合而成 的多球颗粒来近似表达，这些单元球可以相互重叠，并且它们的尺寸可以被改变。一般而 言，一个多球颗粒由若干个具有较小尺寸的单元球组成。在一些多球颗粒的离散元仿真中， 一个复杂的颗粒形状甚至需要由成百上千个单元球来表达。在这种情况下，多球颗粒尺寸 是单元球尺寸的若干倍。近年来，多球颗粒已经被用于描述各类颗粒形状，包括椭圆形颗 粒、胶囊形颗粒、玉米状颗粒、谷物状颗粒、药片状颗粒等。由于发生在多球颗粒间的接触实 际上是单元球间的接触，因此多球颗粒间的接触检测可以使用球形颗粒间的接触检测方法 来处理。

接触检测的内存消耗和时间消耗是衡量一种接触检测方法有效性的两个重要指 标。在接触检测中，网格搜索方法是一种用于球形颗粒邻居搜索的有效方法，它被用于搜索 一个目标颗粒周围一定区域内的候选颗粒。在网格搜索方法中，全局空间被划分为尺寸均 为d的方形或立方体网格，其中d是最大球形颗粒的直径。相比于搜索整个全局空间，在这种 方法中，每个目标颗粒只需要搜索位于它周围的9个(二维)或27个(三维)邻居网格中的邻 居颗粒，因此球形颗粒间的接触检测的时间消耗为O[M²/(L/d)²](二维)或O[M²/(L/d)³](三 维)，其中M是球形颗粒的总数量，L是全局空间的尺寸。

在涉及多球颗粒的离散元仿真中，上述网格搜索方法经常被直接用于多球颗粒的 接触检测之中。全局空间可以被划分为尺寸均为d的网格，其中d是最大单元球的直径；多球 颗粒中的单元球是该方法在全局空间中的搜索对象；如果候选单元球位于一个目标单元球 周围的9个(二维)或27个(三维)邻居网格中，这个目标单元球就需要与那些候选单元球进 行接触检测。在该方法中，接触检测的内存消耗为O[(L/d)²+M·n_s](二维)或O[(L/d)³+M· n_s](三维)；接触检测的时间消耗为(二维)或(三维)， 其中M是多球颗粒的总数量，n_s是单个多球颗粒中单元球的数量。

接触检测是多球颗粒离散元仿真的一个关键步骤。当网格搜索被直接实施在每个 多球颗粒的每个单元球上时，多球颗粒间的接触检测的内存消耗和时间消耗都远大于同等 数量理想球形颗粒间的接触检测的内存消耗和时间消耗。然而，目前已发表的专利或者文 献资料中很少详细地讨论关于提高多球颗粒的接触检测效率的方法。

发明内容

本发明主要是解决现有技术所存在的问题，提供了一种通过对全局空间进行网格 尺寸为D的网格划分，并对多球颗粒的包络球进行第一层网格搜索，然后对局部空间进行网 格尺寸为d的网格划分，并对包络球相交的每对多球颗粒的单元球进行第二层网格搜索，避 免了在全局空间中对每个多球颗粒的每个单元球进行接触检测，从而在实现接触检测的内 存消耗大幅降低的同时，维持接触检测的时间消耗不增加，为多球颗粒离散元仿真在工业 领域中的大规模应用提供了有效的技术手段。

本发明的上述技术问题主要是通过下述技术方案得以解决的：

一种离散元仿真中多球颗粒的两层网格搜索接触检测方法，其特征在于该方法包 括以下步骤：

步骤1，每个多球颗粒各自由一个包络球来表示；全局空间被划分为尺寸均为D的 方形或立方体网格，并供第一层网格搜索所用，其中D是最大包络球的直径；每个包络球被 投影至全局空间的网格中；

步骤2，在第一层网格搜索中，通过网格搜索来查找目标颗粒的包络球i周围的所 有邻居包络球MM(i)，这些邻居包络球中的候选颗粒被挑选出来；如果目标颗粒的包络球i 与某个候选颗粒的包络球j相交，则该候选颗粒的第一层网格搜索完成，进入步骤3；反之则 进行目标颗粒的包络球i与其余候选颗粒的包络球之间的相交检测；

步骤3，当目标颗粒的包络球i与候选颗粒的包络球j相交时，一个尺寸为(D+d)的 方形或立方体局部空间被划分为尺寸均为d的方形或立方体网格，并供第二层网格搜索所 用，其中d是最大单元球的直径；目标颗粒i与候选颗粒j的单元球被投影至局部空间的网格 中；

步骤4，在第二层网格搜索中，通过网格搜索来查找目标单元球u周围的所有邻居 单元球nn(u)，这些邻居单元球为目标单元球u的候选单元球；如果目标单元球u与某个邻居 单元球v接触，则计算多球颗粒的接触力；反之则进行目标单元球u与其余候选单元球之间 的接触检测；

步骤5，重复步骤4直至遍历目标单元球u周围的所有候选单元球nn(u)和所有目标 单元球n_s；

步骤6，重复步骤2到步骤5直至遍历目标颗粒包络球i周围的所有候选包络球MM (i)和所有目标包络球M。

在上述的一种离散元仿真中多球颗粒的两层网格搜索接触检测方法，所述多球颗 粒由单元球组成，用于表达和模拟各种非球形颗粒的几何形状，其中单元球基于以下定义：

所述单元球是组成多球颗粒的最小单元，单元球可以相互重叠，并且它们的尺寸 可以被改变；组成多球颗粒的单元球的数量、尺寸及其组合形式，决定了所表达的多球颗粒 的几何形状精度。

在上述的一种离散元仿真中多球颗粒的两层网格搜索接触检测方法，所述的步骤 1中，所述每个多球颗粒各自由一个包络球来表示，多球颗粒i的包络球的几何中心P_i的全 局坐标和采用如下公式计算：

$q_i^G=\left(\underset{u=1}{\overset{n_s}{\Sigma }}\right(q_u^i·V_u^i\left)\right)/\left(\underset{u=1}{\overset{n_s}{\Sigma }}{V}_u^i\right)$

其中$q_i^G={\left(\begin{array}{ccc}{x}_i^G& y_i^G& z_i^G\end{array}\right)}^T,q_u^i={\left(\begin{array}{ccc}{x}_u^i& y_u^i& z_u^i\end{array}\right)}^T,$和是多球颗粒i的单元球 u的球心的全局坐标，是单元球u的体积，n_s是多球颗粒i中单元球的数量；多球颗粒i的 包络球的半径R_i采用如下公式计算：

$R_i=\underset{u=1}{\overset{n_s}{max}}\left(\right|P_i{C}_u^i|+r_u^i)$

其中是多球颗粒i的包络球与单元球u之间的中心距，是单元球u的半径；

所述全局空间代表了多球颗粒的运动或静止区域，所有的多球颗粒在仿真中都位 于这个空间之中，全局空间的边界由和来定义；

所述全局空间被划分为尺寸均为D的方形或立方体网格，其中网格在x、y和z方向 上的数量分别为和它们采用如下公式计算：

$n_q^G=Int[(q_{max}^G-q_{min}^G)/D]+\psi$

其中$n_q^G={\left(\begin{array}{ccc}{n}_x^G& n_y^G& n_z^G\end{array}\right)}^T,q_{max}^G={\left(\begin{array}{ccc}{x}_{max}^G& y_{max}^G& z_{max}^G\end{array}\right)}^T,q_{min}^G={\left(\begin{array}{ccc}{x}_{min}^G& y_{\min }^G& z_{min}^G\end{array}\right)}^T,$ψ ＝[111]^T；

所述包络球被投影至全局空间的网格中，多球颗粒i的包络球可以获得x、y和z方 向上的网格号码它们采用如下公式计算：

$I_{iq}^G=Int[(q_i^G-q_{min}^G)/D]+\psi$

其中$I_{iq}^G={\left(\begin{array}{ccc}{I}_{ix}^G& I_{iy}^G& I_{iz}^G\end{array}\right)}^T;$多球颗粒的包络球的网格号码被存储在两个数组Q_G和P_G中，其中Q_G是一个用于存储三维网格的、大小为的三维数组，或者是一个用于存 储二维网格的、大小为的二维数组；P_G是一个用于存储三维或二维包络球编号的、大 小为M的一维数组，其中M是多球颗粒的数量；数组Q_G和P_G的总大小为(L/D)²+M(二维)或(L/ D)³+M(三维)，因此第一层网格搜索的内存消耗S_P为：

在上述的一种离散元仿真中多球颗粒的两层网格搜索接触检测方法，所述的步骤 2中，所述候选颗粒可以通过网格搜索来查找目标颗粒i的包络球周围的所有邻居包络球， 并将这些邻居包络球挑选出来，其中邻居包络球位于目标颗粒i的包络球所在的网格以及 此网格周围的8个(二维)或26个(三维)相邻网格之中；

所述判断目标颗粒i的包络球与候选颗粒j的包络球相交的公式如下：

|P_iP_j|≤(R_i+R_j)

其中|P_iP_j|是多球颗粒i的包络球和多球颗粒j的包络球之间的中心距，R_j是多球 颗粒j的包络球的半径；在全局空间中，每个网格中包络球的平均个数为n_CB＝M/(L/D)²(二 维)或n_CB＝M/(L/D)³(三维)；由于每个包络球需要与O(M/(L/D)²)(二维)或O(M/(L/D)³)(三 维)个邻居包络球进行接触检测，因此包络球间的接触检测的时间消耗为：

在上述的一种离散元仿真中多球颗粒的两层网格搜索接触检测方法，所述的步骤 3中，所述局部空间被用于目标多球颗粒的单元球与候选多球颗粒的相应单元球之间的网 格搜索，局部空间的尺寸为(D+d)，多球颗粒i的局部空间的边界坐标和采用如下公式计算：

$\left(\begin{array}{cc}{q}_{min}^i=q_i^G-\left(\right(D+d)/2)\psi & q_{max}^i=q_i^G+\left(\right(D+d)/2)\psi \end{array}\right)$

其中$q_{min}^i={\left(\begin{array}{ccc}{x}_{min}^i& y_{min}^i& z_{min}^i\end{array}\right)}^T,q_{max}^i={\left(\begin{array}{ccc}{x}_{max}^i& y_{max}^i& z_{max}^i\end{array}\right)}^T;$

所述多球颗粒i的局部空间被划分为尺寸均为d的网格，其中网格在x、y和z方向上 的数量分别为和它们采用如下公式计算：

$n_q^L=Int[(q_{max}^i-q_{min}^i)/d]+\psi$

其中$n_q^L={\left(\begin{array}{ccc}{n}_x^L& n_y^L& n_z^L\end{array}\right)}^T;$

所述目标多球颗粒与候选多球颗粒的单元球被投影至局部空间的网格中，判断目 标多球颗粒i的单元球u位于目标多球颗粒i与候选多球颗粒j两者的局部空间的重叠区域 公式如下：

$(x_{min}^j\le x_u^i\le x_{max}^j)\cap (y_{min}^j\le y_u^i\le y_{max}^j)\cap (z_{min}^j\le z_u^i\le z_{max}^j)$

其中和为多球颗粒j的局部空间的边界坐标； 如果多球颗粒i的单元球u位于上述重叠区域内，单元球u被投影至多球颗粒i的局部空间的 网格中，单元球u在多球颗粒i的局部空间中x、y和z方向上的网格号码分别为它们采用如下公式计算：

$I_{Lq}^{i,u}=Int[(q_u^i-q_{min}^i)/d]+\psi$

其中$I_{Lq}^{i,u}={\left(\begin{array}{ccc}{I}_{Lx}^{i,u}& I_{Ly}^{i,u}& I_{Lz}^{i,u}\end{array}\right)}^T;$

判断候选多球颗粒j的单元球v位于上述重叠区域的公式如下：

$(x_{min}^i\le x_v^j\le x_{max}^i)\cap (y_{min}^i\le y_v^j\le y_{max}^i)\cap (z_{min}^i\le z_v^j\le z_{\max }^i)$

其中和是多球颗粒j的单元球v球心的全局坐标；如果多球颗粒j的单元 球v位于上述重叠区域内，单元球v被投影至多球颗粒i的局部空间的网格中，单元球v在多 球颗粒i的局部空间中x、y和z方向上的网格号码分别为它们采用如下公 式计算：

$I_{Lq}^{j,v}=Int[(q_v^j-q_{min}^i)/d]+\psi$

其中$I_{Lq}^{j,v}={\left(\begin{array}{ccc}{I}_{Lx}^{j,v}& I_{Ly}^{j,v}& I_{Lz}^{j,v}\end{array}\right)}^T,q_v^j={\left(\begin{array}{ccc}{x}_v^j& y_v^j& z_v^j\end{array}\right)}^T;$

单元球的网格号码被存储在两个数组Q_L和P_L中，其中Q_L是一个用于存储三维网格 的、大小为的三维数组，或者是一个用于存储二维网格的、大小为的二维 数组，P_L是一个用于存储三维或二维单元球编号的、大小为n_s的一维数组，这两个数组可以 在两个多球颗粒的单元球间的网格搜索中被重复使用；数组Q_L和P_L的总大小为((D+d)/d)²+ n_s(二维)或((D+d)/d)³+n_s(三维)，因此第二层网格搜索的内存消耗S_L为：

在上述的一种离散元仿真中多球颗粒的两层网格搜索接触检测方法，所述的步骤 4中，所述候选单元球v可以通过网格搜索来查找目标单元球u周围的邻居单元球，并将这些 邻居单元球挑选出来，其中邻居单元球位于目标单元球u所在的网格以及此网格周围的8个 (二维)或26个(三维)相邻网格之中，然后进行目标单元球u与候选单元球v之间的接触检 测，判断单元球u与单元球v接触的公式如下：

$\left|C_u^i{C}_v^j\right|\le (r_u^i+r_v^j)$

其中为单元球u与单元球v的球心距，和分别为单元球u和单元球v的半 径；

在二维情况下，多球颗粒i的单元球在每个网格中的平均数量在多球颗粒i与多球颗粒j两者的局部空间的重叠区域中，多球颗粒i的单元球的数量为 其中和L_η分别是两个局部空间的重叠区域在和η方向上的尺寸，它们 可以表达为：

其中|P_iP_j|在d到D之间变化；多球颗粒i中的每个单元球需要与多球颗粒j中的 个邻居单元球进行接触检测，其中因此在二维情况下，多 球颗粒i的单元球与多球颗粒j的单元球之间的接触检测的时间消耗为：

$\left(\begin{array}{c}{T}_E^{i,j}=O[L_{\zeta }*L_{\eta }/{\left(d\right)}^2*n_{\mathrm{CS}}^i*n_{\mathrm{CS}}^j]\\ =O[\left(\frac{(D+d)-\left|P_i{P}_j\right|}{d}\right)*\left(\frac{D+d}{d}\right)\\ =O[\left(\frac{n_s^2}{{((D+d)/d)}^2}\right)*(1-\frac{\left|P_i{P}_j\right|}{(D+d)})]\end{array},*{\left(\frac{n_s}{{((D+d)/d)}^2}\right)}^2]\right)$

以此类推，在三维情况下，多球颗粒i的单元球与多球颗粒j的单元球之间的接触 检测的时间消耗为：

$T_E^{i,j}=O[\left(\frac{n_s^2}{{\left(\right(D+d)/d)}^3}\right)*(1-\frac{\left|P_i{P}_j\right|}{(D+d)})]$

总体来看，在离散元仿真的每一次计算循环中，都要运用如上所述的接触检测步 骤进行多球颗粒与多球颗粒之间的接触检测，最终识别出所有的有效接触信息，供离散元 程序计算多球颗粒当前的接触力。

因此，本发明的有益效果是通过对全局空间进行网格尺寸为D的网格划分，并对多 球颗粒的包络球进行第一层网格搜索，然后对局部空间进行网格尺寸为d的网格划分，并对 包络球相交的每对多球颗粒的单元球进行第二层网格搜索，避免了在全局空间中对每个多 球颗粒的每个单元球进行接触检测，从而在实现接触检测的内存消耗大幅降低的同时，维 持接触检测的时间消耗不增加，为多球颗粒离散元仿真在工业领域的大规模应用提供了有 效的技术手段。

附图说明

图1是本发明的方法流程示意图。

图2是多球颗粒的包络球示意图。

图3是全局空间的二维示意图。

图4是全局空间的二维网格划分示意图。

图5是多球颗粒的包络球的邻居网格搜索示意图。

图6是目标包络球与候选包络球之间的相交检测示意图。

图7是多球颗粒i的局部空间的二维图。

图8是多球颗粒i的局部空间的二维网格划分示意图。

图9是多球颗粒i的局部空间与多球颗粒j的局部空间两者的重叠区域示意图。

图10是多球颗粒的单元球的邻居网格搜索示意图。

图11是目标单元球与候选单元球之间的接触检测示意图。

图12是目标多球颗粒的包络球与候选多球颗粒的包络球的相交情况示意图。

图13a是矩形漏斗与胶囊形多球颗粒的尺寸图。

图13b是矩形漏斗中胶囊形多球颗粒卸落过程的模拟图。

图13c是不同多球颗粒数量下接触检测的内存消耗示意图。

图13d是不同多球颗粒数量下接触检测的时间消耗示意图。

具体实施方式

下面通过实施例，并结合附图，对本发明的技术方案作进一步具体说明，所举实例 只用于解释本发明，并非用于限定本发明的范围。为了表达的清晰，附图以二维形式来显 示，方程以三维形式来表达。

一种离散元仿真中多球颗粒的两层网格搜索接触检测方法，方案实现流程如图1 所示，具体实施例为：

步骤1：首先为每个多球颗粒配置一个相应的包络球，如图2所示，多球颗粒i的包 络球的几何中心P_i的全局坐标和采用如下公式计算：

$q_i^G=\left(\underset{u=1}{\overset{n_s}{\Sigma }}\right(q_u^i·V_u^i\left)\right)/\left(\underset{u=1}{\overset{n_s}{\Sigma }}{V}_u^i\right)$

其中$q_i^G={\left(\begin{array}{ccc}{x}_i^G& y_i^G& z_i^G\end{array}\right)}^T,q_u^i={\left(\begin{array}{ccc}{x}_u^i& y_u^i& z_u^i\end{array}\right)}^T,$和是多球颗粒i的单元 球u的球心的全局坐标，是单元球u的体积，n_s是多球颗粒i中单元球的数量；多球颗粒 i的包络球的半径R_i采用如下公式计算：

$R_i=\underset{u=1}{\overset{n_s}{max}}\left(\right|P_i{C}_u^i|+r_u^i)$

其中是多球颗粒i的包络球与单元球u之间的中心距，是单元球u的半径。

全局空间代表了多球颗粒的运动或静止区域，所有的多球颗粒在仿真中都位于这 个空间之中，如图3所示，全局空间的边界由和来定义。

全局空间被划分为尺寸均为D的方形或立方体网格，其中D是最大包络球的直径， 如图4所示，网格在x、y和z方向上的数量分别为和它们采用如下公式计算：

$n_q^G=Int[(q_{max}^G-q_{min}^G)/D]+\psi$

其中$n_q^G={\left(\begin{array}{ccc}{n}_x^G& n_y^G& n_z^G\end{array}\right)}^T,q_{max}^G={\left(\begin{array}{ccc}{x}_{max}^G& y_{max}^G& z_{max}^G\end{array}\right)}^T,q_{min}^G={\left(\begin{array}{ccc}{x}_{min}^G& y_{\min }^G& z_{min}^G\end{array}\right)}^T,$ψ＝[111]^T。

包络球被投影至全局空间的网格中，如图5所示，多球颗粒i的包络球可以获得x、y 和z方向上的网格号码它们采用如下公式计算：

$I_{iq}^G=Int[(q_i^G-q_{min}^G)/D]+\psi$

其中$I_{iq}^G={\left(\begin{array}{ccc}{I}_{ix}^G& I_{iy}^G& I_{iz}^G\end{array}\right)}^T;$多球颗粒的包络球的网格号码被存储在两个数组Q_G和P_G中，其中Q_G是一个用于存储三维网格的、大小为的三维数组，或者是一个用于存 储二维网格的、大小为的二维数组；P_G是一个用于存储三维或二维包络球编号的、大 小为M的一维数组，其中M是多球颗粒的数量；数组Q_G和P_G的总大小为(L/D)²+M(二维)或(L/ D)³+M(三维)，因此第一层网格搜索的内存消耗S_P为：

步骤2：在第一层网格搜索中，候选颗粒可以通过网格搜索来查找目标颗粒i的包 络球周围的所有邻居包络球，并将这些邻居包络球挑选出来，其中邻居包络球位于目标颗 粒i的包络球所在的网格以及此网格周围的8个(二维)或26个(三维)相邻网格之中。如图5 所示，多球颗粒4、8、9、10、11和14的包络球的几何中心位于目标多球颗粒1的包络球周围的 8个邻居网格之中，这些邻居多球颗粒有一定的概率与目标多球颗粒1发生接触，因此这些 多球颗粒是目标多球颗粒1的候选多球颗粒。

如图6所示，判断目标颗粒i的包络球与候选颗粒j的包络球相交的公式如下：

|P_iP_j|≤(R_i+R_j)

其中|P_iP_j|是多球颗粒i的包络球和多球颗粒j的包络球之间的中心距，R_j是多球 颗粒j的包络球的半径；在全局空间中，每个网格中包络球的平均个数为n_CB＝M/(L/D)²(二 维)或n_CB＝M/(L/D)³(三维)；由于每个包络球需要与O(M/(L/D)²)(二维)或O(M/(L/D)³)(三 维)个邻居包络球进行接触检测，因此包络球间的接触检测的时间消耗为：

如果目标颗粒的包络球i与某个候选颗粒的包络球j相交，则该候选颗粒的第一层 网格搜索完成，进入步骤3；反之则进行目标颗粒的包络球i与其余候选颗粒的包络球之间 的相交检测。

步骤3：目标多球颗粒的局部空间被用于目标多球颗粒的单元球与候选多球颗粒 的相应单元球间的网格搜索，如图7所示，目标多球颗粒i的局部空间的尺寸为(D+d)，多球 颗粒i的局部空间的边界坐标和采用如下公式计算：

$\left(\begin{array}{cc}{q}_{min}^i=q_i^G-\left(\right(D+d)/2)\psi & q_{max}^i=q_i^G+\left(\right(D+d)/2)\psi \end{array}\right)$

其中$q_{min}^i={\left(\begin{array}{ccc}{x}_{min}^i& y_{min}^i& z_{min}^i\end{array}\right)}^T,q_{max}^i={\left(\begin{array}{ccc}{x}_{max}^i& y_{max}^i& z_{max}^i\end{array}\right)}^T.$

多球颗粒i的局部空间被划分为尺寸均为d的网格，如图8所示，其中网格在x、y和z 方向上的数量分别为和它们采用如下公式计算：

$n_q^L=Int[(q_{max}^i-q_{min}^i)/d]+\psi$

其中$n_q^L={\left(\begin{array}{ccc}{n}_x^L& n_y^L& n_z^L\end{array}\right)}^T.$

判断目标多球颗粒i和候选多球颗粒j的单元球是否位于多球颗粒i与多球颗粒j 两者的局部空间的重叠区域之中。如图9所示，判断目标多球颗粒i的单元球u位于上述重叠 区域的公式如下：

$(x_{min}^j\le x_u^i\le x_{max}^j)\cap (y_{min}^j\le y_u^i\le y_{max}^j)\cap (z_{min}^j\le z_u^i\le z_{max}^j)$

其中和为多球颗粒j的局部空间的边界坐标；如 果多球颗粒i的单元球u位于上述重叠区域内，单元球u被投影至多球颗粒i的局部空间的网 格中，单元球u在多球颗粒i的局部空间中x、y和z方向上的网格号码分别为它们采用如下公式计算：

$I_{Lq}^{i,u}=Int[(q_u^i-q_{min}^i)/d]+\psi$

其中$I_{Lq}^{i,u}={\left(\begin{array}{ccc}{I}_{Lx}^{i,u}& I_{Ly}^{i,u}& I_{Lz}^{i,u}\end{array}\right)}^T.$

如图9所示，判断候选多球颗粒j的单元球v位于上述重叠区域的公式如下：

$(x_{min}^i\le x_v^j\le x_{max}^i)\cap (y_{min}^i\le y_v^j\le y_{max}^i)\cap (z_{min}^i\le z_v^j\le z_{max}^i)$

其中和是多球颗粒j的单元球v球心的全局坐标；如果多球颗粒j的单元 球v位于上述重叠区域内，单元球v被投影至多球颗粒i的局部空间的网格中，单元球v在多 球颗粒i的局部空间中x、y和z方向上的网格号码分别为它们采用如下公 式计算：

$I_{Lq}^{j,v}=Int[(q_v^j-q_{min}^i)/d]+\psi$

其中$I_{Lq}^{j,v}={\left(\begin{array}{ccc}{I}_{Lx}^{j,v}& I_{Ly}^{j,v}& I_{Lz}^{j,v}\end{array}\right)}^T,q_v^j={\left(\begin{array}{ccc}{x}_v^j& y_v^j& z_v^j\end{array}\right)}^T.$

单元球的网格号码被存储在两个数组Q_L和P_L中，其中Q_L是一个用于存储三维网格 的、大小为的三维数组，或者是一个用于存储二维网格的、大小为的二 维数组，P_L是一个用于存储三维或二维单元球编号的、大小为n_s的一维数组,这两个数组可 以在两个多球颗粒的单元球间的网格搜索中被重复使用；数组Q_L和P_L的总大小为((D+d)/d )²+n_s(二维)或((D+d)/d)³+n_s(三维)，因此第二层网格搜索的内存消耗S_L为：

步骤4：候选单元球v可以通过网格搜索来查找目标单元球u周围的邻居单元球，并 将这些邻居单元球挑选出来，其中邻居单元球位于目标单元球u所在的网格以及此网格周 围的8个(二维)或26个(三维)相邻网格之中，如图10所示，单元球和的球心位于目标 单元球周围的8个邻居网格中，这两个单元球有一定的概率与目标单元球发生接触。

目标单元球u与候选单元球v之间的接触需要被检测，如图11所示，判断单元球u与 单元球v接触的公式如下：

$\left|C_u^i{C}_v^j\right|\le (r_u^i+r_v^j)$

其中为单元球u与单元球v的球心距，和分别为单元球u和单元球v的半 径。

在二维情况下，多球颗粒i的单元球在每个网格中的平均数量如图12所示，在多球颗粒i与多球颗粒j两者的局部空间的重叠区域中，多球颗粒i的单元球 的数量为其中和L_η分别是两个局部空间的重叠区域在和η方向上的 尺寸，它们可以表达为：

其中|P_iP_j|在d到D之间变化；多球颗粒i中的每个单元球需要与多球颗粒j中的 个邻居单元球进行接触检测，其中因此在二维情况下，多 球颗粒i的单元球与多球颗粒j的单元球之间的接触检测的时间消耗为：

$\left(\begin{array}{c}{T}_E^{i,j}=O[L_{\zeta }*L_{\eta }/{\left(d\right)}^2*n_{\mathrm{CS}}^i*n_{\mathrm{CS}}^j]\\ =O[\left(\frac{(D+d)-\left|P_i{P}_j\right|}{d}\right)*\left(\frac{D+d}{d}\right)\\ =O[\left(\frac{n_s^2}{{((D+d)/d)}^2}\right)*(1-\frac{\left|P_i{P}_j\right|}{(D+d)})]\end{array},*{\left(\frac{n_s}{{((D+d)/d)}^2}\right)}^2]\right)$

以此类推，在三维情况下，多球颗粒i的单元球与多球颗粒j的单元球之间的接触 检测的时间消耗为：

$T_E^{i,j}=O[\left(\frac{n_s^2}{{\left(\right(D+d)/d)}^3}\right)*(1-\frac{\left|P_i{P}_j\right|}{(D+d)})]$

如果目标单元球u与某个邻居单元球v接触，则计算多球颗粒的接触力；反之则进 行目标单元球u与其余候选单元球之间的接触检测。

步骤5，重复步骤4直至遍历目标单元球u周围的所有候选单元球nn(u)和所有目标 单元球n_s。

步骤6，重复步骤2到步骤5直至遍历目标颗粒包络球i周围的所有候选包络球MM (i)和所有目标包络球M。

在离散元仿真的每一次计算循环中，都要运用如上所述的接触检测步骤进行多球 颗粒与多球颗粒之间的接触检测，最终识别出所有的有效接触信息，供离散元程序计算多 球颗粒的接触力。

为了验证本发明提出的两层网格搜索接触检测方法，上述的接触检测方法被嵌入 到离散元仿真程序中，并以矩形漏斗中胶囊形多球颗粒卸落过程的算例来具体说明。使用 的计算机处理器为IntelCorei5-3470CPU,主频为3.20GHz，RAM为4GB，编程工具为 VisualStudio2010。矩形漏斗和胶囊形多球颗粒的尺寸如图13a所示。仿真条件为多球颗 粒的数量M；胶囊形多球颗粒的长宽比为3:1，单个多球颗粒中单元球的数量为12，M分别为 500、1000、1500、2000、2500和3000。

下面为采用本发明的两层网格搜索接触检测方法以及采用网格搜索被直接实施 于全局空间中每个多球颗粒的每个单元球上的全局单层网格搜索接触检测方法，分别对多 球颗粒间的单元球接触进行检测；并比较上述两种方法在接触检测中的内存消耗与时间消 耗。图13b给出了四张矩形漏斗中胶囊形多球颗粒卸落过程的截图。

结果分析：

当多球颗粒的数量M在500到3000之间变化时，两层网格搜索方法与全局单层网格 搜索方法的内存消耗由图13c给出。从图13c中可以看出，对于不同的M，全局单层网格搜索 方法比两层网格搜索方法消耗了更多的内存；两层网格搜索方法与全局单层网格搜索方法 的内存消耗分别随颗粒数量的增加而增加；M对全局单层网格搜索方法的内存消耗变化的 影响要比M对两层网格搜索方法的内存消耗变化的影响更加明显。

图13d给出了两层网格搜索方法与全局单层网格搜索方法的时间消耗。从图13d中 可以看出，对于不同的M，两层网格搜索方法的时间消耗要小于全局单层网格搜索方法的时 间消耗；两层网格搜索方法与全局单层网格搜索方法的时间消耗分别随M的增加而增加。因 此，与全局单层网格搜索方法相比，两层网格搜索接触检测方法在实现接触检测的内存消 耗大幅降低的同时，维持了接触检测的时间消耗不增加且有所降低，从而有效地提高了多 球颗粒间的接触检测的计算效率。

本文中所描述的具体实施例仅仅是对本发明精神作举例说明。本发明所属技术领 域的技术人员可以对所描述的具体实施例做各种各样的修改或补充或采用类似的方式替 代，但并不会偏离本发明的精神或者超越所附权利要求书所定义的范围。