CN103235854A - 一种离散元仿真中球形颗粒与三角网格间的接触判断方法 - Google Patents

Abstract

本发明涉及一种离散元仿真中球形颗粒与三角网格间的接触判断方法，包括以下步骤：搜索目标颗粒的邻居网格以确定与之进行相交检测的边界三角形单元；结合三角形的Voronoi空间和质心坐标的符号组合以确定目标颗粒与边界三角形单元之间的初始接触信息；判定初始接触信息的有效性并排除无效的接触信息，将有效的接触信息加入接触链表，以供离散元法仿真程序计算颗粒的接触力。本发明避免了为每一类基本图形元素分别建立接触判断算法的问题；提供了直接定位发生接触区域的方法以节省接触判断的执行步骤；建立了从多重接触信息中排除无效接触信息的判定条件，最终有效地解决了离散元法仿真中球形颗粒与复杂几何边界之间的接触判断问题。

Description

一种离散元仿真中球形颗粒与三角网格间的接触判断方法

技术领域

本发明涉及基于离散元法的颗粒物料动力学仿真的领域，尤其是涉及一种离散元仿真中球形颗粒与三角网格间的接触判断方法。

背景技术

接触判断在离散元法仿真中具有极其重要的作用。通过接触判断，各目标对象间的接触情况可以通过它们前一时步的相对位置关系推导出来，这一过程对于当前时步的颗粒接触力的计算是必须的。在大型的离散元法仿真中，接触判断通常消耗了大量的计算机资源，所以如何准确并且高效地检测出所有的接触，这个问题已在相关专利和文献资料中进行了着重的阐述。通常，接触判断经历两个阶段：邻居检索和相交检测。

邻居检索的方法通常是将整个仿真区域规则地划分成相对较小的子区域，从而将单个颗粒的接触检测范围缩小到离它最近的子区域，而非整个仿真区域。通过邻居检索，接触判断问题的计算复杂度由原本的O(N²)数量级降低到较低的水平。最常用的邻居检索算法是空间排序算法和网格划分算法。对于这些方法，国际公开号WO2012034176对这些算法的性能进行了评估。

相交检测的方法主要是利用目标对象的几何形体进行空间推导。在邻居检索中判断出来的潜在接触对象将在相交检测阶段中明确地识别出来。实际上，具体的相交检测算法在很大程度上取决于颗粒的形状和几何边界的表达方法。

在三维离散元法仿真中，最常使用到的颗粒形状是球体。在一些经典的离散元法仿真应用中，球形颗粒可以较准确地表达真实的颗粒形体，例如管磨机中的球形研磨介质的运动分析，混合设备中球形物料的混合过程的数值模拟等。而对于非球形或者任意形状的颗粒，连续方程或者参数方程通常被用于来描述颗粒的表面形状。但是，大多数针对非球形颗粒开发的接触判断算法通常比基于球形颗粒的判断方法花费更多的计算时间，特别是在大型的离散元法仿真中，这种消耗是十分巨大的。

根据具体的仿真需求，几何边界表面形状的描述方法将有所不同。在一些特殊的离散元法仿真中，几何边界模型简单地由一组特殊的颗粒来表示，在整个仿真过程中这些颗粒允许具有固定的位置，或者以设定好的运动方式进行运动。基于这样的描述方式，颗粒与边界间的接触被简化成颗粒与颗粒间的接触。但是这种方法不便于描述复杂几何边界模型中的更细微的布局细节。

常规的几何边界形状也可以由一系列的基本图形元素或者它们的组合体进行表达，例如平面、圆柱体或者多面体等等，专利公开号CN102298660A和专利公开号CN1808444说明了这类边界表达方法的具体实施。球形颗粒与这些几何边界间的接触检测算法通常是非常有效的，因为这些算法充分利用了基本图形元素的规则的几何特征。更加复杂的超二次曲面形体则可以由连续方程来描述，对应的检测方法需要联立求解颗粒表面的描述方程和几何边界的描述方程。但是在许多的离散元法仿真中，仍然存在各种几何边界模型，它们具有更加复杂的表面几何形状，这些几何边界通常很难由合适的数学方程或者基本形体组合来加以表达。

三角网格可以在一定的精度范围内近似表达三维几何边界模型的表面形状。事实上，在许多大型的三维离散元法仿真中，三角网格已成为最为常用的几何边界建模方法。有许多基于三角网格的接触判断算法被设计出来。这些接触判断方法详细地讨论了球体与三角形间发生的所有的相交形式及其判定条件。在这些方法中，首先将球形颗粒的球心正交投影到三角形所在的平面上，然后将球心的投影点依次与三角形的内部区域、三角形的各边、各顶点进行相交检测。由于所有的相交形式需要顺序进行检测，所以这类算法在接触被识别出来之前，不得不被执行一系列不必要的针对于边接触或者点接触进行的检测步骤。

基于三角网格的接触判断方法必须发现并排除无效的接触。无效的接触通常发生在三角网格中相邻三角形的公共边或公共顶点上。如果颗粒同时与多个三角形发生碰撞，则多个接触信息（多重接触信息）会被检测出来。颗粒与多个三角形发生的面接触可以准确地被大多数接触判断算法作为有效的接触信息检测出来，然而多重接触中的某些接触信息必须作为无效的接触信息由接触判断方法即时地发现并有效地排除，以保证接下来的接触力的计算准确无误。例如，三角网格中两个三角形由一个公共边连接在一起，某个颗粒可能在这条公共边附近区域发生碰撞，与其中一个三角形发生面接触，而与另一个三角形发生边接触。此时，颗粒与三角形的边发生的接触应该作为一个无效的接触被排除。但是，这些无效的接触情况在已发表的专利或者文献资料中，并没有得到很好的说明。

发明内容

本发明主要是解决现有技术所存在的技术问题；提供了一种通过将颗粒-边界间的接触判断问题与边界的具体几何特征分离开来，避免了为每一类基本图形元素分别建立接触判断算法的问题；提供了直接定位发生接触区域的方法以节省接触判断的执行步骤；建立了从多重接触信息中排除无效接触信息的判定条件，最终有效地解决了离散元法仿真中球形颗粒与复杂几何边界之间的接触判断问题的一种离散元仿真中球形颗粒与三角网格间的接触判断方法。

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

一种离散元仿真中球形颗粒与三角网格间的接触判断方法，其特征在于该方法包括以下步骤：

步骤1，搜索目标颗粒的邻居网格以确定需要进行相交检测的边界三角形单元；

步骤2，结合步骤1中边界三角形单元的Voronoi空间和质心坐标的符号组合以确定目标颗粒与边界三角形单元之间的初始接触信息；

步骤3，判定初始接触信息的有效性并排除无效的接触信息，将有效的接触信息加入接触链表后采用基于离散元法仿真方法计算颗粒的接触力；

步骤4，重复步骤1到步骤3直至识别出所有的有效接触信息。

在上述的一种离散元仿真中球形颗粒与三角网格间的接触判断方法，所述目标颗粒具有球面形状，颗粒尺寸由球体半径表示。

在上述的一种离散元仿真中球形颗粒与三角网格间的接触判断方法，所述的边界三角形单元基于以下定义：

所述边界三角形单元是组成三角形网格的最小单元，其来源包括标准的基于三角网格的三维模型文件；所述边界三角形单元是平直不可弯曲的，组成三角形网格的边界三角形单元的数量和尺寸，决定了所表达的平面或者曲面的网格边界的几何精度。

在上述的一种离散元仿真中球形颗粒与三角网格间的接触判断方法，所述的步骤1中，所述搜索目标颗粒的邻居网格采用了均匀网格划分的空间邻居搜索方法，即：将整个离散元法仿真区域规则地细分成立方体网格，从而将目标颗粒的接触检测范围缩小至离它最近的若干个立方体网格，其中立方体网格的尺寸采用如下公式计算：

$G\≤\sqrt{2}{R}_{\min }$

其中，G是立方体网格的边长，R_min是仿真中使用的颗粒单元的最小半径。

在上述的一种离散元仿真中球形颗粒与三角网格间的接触判断方法，所述的步骤2中，三角形的质心坐标具有如下的定义：

P=αv₁+βv₂+γv₃=Σδ_iv_i    i∈{1,2,3}

其中点P是球形颗粒的中心C在三角形平面内的投影点；δ_i是三角形质心坐标(α,β,γ)的分量表达，它可以通过下面的式子计算：

${\mathrm{\δ}}_i=\frac{A_i}{A}=\frac{A_i}{A}\·\frac{n_T\·n_i}{\left|\right|n_T\left|\right|\left|\right|n_i\left|\right|}=\frac{\left|\right|n_i\left|\right|}{\left|\right|n_T\left|\right|}\·\frac{n_T\·n_i}{\left|\right|n_T\left|\right|\left|\right|n_i\left|\right|}=\frac{n_T\·n_i}{\left|\right|n_T|{|}^2}=\frac{h_i}{\left|\right|n_T|{|}^2}$

其中n_T是三角形的法向量，n_i是用点P与各顶点的连线将三角形划分成的三个子三角形的法向量，n_T=n_i=e₃×e₂；e₃和e₂是三角形的两条相邻的边，A_i是子三角形i的面积，A是整个三角形的面积；δ_i是点P的质心坐标分量，其符号由代表式h_i的符号决定；

三角形的Voronoi空间具有如下定义：

g₂=(C-v₁)·e₂≤0∧g₃=(C-v₁)·e₃≤0    C∈VR(v₁)

g₁=(C-v₂)·e₁≤0∧g₃’=-(C-v₂)·e₃≤0    C∈VR(v₂)

g₁’=-(C-v₃)·e₁≤0∧g₂’=-(C-v₃)·e₂≤0    C∈VR(v₃)

g_i>0∧g_i’>0∧h_i＜0    C∈VR(e_i)

其中，g₂和g₃是用于连接球心C和三角形顶点v₁的向量分别在三角形的边向量e₂和e₃上的投影，其它的g_i变量则有类似的定义；VR(v_i)表示三角的顶点v_i所在的Voronoi空间，所有的点接触VC，将在VR(v_i)中发生，其中，VC表示颗粒与三角形的顶点发生的接触；VR(e_i)则表示三角形的边e_i所在的Voronoi空间，所有的边接触EC将在VR(e_i)中发生，其中，，颗粒与三角形的边发生的接触；如果颗粒中心不位于上述两类Voronoi空间，那么颗粒必然位于三角形的内部区域，它可能与三角形的内面发生面接触FC，其中，FC表示颗粒与三角形内部区域发生的接触。

在上述的一种离散元仿真中球形颗粒与三角网格间的接触判断方法，所述的步骤2中，利用三角形的Voronoi空间和质心坐标的符号组合确定目标颗粒与边界三角形单元之间的初始接触信息的方法涉及到以下步骤：

步骤2.1，分别计算球心C的质心坐标分量δ_i所对应的符号判定表达式h_i，连接球心C和三角形顶点v_i的向量分别在三角形的边e_i上的投影表达式g_i和g_i’三个表达式的符号；

步骤2.2，通过Voronoi定义式计算颗粒中心C所处的区域；

步骤2.3，根据颗粒中心C所处的Voronoi区域，用对应的计算公式求得三角形上距离C点的最近点Q；

步骤2.4，根据公式||Q-C||²＜r²判断颗粒与边界三角形是否接触；

步骤2.5，重复步骤2.1至2.4直至所有目标颗粒遍历完毕。

在上述的一种离散元仿真中球形颗粒与三角网格间的接触判断方法，所述的步骤2.3中，根据颗粒中心C所处的Voronoi区域计算最近点Q的公式如下：

$\left\{\begin{array}{cc}Q=v_{e_i}^o+\frac{g_i}{g_i+g_i^{\′}}{e}_i& C\∈\mathrm{VR}\left(v_i\right)\\ Q=v_i& C\∈\mathrm{VR}\left(e_i\right)\\ Q=\underset{i=1}{\overset{3}{\mathrm{\Σ}}}{h}_i{v}_i& \mathrm{otherwise}\end{array}\right.$

其中，v_ei ⁰是边e_i的起始端点；如果颗粒与边界三角形发生了接触，那么接触平面的单位法向量

可以由下式计算：

$\hat{n}=\frac{Q-C}{\left|\right|Q-C\left|\right|}$

接触的嵌入量l可以由下式计算：

l＝r-||Q-C||

其中，r为目标颗粒的半径，接触点P_contact则可由下式计算：

$P_{\mathrm{contact}}=C+(r-\frac{1}{2}l)\hat{n}$

这些变量的值将作为初始接触信息保存到目标颗粒的初始接触链表中，以便进一步判定初始接触信息的有效性。

在上述的一种离散元仿真中球形颗粒与三角网格间的接触判断方法，所述的步骤2中，考虑到颗粒极可能同时与多个相邻的边界三角形发生接触，这些多重接触中可能包含了无效的接触信息；假设颗粒与两个相邻的边界三角形T₁和T₂同时发生接触，并产生了两个初始接触信息C₁和C₂，基于C₂的接触类型给出判定C₂的有效性的条件：

条件1：若C₂为面接触，则C₂直接判定为有效的面接触信息；

条件1：若C₂为边接触，则在两种情况下，C₂为有效的边接触信息，它们分别是：

A、T₁与T₂共顶；

B、C₁为点接触且T₁与T₂共边；

否则C₂为无效的边接触信息；

条件3：若C₂为点接触，则只有当C₁不是点接触且T₁与T₂共边的情况下，C₂才为有效的点接触信息，否则C₂为无效的点接触信息。

在执行初始接触信息的有效性判定时，首先依次从目标颗粒的初始接触链表中取出初始接触信息，上述判定条件进行有效性验证；被判定为无效的多重接触信息将从初始接触信息中删除，而判定为有效的接触信息将加入到目标颗粒所属的接触链表中。

总体来看，在离散元法仿真的每一次计算循环中，都要针对每一个颗粒单元运用如上所述的接触判断步骤进行处理，最终识别出所有的有效接触信息，以供离散元法程序计算颗粒当前的接触力。

因此，本发明具有如下优点：通过将颗粒-边界间的接触判断问题与边界的具体几何特征分离开来，避免了为每一类基本图形元素分别建立接触判断算法的问题；提供了直接定位发生接触区域的方法以节省接触判断的执行步骤；建立了从多重接触信息中排除无效接触信息的判定条件，最终有效地解决了离散元法仿真中球形颗粒与复杂几何边界之间的接触判断问题。

附图说明

图1是三角形平面内一点的质心坐标表达方法。

图2是三角形质心坐标分量的符号组合。

图3是三角形Voronoi空间的划分。

图4是接触判断方法的数据结构。

图5是接触判断方法的主流程图。

图6是球形颗粒与三角形单元间的相交检测流程图。

图7是初始接触信息的有效性判定流程图。

图8a是拌和装置的三角网格化CAD模型。

图8b是拌和部件的三角网格化CAD模型。

图8c是拌和部件的主视图。

图8d是拌和部件的俯视图。

图8e是拌和部件的左视图。

具体实施方式

下面通过实施例，并结合附图，对本发明的技术方案作进一步具体的说明。

实施例：

下面结合附图进一步说明本发明的实施方式。

首先，介绍一下本发明的具体方法步骤：

本发明提供了一种用于离散元法仿真的球形颗粒与三角网格边界间的接触判断方法，包括以下步骤：

（1）搜索目标颗粒的邻居网格以确定需要进行相交检测的边界三角形单元；

（2）结合三角形的Voronoi空间和质心坐标的符号组合以确定目标颗粒与边界三角形单元之间的初始接触信息；

（3）判定初始接触信息的有效性并排除无效的接触信息，将有效的接触信息加入接触链表，以供离散元法仿真程序计算颗粒的接触力。

本发明提出的接触判断方法中，涉及到的颗粒具有球面形状，颗粒尺寸由球体半径表示。几何边界模型是经过三角网格化的三维CAD模型，这些模型的来源包括标准的基于三角网格的三维模型文件，例如OBJ（Wavefront’s object files）文件、STL（Stereo Lithography files）文件等，或者经过表面网格划分工具处理后得到的三角网格化模型文件。

几何边界的表面由许多平直不可弯曲的边界三角形单元近似表达（几何误差由边界三角形单元的数量和尺寸决定，一般而言，数量越大、尺寸越小，则表达越精确）；几何边界模型可以具有封闭、非封闭的几何形状，可以表达规则的几何实体，或者近似表达复杂的曲面形体。

本发明提供的接触判断方法中，采用了均匀网格划分的空间邻居搜索方法。该方法将整个离散元法仿真区域（求解区域）规则地细分成相对较小的立方体网格，从而将目标颗粒的接触检测范围将缩小至离它最近的若干个立方网格，进而降低了接触判断问题的计算复杂度；其中立方体网格的尺寸通常按如下公式计算：

$G\≤\sqrt{2}{R}_{\min }$

其中，G是立方体网格的边长，R_min是仿真中使用的颗粒单元的最小半径。

本发明使用了三角形的质心坐标，它具有如下的定义：

P=αv₁+βv₂+γv₃=Σδ_iv_i    i∈{1,2,3}

其中点P是球形颗粒的中心C在三角形平面内的投影点，如图1所示；δ_i是三角形质心坐标(α,β,γ)的分量表达，它可以通过下面的式子计算：

${\mathrm{\δ}}_i=\frac{A_i}{A}=\frac{A_i}{A}\·\frac{n_T\·n_i}{\left|\right|n_T\left|\right|\left|\right|n_i\left|\right|}=\frac{\left|\right|n_i\left|\right|}{\left|\right|n_T\left|\right|}\·\frac{n_T\·n_i}{\left|\right|n_T\left|\right|\left|\right|n_i\left|\right|}=\frac{n_T\·n_i}{\left|\right|n_T|{|}^2}=h_i$

其中n_T是三角形的法向量，n_i是用点P与各顶点的连线将三角形划分成的三个子三角形的法向量，n_T=n_i=e₃×e；根据点P与三角形的相对位置（三角形内、三角形边上或者三角形以外区域），点P的质心坐标分量δ_i将具有不同的符号组合，如图2所示。

同时，本发明还涉及到三角形的Voronoi空间，它具有如下定义：

g₂=(C-v₁)·e₂≤0∧g₃=(C-v₁)·e₃≤0    C∈VR(v₁)

g₁=(C-v₂)·e₁≤0∧g₃’=-(C-v₂)·e₃≤0    C∈VR(v₂)

g₁’=-(C-v₃)·e₁≤0∧g₂’=-(C-v₃)·e₂≤0    C∈VR(v₃)

g_i>0∧g_i’>0∧h_i＜0    C∈VR(e_i)

其中，g₂和g₃是用于连接球心C和三角形顶点v₁的向量分别在三角形的边向量e₂和e₃上的投影，其它的g_i变量则有类似的定义；如图3所示，VR(v_i)表示三角的顶点v_i所在的Voronoi空间，所有的点接触（VC，颗粒与三角形的顶点发生的接触）将在VR(v_i)中发生；VR(e_i)则表示三角形的边e_i所在的Voronoi空间，所有的边接触（EC，颗粒与三角形的边发生的接触）将在VR(e_i)中发生；如果颗粒中心不位于上述两类Voronoi空间，那么颗粒必然位于三角形的内部区域，它可能与三角形的内面发生面接触（FC），如图3所示。

本发明提出的结合三角形的Voronoi空间和质心坐标的符号组合确定目标颗粒与边界三角形单元之间的初始接触信息的方法涉及到以下步骤：

（1）计算上述的h_i，g_i和g_i’变量的符号；

（2）通过三角形的Voronoi定义式计算颗粒中心C所处的区域；

（3）根据颗粒中心C所处的Voronoi区域，用对应的计算公式求得三角形上距离C点的最近点Q；

（4）根据公式||Q-C||²＜r²判断颗粒与边界三角形是否接触。

根据颗粒中心C所处的Voronoi区域计算最近点Q的公式如下：

$\left\{\begin{array}{cc}Q=v_{e_i}^o+\frac{g_i}{g_i+g_i^{\′}}{e}_i& C\∈\mathrm{VR}\left(v_i\right)\\ Q=v_i& C\∈\mathrm{VR}\left(e_i\right)\\ Q=\underset{i=1}{\overset{3}{\mathrm{\Σ}}}{h}_i{v}_i& \mathrm{otherwise}\end{array}\right.$

其中，v_ei ⁰是边e_i的起始端点；如果颗粒与边界三角形发生了接触，那么接触平面的单位法向量

可以由下式计算：

$\hat{n}=\frac{Q-C}{\left|\right|Q-C\left|\right|}$

接触的嵌入量l可以由下式计算：

l＝r-||Q-C||

其中，r为目标颗粒的半径，接触点P_contact则可由下式计算：

$P_{\mathrm{contact}}=C+(r-\frac{1}{2}l)\hat{n}$

这些变量的值将作为初始接触信息保存到目标颗粒的初始接触链表中，以便进一步判定初始接触信息的有效性。

本发明提出的接触判断方法考虑了颗粒的相邻边界三角形间的相对位置关系，其中包括：共面共边（CP&ES）、共面共点（CP&VS）、异面共边（NP&ES）、异面共点（NP&VS）、异面无关联或者共面无关联（Unlinked）六种情况。

由于颗粒极可能同时与多个相邻的边界三角形发生接触，这些多重接触中可能包含了无效的接触信息；假设颗粒与两个相邻的边界三角形T₁和T₂同时发生接触，并产生了两个初始接触信息C₁和C₂，下面的分析方法给出了判定C₂的有效性的条件：

（1）如果C₂为面接触，则有效性判定条件如下表所示，表中“-”代表“不可能”或者“无影响”；“√”代表“有效”而“×”代表“无效”。

表1C₂为面接触时的有效性判定

（2）如果C₂为边接触，则有效性判定条件如下表所示：

表2C₂为边接触时的有效性判定

（3）如果C₂为点接触，则有效性判定条件如下表所示：

表3C₂为点接触时的有效性判定

下面，结合上述的本发明的方法步骤，应用到具体的实施例中。在实验中，一旦颗粒同时与网格边界中的多个三角形单元发生碰撞，多重接触信息将作为初始接触信息被记录下来。各种类型的初始接触信息（FC,EC和VC）将通过表1至表3中所列出的有效性判定条件进行有效性检测，从而产生一组有效的接触信息。通常，这些有效的接触信息被添加到一个链表类型的数据结构中，这个链表被称为“接触链表”。这个接触信息的有效性的检测过程将一直循环下去，直到所有的有效接触信息都被检测出来。在实现这个检测过程中，需要解决以下几个主要问题：

（1）需要设计一套数据结构以便于高效地执行表1至表3中列出的判定条件。因为表中的有效性判定条件与待检测的接触信息的接触类型相关，所以整个数据结构可以根据不同的接触类型划分成三个部分。

（2）离散元法仿真中，对象间允许有少量的交叠，所以对象间的接触通常需要经历数个仿真时间步长。在目前的接触判断算法中，一个接触信息在它第一次被判定为有效的接触信息时，被作为一个“新接触”加载到接触链表中，然后作为一个“已存在的接触”记录在接触链表中存在若干时间步长，最后随着这个接触的消失而从链表中删除。由于对新接触与已存在的接触的处理方式有所不同，所以接触链表中的接触信息需要在被加载到链表中之前标明它是否为新接触。

（3）颗粒与边界发生碰撞并反弹出去，两者间的接触消失。这些消失的接触对应的接触信息也需要即时地从接触链表中清除。

为了解决第一个问题，被判定为有效的接触信息可以被加载到如图4所示的数据结构中。其中，三个临时链表被创建出来用于保存三种类型的初始接触信息，它们被称为“初始接触链表”。在图4中，这些初始接触链表中的元素由相互联接的正方形表达（用方格填充），标记为V_Pool，E_Pool和F_Pool的初始接触接触链表分别用于保存初始的点接触信息、初始的边接触信息和初始的面接触信息。被判定为有效的接触信息将被加载到对应颗粒的接触链表中，图4中空白矩形表示颗粒的接触链表中的有效的接触信息。在全局范围内，离散元法仿真框架为仿真中的所有发生了接触的颗粒单元维护了一个“全局接触链表”，图4中斜方格填充的矩形代表了全局接触链表中的元素，其中保存了各颗粒的接触链表的表头地址。

在上述数据结构的基础上，第二个问题很容易得到解决。一个接触信息至少包含相互接触的两个对象（颗粒与三角形单元，或者颗粒与颗粒）的索引信息。如果每个对象使用两个参数来标记自己，例如对象的类型和对象的序号，那么一个接触信息可以由两个对象的四个索引参数确定。如果两个接触信息C_i和C_j中的对象类型和序号分别相同，那么这两个接触是同一个接触信息。按照这样的方法，通过遍历某颗粒的接触链表可以查找出该链表中是否存在与给定接触信息相同的接触。如果可以找到相同的接触，那么直接将当前的接触信息更新到这个已存在的接触信息中，否则，需要建立一个新的接触信息并加载到颗粒的接触链表中。

至于第三个问题，可以在接触信息中添加一个布尔类型的变量“OUT_OF_DATA”。在仿真的每个计算循环中，通过检查这个变量的值可以确定对应的接触是否已在当前时步消失。具体的过程可以由图5表示。

图5的环节A中，所有记录的接触信息中的OUT_OF_DATA变量将被置为真值。在环节B中，只有那些被判定为有效的接触信息中的OUT_OF_DATA变量的值被修改为假值。当算法流程到达环节C，任何OUT_OF_DATA仍然为真值的接触信息将作为已消失的接触从接触链表中删除。通过这样的方法，当前已消失的接触记录将准确地清除。其中环节B中详细的流程由图6说明。

接触判断算法中采用的邻居检索算法是均匀网格划分方法。当算法流程到达环节D，如图6所示，可以通过前面描述的相交检测条件识别当前仿真计算循环中，已发生的所有初始接触信息。当颗粒的所有初始接触信息被加载到对应的初始接触链表中之后，它们将在环节E中，通过表1至表3中列出的有效性判定条件进行有效性检测。具体的判定流程由图7表示。在环节F，有效的接触信息被加载到颗粒的接触链表中，同时，有效接触信息的OUT_OF_DATA变量将被修改为假值。最后，三种类型的初始接触链表将在环节G中被清空。

在图7中，标记V_C，E_C和F_C分别代表点接触VC，边接触EC和面接触FC；W则表示与颗粒发生接触的边界三角形单元。

为了验证本发明提出的接触判断方法，上述的接触判断方法被写成计算机程序并嵌入到离散元法仿真程序中，下面使用一个物料混合的算例来具体说明。

所图8a所示的静态拌和装置由四层拌和部件以相互交错的方式搭建而成，其中，拌和部件的结构如图8b所示，其尺寸由图8c、图8d和图8e给出。拌和部件的CAD模型表面由多个边界三角形单元组成，其中各滑坡（slope）由三个切面（facet）组成；各切面又由两个边界三角形组成。拌和装置的四周由挡板（baffle board）包围，各挡板分别由两个边界三角形单元组成。这个边界模型是在Pro/E软件中进行建模的，然后将模型保存为OBJ文件，通过离散元法仿真程序导入这个OBJ文件并在仿真场景中重绘出这个边界模型。此边界模型的表面共由56片三角形单元组成，其中有52条公共边，16条自由边和12个公共点，因此选择它作为验证接触判断方法的几何边界模型是合适的。另外，算例中采用的物料样本是由蓝白两种球形颗粒各1500颗组成。

仿真初始时刻，颗粒物料沉降于拌和装置顶端的料斗中；仿真开始后，离散元法仿真程序进行循环迭代计算，在每一个计算循环中，首先调用了本发明提出的接触判断方法，根据前一迭代时步中颗粒的位置来识别颗粒与三角网格边界间发生的接触，然后根据各颗粒所属的接触链表中包含的有效接触信息进行颗粒的接触力的更新，进而计算出各颗粒当前时步的速度与位置，如此往复；当仿真结束，颗粒被均匀地混合在一起并沉降在混合设备的底部箱体。整个仿真过程能正确地反映物料在混合设备中的运动过程，本发明提出的接触判断方法正确地识别了球形颗粒与几何边界间发生的各类接触。

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

Claims (8)

1.一种离散元仿真中球形颗粒与三角网格间的接触判断方法，其特征在于该方法包括以下步骤：

步骤1，搜索目标颗粒的邻居网格以确定需要进行相交检测的边界三角形单元；

步骤2，结合步骤1中边界三角形单元的Voronoi空间和质心坐标的符号组合以确定目标颗粒与边界三角形单元之间的初始接触信息；

步骤3，判定初始接触信息的有效性并排除无效的接触信息，将有效的接触信息加入接触链表后采用基于离散元法仿真方法计算颗粒的接触力；

步骤4，重复步骤1到步骤3直至识别出所有的有效接触信息。

2.根据权利要求1所述的一种离散元仿真中球形颗粒与三角网格间的接触判断方法，其特征在于，所述目标颗粒具有球面形状，颗粒尺寸由球体半径表示。

3.根据权利要求1所述的一种离散元仿真中球形颗粒与三角网格间的接触判断方法，其特征在于，所述的边界三角形单元基于以下定义：

所述边界三角形单元是组成三角形网格的最小单元，其来源包括标准的基于三角网格的三维模型文件；所述边界三角形单元是平直不可弯曲的，组成三角形网格的边界三角形单元的数量和尺寸，决定了所表达的平面或者曲面的网格边界的几何精度。

4.根据权利要求1所述的一种离散元仿真中球形颗粒与三角网格间的接触判断方法，其特征在于，所述的步骤1中，所述搜索目标颗粒的邻居网格采用了均匀网格划分的空间邻居搜索方法，即：将整个离散元法仿真区域规则地细分成立方体网格，从而将目标颗粒的接触检测范围缩小至离它最近的若干个立方体网格，其中立方体网格的尺寸采用如下公式计算：

$G\≤\sqrt{2}{R}_{\min }$

其中，G是立方体网格的边长，R_min是仿真中使用的颗粒单元的最小半径。

5.根据权利要求1所述的一种离散元仿真中球形颗粒与三角网格间的接触判断方法，其特征在于，所述的步骤2中，三角形的质心坐标具有如下的定义：

P=αv₁+βv₂+γv₃=Σδ_iv_i    i∈{1,2,3}

其中点P是球形颗粒的中心C在三角形平面内的投影点；δ_i是三角形质心坐标(α,β,γ)的分量表达，它可以通过下面的式子计算：

${\mathrm{\δ}}_i=\frac{A_i}{A}=\frac{A_i}{A}\·\frac{n_T\·n_i}{\left|\right|n_T\left|\right|\left|\right|n_i\left|\right|}=\frac{\left|\right|n_i\left|\right|}{\left|\right|n_T\left|\right|}\·\frac{n_T\·n_i}{\left|\right|n_T\left|\right|\left|\right|n_i\left|\right|}=\frac{n_T\·n_i}{\left|\right|n_T|{|}^2}=\frac{h_i}{\left|\right|n_T|{|}^2}$

其中n_T是三角形的法向量，n_i是用点P与各顶点的连线将三角形划分成的三个子三角形的法向量，n_T=n_i=e₃×e₂；e₃和e₂是三角形的两条相邻的边，A_i是子三角形i的面积，A是整个三角形的面积；δ_i是点P的质心坐标分量，其符号由代表式h_i的符号决定；

三角形的Voronoi空间具有如下定义：

g₂=(C-v₁)·e₂≤0∧g₃=(C-v₁)·e₃≤0    C∈VR(v₁)

g₁=(C-v₂)·e₁≤0∧g₃’=-(C-v₂)·e₃≤0    C∈VR(v₂)

g₁’=-(C-v₃)·e₁≤0∧g₂’=-(C-v₃)·e₂≤0    C∈VR(v₃)

g_i>0∧g_i’>0∧h_i＜0    C∈VR(e_i)

其中，g₂和g₃是用于连接球心C和三角形顶点v₁的向量分别在三角形的边向量e₂和e₃上的投影，其它的g_i变量则有类似的定义；VR(v_i)表示三角的顶点v_i所在的Voronoi空间，所有的点接触VC，将在VR(v_i)中发生，其中，VC表示颗粒与三角形的顶点发生的接触；VR(e_i)则表示三角形的边e_i所在的Voronoi空间，所有的边接触EC将在VR(e_i)中发生，其中，，颗粒与三角形的边发生的接触；如果颗粒中心不位于上述两类Voronoi空间，那么颗粒必然位于三角形的内部区域，它可能与三角形的内面发生面接触FC，其中，FC表示颗粒与三角形内部区域发生的接触。

6.根据权利要求5所述的一种离散元仿真中球形颗粒与三角网格间的接触判断方法，其特征在于，所述的步骤2中，利用三角形的Voronoi空间和质心坐标的符号组合确定目标颗粒与边界三角形单元之间的初始接触信息的方法涉及到以下步骤：

步骤2.1，分别计算球心C的质心坐标分量δ_i所对应的符号判定表达式h_i，连接球心C和三角形顶点v_i的向量分别在三角形的边e_i上的投影表达式g_i和g_i’三个表达式的符号；

步骤2.2，通过Voronoi定义式计算颗粒中心C所处的区域；

步骤2.3，根据颗粒中心C所处的Voronoi区域，用对应的计算公式求得三角形上距离C点的最近点Q；

步骤2.4，根据公式||Q-C||²＜r²判断颗粒与边界三角形是否接触；

步骤2.5，重复步骤2.1至2.4直至所有目标颗粒遍历完毕。

7.根据权利要求6所述的一种离散元仿真中球形颗粒与三角网格间的接触判断方法，其特征在于，所述的步骤2.3中，根据颗粒中心C所处的Voronoi区域计算最近点Q的公式如下：

$\left\{\begin{array}{cc}Q=v_{e_i}^o+\frac{g_i}{g_i+g_i^{\′}}{e}_i& C\∈\mathrm{VR}\left(v_i\right)\\ Q=v_i& C\∈\mathrm{VR}\left(e_i\right)\\ Q=\underset{i=1}{\overset{3}{\mathrm{\Σ}}}{h}_i{v}_i& \mathrm{otherwise}\end{array}\right.$

其中，v_ei ^o是边e_i的起始端点；如果颗粒与边界三角形发生了接触，那么接触平面的单位法向量

可以由下式计算：

$\hat{n}=\frac{Q-C}{\left|\right|Q-C\left|\right|}$

接触的嵌入量l可以由下式计算：

l＝r-||Q-C||

其中，r为目标颗粒的半径，接触点P_contact则可由下式计算：

$P_{\mathrm{contact}}=C+(r-\frac{1}{2}l)\hat{n}$

这些变量的值将作为初始接触信息保存到目标颗粒的初始接触链表中，以便进一步判定初始接触信息的有效性。

8.根据权利要求6所述的一种离散元仿真中球形颗粒与三角网格间的接触判断方法，其特征在于，所述的步骤2中，考虑到颗粒极可能同时与多个相邻的边界三角形发生接触，这些多重接触中可能包含了无效的接触信息；假设颗粒与两个相邻的边界三角形T₁和T₂同时发生接触，并产生了两个初始接触信息C₁和C₂，基于C₂的接触类型给出判定C₂的有效性的条件：

条件1：若C₂为面接触，则C₂直接判定为有效的面接触信息；

条件1：若C₂为边接触，则在两种情况下，C₂为有效的边接触信息，它们分别是：

A、T₁与T₂共顶；

B、C₁为点接触且T₁与T₂共边；

否则C₂为无效的边接触信息；

条件3：若C₂为点接触，则只有当C₁不是点接触且T₁与T₂共边的情况下，C₂才为有效的点接触信息，否则C₂为无效的点接触信息。

Images