离散元仿真中挖掘机铲斗结构表面的单元球表示方法探究

方自强

离散元仿真中挖掘机铲斗结构表面的单元球表示方法探究

方自强

（江汉大学 机电与建筑工程学院，湖北 武汉 430056）

提出了一种基于单元球的铲斗结构表面的多精度几何表示方法。该方法针对铲斗的斗底板、侧腹板、侧刀板、前刀板和斗齿5个部位的几何特征和功能要求，实施基于单元球的多精度几何表示。对比分析不同精度下单元球表达的铲斗结构表面，发现离散元仿真时间会随铲斗结构表面的几何表示精度的提高而增加。因而，铲斗结构表面的单元球表示需综合考虑仿真精准性与仿真时间的需求。对比分析不同精度下单元球与三角面片表示方法，证实了单元球表示方法具备有效表示铲斗结构表面的局部几何细节的能力。相比于三角面片表示方法，单元球表示方法有望在铲斗挖掘过程离散元仿真中花费更少的计算时间。

挖掘机铲斗；结构表面；几何表示；离散元仿真；单元球

颗粒离散元法[1]主要用于模拟工业过程中颗粒与结构表面间的碰撞行为，如挖掘机铲斗对颗粒散料的挖掘过程[2-3]。通过在一个足够小的时步内更新颗粒的位移、速度、加速度和接触力等参数[4]，离散元法可以获取颗粒与铲斗表面间相互碰撞过程中的接触速度、碰撞能量等信息[5]。

挖掘机铲斗的结构表面表示是离散元仿真的一项关键技术[6]，它直接影响着散料挖掘过程离散元仿真的接触检测效率以及接触力求解的精准性。目前，三角面片是离散元仿真中使用最为广泛的一类结构表面表示方法[7]。利用ICEM CFD等网格划分工具对铲斗结构表面进行三角面片划分，这些三角面片可有效描述铲斗结构表面的几何形状。经三角面片表示，颗粒与铲斗表面间的接触检测可以转化为颗粒与三角面片间的接触检测。于是，在一个颗粒与一个三角面片间的相交求解中，该颗粒需要与此三角面片的面、边和点相继进行相交求解，从而导致了巨大的仿真计算开销。

由于单元球与单元球间的接触检测仅涉及一个中心距和两个半径和之间的比较[6]，因而单元球与单元球间的接触检测效率有望远高于单元球与三角面片之间的接触检测效率。因此，本文以离散元仿真中常见的挖掘机铲斗结构表面为对象[8-9]，开展不同分辨率下的单元球表示方法的探究。

1 基于单元球的铲斗结构表面的多精度几何表示方法

1.1 铲斗结构表面几何表示精度的解析

在铲斗结构表面的几何表示之前，首先需要根据离散元仿真对铲斗结构表面的几何特征和功能要求[10]，进行几何表示精度的解析。如图1所示，挖掘机铲斗由斗底板、侧腹板、侧刀板、前刀板和斗齿共5个部件组成[11-12]。铲斗的挖掘作业包括挖掘、回转、卸料和返回四个过程。在挖掘与回转过程中，斗齿[13]对颗粒堆表面进行剥离作业，因而承受了颗粒的绝大部分接触力，而侧刀板和前刀板也与颗粒发生直接碰撞。相对而言，斗底板和侧腹板主要用于颗粒散料的装载[14-15]，因而较少受到颗粒的冲击和磨擦。

图1 铲斗的几何模型与尺寸

为了利用离散元仿真有效分析挖掘机作业中颗粒对铲斗结构表面的冲击和磨擦[16-17]，需要对铲斗结构表面的不同部位进行不同精度的单元球几何表示。例如，斗齿表面与颗粒的接触最为频繁，该部位需要进行高精度的几何表示；侧刀板和前刀板也要长时间承受颗粒的作用力，它们的几何表示精度要求也相对较高；斗底板和侧腹板较少受到颗粒的冲击和磨擦，对这两个部位进行低精度的几何表示。于是，在本文的铲斗结构表面的几何表示中，斗底板和侧腹板的表示精度最低，侧刀板和前刀板的表示精度是它们的两倍，斗齿的表示精度是它们的三倍。

在采用单元球对铲斗结构表面进行几何表示的过程中，决定铲斗结构表面几何表示精度的参数主要包括相邻单元球之间的中心距离和单元球的直径d。因而在相邻单元球之间的保持一定的前提下，有：

式中：1为斗底板和侧腹板的单元球直径，mm；2为斗底板和侧腹板的单元球直径，mm；1为斗齿的单元球直径，mm。

此外，在d保持一定的前提下，有：

式中：1为斗底板和侧腹板的相邻单元球之间的中心距，mm；2为侧刀板和前刀板的相邻单元球之间的中心距，mm；3为斗齿的相邻单元球之间的中心距，mm。

1.2 铲斗结构表面的单元球多精度几何表示

根据上述解析，依据式（1）所给出的通过变化d来改变铲斗结构表面的不同部位几何表示精度的方式，分别对图1（a）中的斗底板、侧腹板、侧刀板、前刀板和斗齿五个部位进行几何重构和三角面片划分[18]。

首先，利用三维CAD软件，对上述结构表面进行表面内偏移操作，完成几何重构。当表示斗底板和侧腹板、侧刀板和前刀板、斗齿结构表面的单元球直径分别被设置为15 mm、10 mm、5 mm时，它们各自对应的内偏移距离分别为7.5 mm、5 mm、2.5 mm。如图2（a）和（b）所示，重构后的斗底板和侧腹板分别由曲面和平面来表示。如图2（c）、（d）和（e）所示，重构后的侧刀板、前刀板和斗齿几何形状则由相应几何体的外壳表面来表示。

然后，利用ICEM CFD或Gambit等软件，对几何重构后的结构表面进行三角面片划分。在结构表面的三角面片划分中，三角面片的最大尺寸依据各自单元球的直径来确定。在结构表面的单元球几何表示中，要求单元球之间不存在缝隙，能够相互重叠，于是三角面片划分的最大尺寸需要设定得小于或等于单元球的直径。如图2（a-1）、（b-1）、（c-1）、（d-1）、（e-1）所示，重构后的斗底板和侧腹板、侧刀板和前刀板、斗齿的结构表面被分别划分为最大尺寸为12 mm、8 mm、4 mm的三角面片。

最后，以各个三角面片的顶点为依据，在斗底板和侧腹板、侧刀板和前刀板、斗齿对应的三角面片顶点上分别添加直径为15 mm、10 mm、5 mm的单元球，可实现如图3（a）所示的铲斗结构表面单元球多精度几何表示[6]。

需要说明的是，依据式（2）所给出的通过变化相邻单元球之间的来改变铲斗结构表面的不同部位几何表示精度的方式，同样可以实现铲斗结构表面的单元球多精度几何表示。其中，相邻单元球之间的对应于几何重构后结构表面三角面片划分中三角面片的尺寸大小。

2 铲斗结构表面单元球表示方法的性能分析

2.1 三种不同的最大中心距下单元球表示方法的性能分析

通过上述基于单元球的铲斗结构表面多精度几何表示方法，本文在采用相同单元球粒径d的前提下，通过变换单元球间的中心距离，开展三种不同最大中心距max下单元球表示方法的性能分析。

依据式（1），用于表示斗底板和侧腹板、侧刀板和前刀板、斗齿结构表面的d分别为15 mm、10 mm、5 mm。在三种单元球表示精度中，单元球间的max分别为0.8d、0.6d和0.4d，具体相关尺寸由表1列出。

图2 铲斗模型的几何重构与三角面片划分

表1 三种不同表示精度下单元球间的最大中心距zimax

图3分别给出了max为0.8d、0.6d和0.4d的条件下，采用多粒径单元球表示的铲斗结构表面的几何形状。随着max的逐渐减小，铲斗结构表面的几何表面更加平滑。这说明铲斗结构表面的几何表示精度随着max的减小而增加，更有益于铲斗对颗粒挖掘过程离散元仿真精度的提升。

此外，上述三种max所表示的铲斗结构表面的单元球数量分别为91342、161317、354328。在离散元仿真中，颗粒与结构表面间接触检测的邻居检索时间复杂度[18]可表示为：

式中：N为颗粒的总数量；N为单元球的总数量；为一个与颗粒尺寸和分布相关的经验参数；{1,2}Î[1, 2]。

于是，随着单元球数量的增加，颗粒与结构表面间接触检测的时间消耗会呈几何级数的增长，不利于离散元仿真效率的提高。因此，在进行铲斗结构表面的单元球表示时，需综合考虑仿真时间与仿真精准性的需求，以确定单元球间的max。

2.2 三种不同的单元球直径下单元球表示方法的性能分析

类似地，在保持单元球间max一致的前提下，通过改变d，开展三种不同d下单元球表示方法的性能分析。依据式（1），用于表示铲斗结构表面的单元球之间的max均为0.8d，三组用于表示铲斗结构表面的单元球直径1、2和3分为等于D、0.8D和0.6D，其中D为图2中铲斗结构边界的内偏移距离。于是，得出d和max值如表2所示。

在单元球直径d为D、0.8D和0.6D条件下，采用多粒径单元球表示的铲斗结构表面的几何形状分别如图4所示。当单元球之间的max为0.8d时，铲斗结构表面的几何表示精度随着d的减小而显著提高。因此， d的减小有益于铲斗结构表面的几何表示精度的提升。

如图4所示，用于表示铲斗结构表面的单元球数量分别为91342、148957、279848。根据式（3）给出的颗粒与结构表面间接触检测的时间复杂度，离散元仿真中颗粒与结构表面间接触检测的时间消耗会随d的降低而显著增大。因此，在铲斗结构表面的单元球表示中，d的选择也是决定离散元仿真效率与铲斗结构表面几何表示精度的重要参数之一。

图3 三种不同最大中心距下单元球表示的铲斗结构表面

表2 三种不同表示精度下单元球的粒径di和最大中心距zimax

图4 三种不同di下单元球表示的铲斗结构表面

3 单元球与三角面片表示方法的性能比较及分析

3.1 几何表示精度的比较

为了对比分析单元球与三角面片[7]表示方法的性能，本文以图1所示的铲斗各部位的壳体表面为对象，利用ICEM CFD工具，分别对斗底板和侧腹板、侧刀板和前刀板、斗齿进行三角面片划分，其中三角面片尺寸以表1所列出的max为依据。于是，对应于图3中由单元球表示的铲斗结构表面，图5给出了由三角面片表示的铲斗结构表面。

在图5中，三角面片的数量分别为227844、403008和909658。随着三角面片尺寸的减小，三角面片所表示的铲斗结构表面的几何精度显著提升，这与单元球表示的铲斗结构表面的几何精度随着的变化规律一致。

此外，通过对比同一精度下由单元球和由三角面片表示的铲斗结构表面的几何形状，可发现单元球能够有效表示铲斗结构表面的局部结构细节，因而能够有效表示上述不同精度下铲斗结构表面的几何形状。

图5 三种不同精度下三角面片表示的铲斗结构表面

3.2 接触检测时间消耗的理论分析

根据式（3）所给出的时间复杂度，在颗粒数量相同的前提下，单元球与三角面片表示方法对应的颗粒与结构表面间接触检测的邻居检索的时间复杂度在同一个数量级。

对于相交求解，一个单元球与一个球形颗粒间的相交求解条件为：

参照文献[7,18]，一个三角面片与一个球形颗粒间的相交求解涉及三角面片的一个面、三条边和三个角与球形颗粒间的连续求解操作。于是，本文将一个单元球与一个球形颗粒间的相交求解以及一个三角面片与一个球形颗粒间的相交求解[7]的计算操作符号数量进行对比，如表3所示。

在以C++为语言的离散元仿真的算术运算中，“+”、“-”、“´”、“/”等计算操作符号为基本符号，其计算时间消耗较小，而“||-||”、“|-|”、“Ö”、“acos(-)”、“min(-)”等计算操作符号是以“+”、“-”、“´”、“/”为基础的复合计算操作符号[19]，所消耗的计算时间十分巨大。

如表3所示，一个单元球与一个球形颗粒间的相交求解仅涉及一个“+”、一个“||-||”和一个“>”。相较而言，一个三角面片与一个球形颗粒间的相交求解所涉及的“+”、“-”、“´”、“/”等计算操作符号是一个单元球与一个球形颗粒间的相交求解对应操作符号数量的若干倍，并且还涉及多个“||-||”、“|-|”、“Ö”、“acos(-)”、“min(-)”等复合计算操作符号。因而，单元球与颗粒间的相交求解时间消耗远小于三角面片[7,18]与颗粒间的相交求解时间消耗。可以期望的是，采用单元球表示铲斗结构表面的几何形状会在离散元仿真中花费更少的计算时间，更有利于离散元仿真计算效率的提高。

表3 相交求解的计算操作符号数量

4 结论

（1）本文提出了一种基于单元球的铲斗结构表面的多精度几何表示方法，该方法针对铲斗不同部位的几何特征和功能要求，实施不同的单元球粒径和不同最大中心距的多精度几何表示。

（2）对比分析三种不同精度下单元球表达的铲斗结构表面，发现铲斗结构表面的几何表示精度以及颗粒与结构表面间接触检测的时间消耗均随着的减小而增加。因而，铲斗结构表面的单元球表示需综合考虑仿真精准性与仿真时间的需求，以确定单元球间的max。

（3）对比分析三种不同精度下单元球与三角面片表示方法，可知单元球具备有效表示铲斗结构表面的局部几何细节的能力。相比于三角面片方法，单元球表示方法有望在离散元仿真中花费更少的计算时间，更有利于仿真效率的提高。

[1]方自强，吴洪亮，周享楠. 球磨机DEM 仿真中接触模型的精准性分析[J]. 机械，2019，46（2）：17-24.

[2]C. J. Coetzee，A. H. Basson，P. A. Vermeer. Discrete and continuum modelling of excavator bucket filling [J]. J. Terramechanics，2007（44）：177-186.

[3]滕忠斌. 双斗轮取料机改造分析[J]. 机械，2016，43（11）：23-29.

[4]胡国明. 颗粒系统的离散元素法分析仿真[M]. 武汉：武汉理工大学出版社，2010.

[5]J. Maciejewski，A. Jarzebowski. Laboratory optimization of the soildigging process [J]. J. Terramechanics，2002（39）：161-179.

[6]胡励. 机械产品设计中的颗粒离散元仿真技术及实现[D]. 武汉：武汉大学，2013.

[7]M. Kremmer，J. F. Favier. A method for representing boundaries in discrete element modeling-part I: geometry andcontact detection [J]. Int. J. Numer. Meth. Eng.，2001（51）：1407-1421.

[8]P. W. Cleary. DEM prediction of industrial and geophysical particle flows [J]. Particuology，2010（8）：106-118.

[9]侯亚娟，王吉生，李爱峰，等. 离散元法及其在铲斗设计中的应用[J]. 工程机械，2014（45）：23-29.

[10]涂宇，王怡，刘爱云，等. 挖掘机工作装置机械-液压联合仿真系统的建立与优化[J]. 内燃机与配件，2019（23）：63-65.

[11]冯豪，杜群贵，于树栋. 液压挖掘机挖掘动力学建模研究[J].振动与冲击，2019（13）：242-248.

[12]阎书文. 机械式挖掘机设计[M]. 北京：机械工业出版社，1982：12-75.

[13]靳海军. 大型挖掘机斗齿定位监测系统研究[J]. 工程机械，2019（11）：95-98.

[14]M. B. Srđan，A. A. Miodrag，B. G. Nebojša，L. J. M. Ivan，M. A. Dušan. Failure of the bucket wheel excavator buckets [J]. Eng. Fail. Anal.，2018（84）：247-261.

[15]A. Dušan，G. Nebojša，S. Simon，A. Aleksandra，U. Milan. Integrity assessment and determination of residual fatigue life ofvital parts of bucket-wheel excavator operating under dynamicloads[J]. Eng. Fail. Anal.，2019（105）：182-195.

[16]C. J. Coetzee，D. N. J. Els. The numerical modelling of excavator bucket filling using DEM [J]. J. Terramechanics，2009（46）： 217-227.

[17]R. Eugeniusz，C. Lesław，M. Andrzej，M. Przemysław，O. Jerzy，P. Damian. Investigation and modernization of buckets of surface mining machines [J]. Eng. Struct.，2015（90）：29-37.

[18]Z. Q. Fang，G. M. Hu，J. Du，et al. A method for representing structural boundaries by mathematical equations and triangular meshes in DEM simulations [J]. Int. J. Numer. Meth. Eng.，2017，111（3）：218-246.

[19]胡国明，方自强，李婉婉，等. 多球颗粒离散元仿真计算软件[P]. 中国：2016SR154869，2016-06-24.

An Investigation on the Geometric Modeling for Structure Surface of Digging Bucket Using Elemental Spheres in the DEM Simulation

FANG Ziqiang

( School of Electromechanical and Architectural Engineering, Jianghan University, Wuhan 430056,China)

A multi-precision geometric modeling is proposed for the structure surface of digging bucket using elemental spheres. Based on the geometrical features and functionality requirements, this multi-precision geometric modeling is applied in analyzing the structure surfaces of the bottom plate, web, lateral plates, front plates and teeth in a digging bucket. Comparing the geometric models of the digging bucket with different precision, the result shows that the time consumed in the DEM simulation increases with the rise of the precision of geometric modeling for structure surface. Therefore, both the efficiency and veracity of DEM simulation should be taken into account in the geometric modeling. The comparison between our method and triangular facet method shows that the elemental spheres are able to express the local geometry characteristic of the digging bucket. It can be expected that our method takes less computational time than triangular facet method in the DEM simulation of digging process.

digging bucket；structure surface；geometric modeling；DEM simulation；elemental spheres

U294.8+8

A

10.3969/j.issn.1006-0316.2020.04.003

1006－0316 (2020) 04－0013－07

2019－12－16

湖北省自然科学基金（2019CFB208）；江汉大学高层次人才科研启动经费（1006－06640001）

方自强（1987－），男，湖北武汉人，博士，讲师，主要研究方向为颗粒系统离散元仿真技术研究与应用，fangziqiang1202@sina.com。