CN111104728B - 结构表面磨损的仿真预测方法和装置 - Google Patents

Abstract

本公开公开了一种结构表面磨损的仿真预测方法和装置，属于数据处理领域。所述仿真预测方法包括：形成连续排列的多个球形单元，所述多个球形单元的包络曲面用于表示结构表面；获取各个所述球形单元中所述结构表面的表示区域的法向量和面积；根据各个所述球形单元中所述结构表面的表示区域的法向量和面积，对各个球形颗粒与所述多个球形单元之间的碰撞进行离散元仿真，得到各个所述球形单元吸收的来自各个所述球形颗粒的碰撞能量；根据各个所述球形单元吸收的来自各个所述球形颗粒的碰撞能量，改变各个所述球形单元的球心位置，所述多个球形单元的球心位置的改变用于表示所述结构表面的磨损。本公开可有效减少仿真预测的计算时间。

Description

结构表面磨损的仿真预测方法和装置

技术领域

本公开涉及数据处理领域，特别涉及一种结构表面磨损的仿真预测方法和装置。

背景技术

在矿业、冶金、化工、食品加工等众多采用机械结构制备、处理散体(即固体颗粒的集合体)的工业过程中，颗粒撞击导致机械结构表面磨损是一种常见且复杂的物理现象。机械结构表面磨损预测对机械结构使用寿命的延长有着极为重要的意义。

离散元(英文：Distinct Element Method，简称：DEM)是专门用来解决不连续介质问题的数值模拟方法，可用于模拟颗粒与结构表面之间的碰撞行为，提供两者之间的碰撞能量、接触力、接触速度等信息。根据这些信息，可以有效预测受磨损结构表面在其整个生命周期中的几何演变过程。因此，DEM仿真已经广泛应用于众多机械结构表面的磨损预测。

在磨损预测中，结构表面会经历一系列的磨损演变。在每一步的磨损演变中，先采集DEM仿真的碰撞过程中颗粒作用于结构表面的磨损数据，再根据采集的磨损数据，更改结构表面的局部几何形状，响应颗粒的磨损作用。由于结构表面的几何形状一般采用许多小尺寸的边界单元表达，因此结构表面的磨损数据可以通过逐个搜集这些边界单元的接触情况来实现，结构表面的局部几何形状也可以通过逐个更新这些边界单元的位置来更改。

接触判断在DEM仿真中具有极其重要的作用。通过接触判断，各目标对象间的接触情况可以通过它们在之前一步磨损演变中的相对位置关系推导出来，这一过程对于当前一步磨损演变中的颗粒接触力计算是必须的。通常接触判断经历两个阶段：邻居检索和相交检测。邻居检索时将整个仿真区域规则地划分成相对较小的子区域，找出颗粒与边界单元之间的潜在接触对，从而将单个颗粒的接触检测范围缩小到离它最近的子区域，而非整个仿真区域。相交检测主要是利用目标对象的几何形体进行空间推导，对邻居检索中找到的潜在接触对象进行识别，排除无效的接触。相交检测的算法与颗粒的形状和结构表面的几何表达方式密切相关。

在DEM仿真中，最常使用到的颗粒形状是球体。而三角形单元是目前使用最为广泛的边界单元，已经应用在滚磨机的衬板、破碎机的破碎板、农业机械的触土部件等机械结构表面的几何形状的表达中。通过网格划分工具将机械结构表面划分成许多精度可控的三角形单元，将球形颗粒与任意形状的结构表面之间的接触检测转化为球形颗粒与三角形单元之间的接触检测；在每一步的磨损演变时，分别搜集各个三角形单元的磨损数据，并沿着磨损的深度方向逐个移动相应三角形单元的顶点，即可实现结构表面局部几何形状的更改。

在实现本公开的过程中，发明人发现现有技术至少存在以下问题：

在球形颗粒与三角形单元之间的相交检测中，球形颗粒需要与三角形单元的内部区域、三条边和三个顶点逐个进行相交检测。显然，这种相交检测的方式过于繁杂，影响接触检测时间和效率。特别是在DEM仿真由大量边界单元表达的机械结构表面时，接触检测会产生巨大的计算开销，消耗大量的计算机资源，每一步磨损演变对应的DEM仿真都可能会花费至少数周的计算时间，严重阻碍DEM仿真在磨损预测中的应用和发展。

发明内容

本公开实施例提供了一种结构表面磨损的仿真预测方法和装置，能够解决现有技术检测时间长、检测效率低的问题。所述技术方案如下：

一方面，本公开实施例提供了一种结构表面磨损的仿真预测方法，所述仿真预测方法包括：

形成连续排列的多个球形单元，所述多个球形单元的包络曲面用于表示结构表面；

获取各个所述球形单元中所述结构表面的表示区域的法向量和面积；

根据各个所述球形单元中所述结构表面的表示区域的法向量和面积，对各个球形颗粒与所述多个球形单元之间的碰撞进行离散元仿真，得到各个所述球形单元吸收的来自各个所述球形颗粒的碰撞能量；

根据各个所述球形单元吸收的来自各个所述球形颗粒的碰撞能量，改变各个所述球形单元的球心位置，所述多个球形单元的球心位置的改变用于表示所述结构表面的磨损。

可选地，所述形成连续排列的多个球形单元，所述多个球形单元的包络曲面用于表示结构表面，包括：

绘制所述结构表面；

将所述结构表面沿所述结构表面的内法线方向偏移设定距离，生成内偏移曲面；

将所述内偏移曲面划分为多个三角形单元；

以各个所述三角形单元的顶点为球心，设定距离为半径，形成所述球形单元。

进一步地，所述获取各个所述球形单元中所述结构表面的表示区域的法向量和面积，包括：

根据以所述球形单元的球心为顶点的所有三角形单元的法向量之和，得到所述球形单元中所述结构表面的表示区域的法向量；

根据以所述球形单元的球心为顶点的所有三角形单元中所述结构表面的表示区域的面积之和，得到所述球形单元中所述结构表面的表示区域的面积。

进一步地，所述根据各个所述球形单元中所述结构表面的表示区域的法向量和面积，对各个球形颗粒与所述多个球形单元之间的碰撞进行离散元仿真，得到各个所述球形单元吸收的来自各个所述球形颗粒的碰撞能量，包括：

在离散元仿真中每隔设定时长判定各个所述球形颗粒与各个所述球形单元是否接触；

当本次判定所述球形颗粒与所述球形单元发生接触，且在上一次判定中与所述球形颗粒发生接触的球形单元的数量等于0时，确定所述多个球形单元吸收的来自所述球形颗粒的碰撞能量；

当本次判定所述球形颗粒与所述球形单元发生接触时，记录与所述球形颗粒发生接触的球形单元；

当上一次判定所述球形颗粒与所述球形单元发生接触，且在本次判定中与所述球形颗粒发生接触的球形单元的数量等于0时，将所述多个球形单元吸收的来自所述球形颗粒的碰撞能量分配给与所述球形颗粒发生接触的球形单元。

更进一步地，所述根据各个所述球形单元吸收的来自各个所述球形颗粒的碰撞能量，改变各个所述球形单元的球心位置，所述多个球形单元的球心位置的改变用于表示所述结构表面的磨损，包括：

根据各个所述球形单元吸收的来自各个所述球形颗粒的碰撞能量，预测各个所述球形单元的磨损体积；

根据各个所述球形单元的磨损体积，确定各个所述球形单元的球心位置；

按照确定的球心位置移动所述球形单元。

可选地，所述仿真预测方法还包括：

当相邻两个所述球形单元之间的缝隙长度达到阈值时，在相邻两个所述球形单元之间增设一个所述球形单元，增设的所述球形单元的球心与相邻两个所述球形单元的球心连线的中点重合。

进一步地，所述仿真预测方法还包括：

确定所述阈值。

可选地，所述仿真预测方法还包括：

在所述球形单元的球面上选择采样点；

将所述多个球形单元的采样点进行拟合，形成多个球形单元的包络曲面。

另一方面，本公开实施例提供了一种结构表面磨损的仿真预测装置，所述仿真预测装置包括：

形成模块，用于形成连续排列的多个球形单元，所述多个球形单元的包络曲面用于表示结构表面；

获取模块，用于获取各个所述球形单元中所述结构表面的表示区域的法向量和面积；

仿真模块，用于根据各个所述球形单元中所述结构表面的表示区域的法向量和面积，对各个球形颗粒与所述多个球形单元之间的碰撞进行离散元仿真，得到各个所述球形单元吸收的来自各个所述球形颗粒的碰撞能量；

改变模块，用于根据各个所述球形单元吸收的来自各个所述球形颗粒的碰撞能量，改变各个所述球形单元的球心位置，所述多个球形单元的球心位置的改变用于表示所述结构表面的磨损。

可选地，所述仿真预测装置还包括：

采样模块，用于在所述球形单元的球面上选择采样点；

表示模块，用于将所述多个球形单元的采样点进行拟合，形成多个球形单元的包络曲面。

本公开实施例提供的技术方案带来的有益效果是：

通过采用多个球形单元的包络曲面表示结构表面，使得球形颗粒与结构表面之间的相交检测仅涉及一个球心距离与两个半径和之间的比较，相交检测方式得到有效简化，可以大幅减少接触检测产生的计算开销，降低接触检测的时间消耗，有利于DEM仿真在磨损预测中的应用和发展。

附图说明

为了更清楚地说明本公开实施例中的技术方案，下面将对实施例描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本公开的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。

图1是本公开实施例提供的一种结构表面磨损的仿真预测方法的流程图；

图2是本公开实施例提供的内偏移曲面形成过程的示意图；

图3是本公开实施例提供的内偏移曲面划分三角形单元的示意图；

图4是本公开实施例提供的球形单元表达结构表面的示意图；

图5是本公开实施例提供的球形单元的相邻内三角形的示意图；

图6是本公开实施例提供的相邻球形单元的示意图；

图7是本公开实施例提供的球形单元中结构表面的表示区域的法向量的示意图；

图8是本公开实施例提供的三角形单元划分子三角形的示意图；

图9是本公开实施例提供的三角形单元划分子三角形的示意图；

图10是本公开实施例提供的立方体网格的示意图；

图11是本公开实施例提供的相交检测条件的示意图；

图12是本公开实施例提供的球形颗粒与球形单元碰撞的示意图；

图13是本公开实施例提供的球形单元的球心位置变化之后的示意图；

图14是本公开实施例提供的识别相邻球形单元之间缝隙的示意图；

图15是本公开实施例提供的处理相邻球形单元之间缝隙的示意图；

图16是本公开实施例提供的采样点拟合形成包络曲面的示意图；

图17是本公开实施例提供的工业滚磨机上曲面提升条表面各个区域经历磨损的变化情况图；

图18是本公开实施例提供的曲面提升条表面采用不同球心距离形成的球形单元的示意图；

图19是本公开实施例提供的不同球心距离下滚磨机中颗粒运动形态的DEM截图；

图20是本公开实施例提供的不同球心距离下预测的提升条表面的三维几何形状的示意图；

图21是本公开实施例提供的不同球心距离下预测的提升条表面的二维平均轮廓及准确率的示意图；

图22是本公开实施例提供的不同球心距离下本预测方法与不同三角网格尺寸下现有预测方法的接触检测时间的对比图；

图23是本公开实施例提供的一种结构表面磨损的仿真预测的结构示意图。

具体实施方式

为使本公开的目的、技术方案和优点更加清楚，下面将结合附图对本公开实施方式作进一步地详细描述。

本公开实施例提供了一种结构表面磨损的仿真预测方法。图1为本公开实施例提供的一种结构表面磨损的仿真预测方法的流程图。参见图1，仿真预测方法包括：

步骤101：形成连续排列的多个球形单元，多个球形单元的包络曲面用于表示结构表面。

在本实施例中，球形单元是表达结构表面的最小边界单元；表达结构表面的球形单元的数量、尺寸及其组合形式，决定了所表达的表达结构表面的几何形状精度。

相邻两个球形单元之间可以相互重叠或者存在缝隙。各个球形单元的尺寸均由相同大小的直径d_i定义，直径d_i可以小于或等于球形颗粒的直径D_m的三分之一。

可选地，该步骤101可以包括：

第一步，绘制结构表面；

第二步，将结构表面沿结构表面的内法线方向偏移设定距离，生成内偏移曲面；

第三步，将内偏移曲面划分为多个三角形单元；

第四步，以各个三角形单元的顶点为球心，设定距离为半径，形成球形单元。

在实际应用中，在第一步中，可以以SolidWorks或ProE等三维CAD(英文：ComputerAided Design，中文：计算机辅助设计)软件为平台，设计和绘制结构表面Σ_O。

在第二步中，可以利用上述三维CAD平台，沿着结构表面Σ_O的内法线方向，向内整体偏移结构表面Σ_O，从而生成一个内偏移曲面Σ_I。其中，偏移距离等于球形单元的半径d_i/2。图2为本公开实施例提供的内偏移曲面形成过程的示意图。如图2所示，将结构表面Σ_O中除边缘之外的所有区域，沿着结构表面Σ_O的外法线方向n_i的反方向，向内整体偏移d_i/2，得到内偏移曲面Σ_I。例如，结构表面Σ_O上的中心Q_i沿着结构表面Σ_O上在中心Q_i处的外法线方向n_i的反方向偏移d_i/2，变成内偏移曲面Σ_I上的中心S_i。

在第三步中，可以利用ICEM CFD或Gambit等网格划分工具，将内偏移曲面Σ_I分解为多个三角形单元，相邻两个三角形单元的一条边重合。其中，三角形单元的最大边长△_max小于或等于球形单元的半径d_i/2，即三角形单元的每条边长均小于或等于球形单元的半径d_i/2，以确保形成的各个球形单元之间相互重叠，不存在缝隙。图3为本公开实施例提供的内偏移曲面划分三角形单元的示意图。参见图3，内偏移曲面Σ_I由多个三角形单元k组成。

在第四步中，可以将直径同为d_i的球形单元的球心逐一添加至三角形单元的各个顶点位置上，所有球形单元整体的包络曲面表达结构表面Σ_O的物理边界。图4为本公开实施例提供的球形单元表达结构表面的示意图。参见图4，以内偏移曲面Σ_I上的三角形单元的顶点S_i为球心、d_i/2为半径形成球形单元i，多个球形单元i整体的包络曲面表示结构表面Σ_O。

步骤102：获取各个球形单元中结构表面的表示区域的法向量和面积。

可选地，该步骤102可以包括：

根据以球形单元的球心为顶点的所有三角形单元的法向量之和，得到球形单元中结构表面的表示区域的法向量；

根据以球形单元的球心为顶点的所有三角形单元中结构表面的表示区域的面积之和，得到球形单元中结构表面的表示区域的面积。

在实际应用中，当一个三角形单元的三个顶点中的任意一个顶点与一个球形单元的球心重合时，这个三角形单元为这个球形单元的相邻内三角形。图5为本公开实施例提供的球形单元的相邻内三角形的示意图。参见图5，以三角形单元的顶点S_i为球心的球形单元的相邻内三角形包括：以S_i、S_j-2和S_j-1为顶点的三角形单元(k-2)，以S_i、S_j-1和S_j为顶点的三角形单元(k-1)，以S_i、S_j和S_j+1为顶点的三角形单元k，以S_i、S_j+1和S_j+2为顶点的三角形单元(k+1)，以S_i、S_j+2和S_j+3为顶点的三角形单元(k+2)，以S_i、S_j+2和S_j+3为顶点的三角形单元(k+3)。

当两个球形单元的球心分别与同一个三角形单元的不同顶点重合时，这两个球形单元相邻。图6为本公开实施例提供的相邻球形单元的示意图。参见图6，以S_i为球心的球形单元i的相邻内三角形包括：以S_j-2为球心的球形单元(j-2)，以S_j-1为球心的球形单元(j-1)，以S_j为球心的球形单元j，以S_j+1为球心的球形单元(j+1)，以S_j+2为球心的球形单元(j+2)，以S_j+3为球心的球形单元(j+3)。

遍历内偏移曲面Σ_I上所有三角形单元的顶点，可以得到各个球形单元i的相邻内三角形；遍历结构表面Σ_O上所有球形单元的球心，可以得到各个球形单元i的相邻球形单元。

球形单元中结构表面的表示区域的法向量，为结构表面Σ_O在球形单元i处的合法向量(n_i)^o，由这个球形单元i的所有相邻内三角形的法向量的线性加权和得到。具体可以采用如下公式(1)计算：

其中，K是球形单元i的相邻内三角形的数量，n_k是三角形单元k的外法线单位向量。图7为本公开实施例提供的球形单元中结构表面的表示区域的法向量的示意图。参见图7，结构表面Σ_O在球形单元i处的合法向量(n_i)^o，由以三角形单元的顶点S_i为球心的球形单元i的所有相邻内三角形的法向量的线性加权和得到，即由以S_i、S_j-2和S_j-1为顶点的三角形单元(k-2)的法向量n_k-2，以S_i、S_j-1和S_j为顶点的三角形单元(k-1)的法向量n_k-1，以S_i、S_j和S_j+1为顶点的三角形单元k的法向量n_k，以S_i、S_j+1和S_j+2为顶点的三角形单元(k+1)的法向量n_k+1，以S_i、S_j+2和S_j+3为顶点的三角形单元(k+2)的法向量n_k+2，以S_i、S_j+2和S_j+3为顶点的三角形单元(k+3)的法向量n_k+3加权求和得到。

球形单元中结构表面的表示区域的面积Ω_i，可以采用如下公式(2)计算：

其中Ω_i,k,j和Ω_i,k,(j+1)分别为球形单元i的相邻内三角形k中子三角形△i,k,j和△i,k,(j+1)的面积；ψ_i,k为球形单元i中结构表面Σ_O的表示区域在三角形单元k内的面积与Ω_i,k,j和Ω_i,k,(j+1)之和的比值；子三角形△i,k,j和△i,k,(j+1)为三角形单元k的六个子三角形△i,k,j、△i,k,(j+1)、△j,k,i、△j,k,(j+1)、△(j+1),k,i和△(j+1),k,j中，位于球形单元i内的两个；上述六个子三角形由三角形单元k分解而来，分解点包括三角形单元k的三个顶点S_i、S_j和S_(j+1)，内心C_k，内心C_k在三角形单元k的三条边S_iS_j、S_iS_(j+1)和S_jS_(j+1)上的投影点C_i,j、C_i,(j+1)和C_j,(j+1)。

图8和图9为本公开实施例提供的三角形单元划分子三角形的示意图。参见图8和图9，三角形单元k按照三个顶点S_i、S_j和S_(j+1)，内心C_k，内心C_k在三角形单元k的三条边S_iS_j、S_iS_(j+1)和S_jS_(j+1)上的投影点C_i,j、C_i,(j+1)和C_j,(j+1)划分成六个子三角形△i,k,j、△i,k,(j+1)、△j,k,i、△j,k,(j+1)、△(j+1),k,i和△(j+1),k,j，阴影表示的两个子三角形△i,k,j和△i,k,(j+1)为六个子三角形中位于球形单元i内的子三角形。图8和图9的不同之处在于，图8中相邻球形单元之间没有缝隙，图9中相邻球形单元之间存在缝隙δ_ij。

步骤103：根据各个球形单元中结构表面的表示区域的法向量和面积，对各个球形颗粒与多个球形单元之间的碰撞进行离散元仿真，得到各个球形单元吸收的来自各个球形颗粒的碰撞能量。

可选地，该步骤103可以包括：

在离散元仿真中每隔设定时长判定各个球形颗粒与各个球形单元是否接触；

当本次判定球形颗粒与球形单元发生接触，且在上一次判定中与球形颗粒发生接触的球形单元的数量等于0时，确定多个球形单元吸收的来自球形颗粒的碰撞能量；

当本次判定球形颗粒与球形单元发生接触时，记录与球形颗粒发生接触的球形单元；

当上一次判定球形颗粒与球形单元发生接触，且在本次判定中与球形颗粒发生接触的球形单元的数量等于0时，将多个球形单元吸收的来自球形颗粒的碰撞能量分配给与球形颗粒发生接触的球形单元。

在实际应用中，可以先进行网络搜索，如采用基于均匀立方体网格划分的空间邻居搜索方法，即将整个DEM仿真区域规则地划分成立方体网格，从而将球形颗粒的接触检测范围缩小至离它最近的多个立方体网格。图10为本公开实施例提供的立方体网格的示意图。参见图10，整个DEM仿真区域规则地划分成多个立方体网格。其中，立方体网格的尺寸L_c由球形颗粒的直径D_m决定。

图11为本公开实施例提供的相交检测条件的示意图。参见图11，球形颗粒m与球形单元i之间的相交检测条件，可以如下式子(3)所示：

其中，|S_iP_m|为球形颗粒m与球形单元i的球心距离；P_m为球形颗粒m的球心，S_i为球形单元i的球心。

结构表面Σ_O吸收的来自球形颗粒m冲击磨损的碰撞能量

可以采用如下公式(4)计算：

其中，v_Threshold为临界磨损速度(通常为0.1m/s)；M_m为球形颗粒m的质量；v_Σ,m为颗粒m与结构表面Σ_O在球形单元i处的相对速度，它统计于颗粒m与球形单元i刚接触时，即本次判定球形颗粒与球形单元发生接触，且在上一次判定中与球形颗粒发生接触的球形单元的数量等于0；相对速度v_Σ,m采用如下公式计算：v_Σm＝v_m-v_Σ，v_Σ为结构表面Σ_O在球形单元i处的运动速度；v_m为球形颗粒m的运动速度，它们均统计于颗粒m与球形单元i刚接触时。

结构表面Σ_O吸收的来自球形颗粒m磨粒磨损的碰撞能量

可以采用如下公式(5)计算：

其中，α_max为结构表面Σ_O产生最大磨粒磨损的颗粒入射角(一般为22°)；M_m为球形颗粒m的质量；v_Σ,m为颗粒m与结构表面Σ_O在球形单元i处的相对速度，它统计于颗粒m与球形单元i刚接触时，即本次判定球形颗粒与球形单元发生接触，且在上一次判定中与球形颗粒发生接触的球形单元的数量等于0；相对速度vΣ,m采用如下公式计算：v_Σm＝v_m-v_Σ，v_Σ为结构表面Σ_O在球形单元i处的运动速度；v_m为球形颗粒m的运动速度，它们均统计于颗粒m与球形单元i刚接触时，如图12所示；α_Σ,m为球形颗粒m与结构表面Σ_O在球形单元i处的相对入射角，也统计于颗粒m和球形单元i刚接触时，相对入射角α_Σ,m采用如下公式计算：

(n_i)^o为结构表面Σ_O在球形单元i处的合法向量。

将多个球形单元吸收的来自球形颗粒的碰撞能量分配给与球形颗粒发生接触的球形单元，即将结构表面Σ_O吸收的来自球形颗粒m冲击磨损的碰撞能量

和磨粒磨损的碰撞能量

分配给球形单元i、以及与球形单元i在同一次碰撞中、与球形颗粒m发生接触的其它球形单元，可以采用这些球形单元中结构表面Σ_O的表示区域的面积比值为依据进行量化计算。具体地，球形单元i吸收的来自球形颗粒m冲击磨损的碰撞能量

和磨粒磨损的碰撞能量

可以分别采用如下公式(6)和(7)计算：

其中，Ω_i为球形单元i中结构表面的表示区域的面积；g为在同一次碰撞中，与同一个球形颗粒m接触的球形单元的序列号码；G为在同一次碰撞中，与同一个球形颗粒m接触的球形单元的数量；Ω_g为球形单元g中结构表面Σ_O的表示区域的面积。

由于碰撞是由球形颗粒产生的，因此在执行步骤103时，需要遍历所有的球形颗粒m作用在结构表面上的碰撞能量。另外，离散元仿真将仿真时间每隔设定时长分为一个阶段，每个阶段的持续时长即为前述每次判定球形颗粒与球形单元是否接触相隔的设定时长。因此在执行步骤103时，还需要遍历所有的仿真阶段，并在每个仿真阶段遍历所有的球形颗粒m作用在结构表面上的碰撞能量。

步骤104：根据各个球形单元吸收的来自各个球形颗粒的碰撞能量，改变各个球形单元的球心位置，多个球形单元的球心位置的改变用于表示结构表面的磨损。

可选地，该步骤104可以包括：

根据各个球形单元吸收的来自各个球形颗粒的碰撞能量，预测各个球形单元的磨损体积；

根据各个球形单元的磨损体积，确定各个球形单元的球心位置；

按照确定的球心位置移动球形单元。

在实际应用中，结构表面Σ_O在球形单元i处的磨损体积V_i ^e，可以采用如下公式(8)计算：

其中，W_i是磨损率，单位是[mm³/J]；aⁱ _Impact是碰撞过程中冲击磨损作用的影响系数，aⁱ _Abrasion颗粒碰撞过程中磨粒磨损作用的影响系数；c为单个DEM仿真数据搜集周期τd内球形单元i与球形颗粒m之间接触的序列号码；N_i,m为该数据搜集周期τd内球形单元i与球形颗粒m之间接触的总次数；N_i为在该数据搜集周期τd内，与球形单元i发生接触的球形颗粒m的总数量；τe为与该数据搜集周期τd对应的磨损演变周期时长。

结构表面Σ_O在球形单元i处几何形状的局部变化是通过移动球形单元i的球心来完成；球形单元i的球心移动后的位置S_i ^ξ，可以采用如下公式(9)计算：

其中，S_i为球形单元i的球心移动前的位置，(n_i)^o为结构表面Σ_O在球形单元i处的合法向量，V_i ^e为结构表面Σ_O在球形单元i处的磨损体积，Ω_i为球形单元i中结构表面的表示区域的面积。

图13为本公开实施例提供的球形单元的球心位置变化之后的示意图。参见图13，球形单元i的球心S_i从沿结构表面Σ_O在球形单元i处的合法向量(n_i)^o的反方向移动到S_i ^ξ。

本公开实施例通过采用多个球形单元的包络曲面表示结构表面，使得球形颗粒与结构表面之间的相交检测仅涉及一个球心距离与两个半径和之间的比较，相交检测方式得到有效简化，可以大幅减少接触检测产生的计算开销，降低接触检测的时间消耗，有利于DEM仿真在磨损预测中的应用和发展。

在实际应用中，结构表面在磨损预测中会经历一系列的磨损演变，每一步的磨损演变都需要执行一次步骤102～步骤104。因此，整个仿真预测的过程中，会循环执行步骤102～步骤104，直到完成所有步的磨损演变。

可选地，仿真预测方法还可以包括：

当相邻两个球形单元之间的缝隙长度达到阈值时，在相邻两个球形单元之间增设一个球形单元，增设的球形单元的球心与相邻两个球形单元的球心连线的中点重合。

图14为本公开实施例提供的识别相邻球形单元之间缝隙的示意图。参见图14，在实际应用中，当球形单元i与其两个相邻球形单元j和(j+1)之间的球面间距δ_ij、δ_i(j+1)、δ_j(j+1)满足如下式子(10)时，确定相邻球形单元之间的缝隙长度达到阈值δ_max：

(δ_ij≥δ_max)∪(δ_i(j+1)≥δ_max)∪(δ_j(j+1)≥δ_max)；  (10)

其中，球面间距δ_ij、δ_i(j+1)、δ_j(j+1)可以采用如下公式(11)计算：

δ_ij＝|S_iS_j|-(d_i+d_j)/2,

δ_i(j+1)＝|S_iS_(j+1)|-(d_i+d_(j+1))/2,；  (11)

δ_j(j+1)＝|S_jS_(j+1)|-(d_j+d_(j+1))/2

|S_iS_j|为球形单元i与球形单元j的球心距离，|S_iS_(j+1)|为球形单元i与球形单元(j+1)的球心距离，|S_jS_(j+1)|为球形单元j与球形单元(j+1)的球心距离；d_i为球形单元i的直径，d_j为球形单元j的直径，d_(j+1)为球形单元(j+1)的直径。

图15为本公开实施例提供的处理相邻球形单元之间缝隙的示意图。参见图15，当δ_ij＞δ_max、δ_i(j+1)＜δ_max、δ_j(j+1)＜δ_max时，在球形单元i和球形单元j之间增加一个球形单元New，球形单元New与球形单元i之间的球面间距δ_iNew小于阈值δ_max，球形单元New与球形单元j之间的球面间距δ_jNew小于阈值δ_max，即。

进一步地，仿真预测方法还可以包括：

确定阈值。

在实际应用中，阈值δ_max可以采用如下公式(12)计算：

D_m为球形颗粒m的直径，d_i为球形单元i的直径。

可选地，仿真预测方法还可以包括：

在球形单元的球面上选择采样点；

将多个球形单元的采样点进行拟合，形成多个球形单元的包络曲面。

在实际应用中，采用多个球形单元的包络曲面描述结构表面Σ_O的外观几何形貌，包络曲面是一个基于结构表面Σ_O中球形单元i球面上采样点S_i ^E的拟合曲面，可以通过调用MATLAB工具箱中的Griddata方程来生成。具体地，球形单元i球面上的采样点S_i ^E，可以采用如下公式(13)计算：

其中，S_i ^ξ为球形单元i的球心移动后的位置，(n_i)^o为结构表面Σ_O在球形单元i处的合法向量，d_i为球形单元i的直径。

图16为本公开实施例提供的采样点拟合形成包络曲面的示意图。参见图16，将球形单元i上的取样点S_i ^E，球形单元(j-2)上的取样点S_(j-2) ^E，球形单元(j-1)上的取样点S_(j-1) ^E，球形单元j上的取样点S_j ^E，球形单元(j+1)上的取样点S_(j+1) ^E，球形单元(j+2)上的取样点S_(j+2) ^E，球形单元(j+3)上的取样点S_(j+3) ^E等拟合形成包络面Σ_O ^E。

将本公开实施例提供的仿真预测方法写成VC++计算机仿真程序，应用在一个工业滚磨机中颗粒研磨运动上进行检验。使用的计算机处理器为Intel Core i7-6700K CPU，主频为4.00GHz，RAM(英文：Random Access Memory，中文：随机存取存储器)为16GB，编程工具为Visual Studio 2010。

图17为本公开实施提供的工业滚磨机上曲面提升条表面各个区域经历磨损的变化情况图。参见图17，工业滚磨机包括弯成圆环形的曲面提升条100和分别设置曲面提升条100两侧的圆形端面平板200。截取曲面提升条100的部分区域进行物理测量，发现各个区域在表面磨损之前的厚度相同，在经历57726h磨损之后轮廓的变化情况不同。根据轮廓变化的大小分为磨损区域1、磨损区域2和磨损区域3。参照实际工作条件，磨机的转速和填充率分别为15.7rpm和24.2％(也就是1600个直径D_m为φ50mm的球形钢质颗粒)，提升条表面的真实生命时长为57726h。

图18为本公开实施例提供的曲面提升条表面采用不同球心距离形成的球形单元的示意图。参见图18，曲面提升条表面的凸起部和凹陷部沿周期性交替出现，截取一个凸起部和两侧的凹陷部分别采用不同球形距离的多个球形单元进行表示。例如，球心距离ζ＝d_i/2＝(D_m/4)/2＝6.25mm；又如，球心距离ζ＝3*d_i/4＝3*(D_m/4)/4＝9.38mm；再如，球心距离ζ＝d_i＝D_m/4＝12.5mm。其中，d_i为球形单元的直径，D_m为球形颗粒的直径。

将磨损演变的步数Esp＝150，数据搜集周期τ_d＝11.47s，在DEM仿真中，每个阶段的持续时长为1×10-5s。基于不同球形距离的多个球形单元分别预测工业滚磨机中受磨损提升条表面的磨损几何演变结果；并在三种不同ζ的条件下，比较本公开提出的预测方法所预测的提升条表面二维平均轮廓的准确率、比较本公开提出的预测方法与现有基于三角形单元的预测方法在接触检测中的时间消耗。

图19为本公开实施例提供的不同球心距离下滚磨机中颗粒运动形态的DEM截图。参见图19，三种不同球心距离ζ对应的提升条表面的颗粒运动形态基本一致。

图20为本公开实施例提供的不同球心距离下预测的提升条表面的三维几何形状的示意图。参见图20，从外观形貌上看，三种不同球心距离ζ对应的提升条表面的三维几何形状基本一致。

图21为本公开实施例提供的不同球心距离下预测的提升条表面的二维平均轮廓及准确率的示意图。参见图21，经过150步的磨损演变，三种不同ζ对应的提升条表面的二维平均轮廓与物理测量获取的二维平均轮廓均基本保持一致。而且三种ζ对应的二维轮廓的平均准确率中的绝大多数在90％以上，有效证实了本方法在磨损预测中的准确性和稳定性。

当ζ为d_i/2、d_i*(3/4)和d_i时，本预测方法中用于表达提升条表面的球形单元的数量分别为26160、13800和10560。相应地，当基于三角形单元的预测方法中三角网格的最大尺寸△_max为11.6、14.9和16.6mm时，用于表达提升条表面的三角网格的数量分别为25800、13440和10320。

图22为本公开实施例提供的不同球心距离下本预测方法与不同三角网格尺寸下现有预测方法的接触检测时间的对比图。其中，现有预测方法中三角形单元的接触检测时间T_c，为磨机每旋转一圈所消耗的DEM仿真中颗粒与提升条表面间的接触检测的计算机程序运行时间。参见图22，当ζ和△_max分别减小时，本预测方法和现有预测方法的T_c均随之分别增加。对于不同的ζ和△_max，本预测方法的T_c均远小于现有预测方法的T_c；当ζ或△_max较小时，本预测方法的T_c的大小仅等于现有预测方法的T_c大小的1.5％。

上述结果表明，相较于现有基于三角形单元的预测方法，本预测方法极大地提升了DEM仿真中颗粒与受磨损结构表面之间接触检测的时间效率，解决了接触检测在结构表面磨损的DEM仿真中计算开销过大的问题；此外，本预测方法还准确且稳定地预测了不同几何表达分辨率下工业滚磨机中受磨损提升条表面的几何演变结果，为DEM仿真在结构表面磨损预测中的大规模应用提供了有效技术手段。

本公开实施例提供了一种结构表面磨损的仿真预测装置，适用于实现图1所示的一种结构表面磨损的仿真方法。图23为本公开实施例提供的一种结构表面磨损的仿真预测的结构示意图。参见图23，仿真预测装置包括：

形成模块201，用于形成连续排列的多个球形单元，多个球形单元的包络曲面用于表示结构表面；

获取模块202，用于获取各个球形单元中结构表面的表示区域的法向量和面积；

仿真模块203，用于根据各个球形单元中结构表面的表示区域的法向量和面积，对各个球形颗粒与多个球形单元之间的碰撞进行离散元仿真，得到各个球形单元吸收的来自各个球形颗粒的碰撞能量；

改变模块204，用于根据各个球形单元吸收的来自各个球形颗粒的碰撞能量，改变各个球形单元的球心位置，多个球形单元的球心位置的改变用于表示结构表面的磨损。

可选地，形成模块201可以包括：

绘制单元，用于绘制结构表面；

生长单元，用于将结构表面沿结构表面的内法线方向偏移设定距离，生成内偏移曲面；

划分单元，用于将内偏移曲面划分为多个三角形单元；

形成单元，用于以各个三角形单元的顶点为球心，设定距离为半径，形成球形单元。

进一步地，获取模块202可以用于，

根据以球形单元的球心为顶点的所有三角形单元的法向量之和，得到球形单元中结构表面的表示区域的法向量；

根据以球形单元的球心为顶点的所有三角形单元中结构表面的表示区域的面积之和，得到球形单元中结构表面的表示区域的面积。

进一步地，仿真模块203可以用于，

在离散元仿真中每隔设定时长判定各个球形颗粒与各个球形单元是否接触；

当本次判定球形颗粒与球形单元发生接触，且在上一次判定中与球形颗粒发生接触的球形单元的数量等于0时，确定多个球形单元吸收的来自球形颗粒的碰撞能量；

当本次判定球形颗粒与球形单元发生接触时，记录与球形颗粒发生接触的球形单元；

当上一次判定球形颗粒与球形单元发生接触，且在本次判定中与球形颗粒发生接触的球形单元的数量等于0时，将多个球形单元吸收的来自球形颗粒的碰撞能量分配给与球形颗粒发生接触的球形单元。

更进一步地，改变模块204可以包括：

预测单元，用于根据各个球形单元吸收的来自各个球形颗粒的碰撞能量，预测各个球形单元的磨损体积；

确定单元，用于根据各个球形单元的磨损体积，确定各个球形单元的球心位置；

移动单元，用于按照确定的球心位置移动球形单元。

可选地，仿真预测装置还可以包括：

增设模块，用于当相邻两个球形单元之间的缝隙长度达到阈值时，在相邻两个球形单元之间增设一个球形单元，增设的球形单元的球心与相邻两个球形单元的球心连线的中点重合。

进一步地，仿真预测装置还可以包括：

确定模块，用于确定阈值。

可选地，仿真预测装置还可以包括：

采样模块，用于在球形单元的球面上选择采样点；

表示模块，用于将多个球形单元的采样点进行拟合，形成多个球形单元的包络曲面。

需要说明的是：上述实施例提供的结构表面磨损的仿真预测装置在仿真预测结构表面磨损时，仅以上述各功能模块的划分进行举例说明，实际应用中，可以根据需要而将上述功能分配由不同的功能模块完成，即将装置的内部结构划分成不同的功能模块，以完成以上描述的全部或者部分功能。另外，上述实施例提供的结构表面磨损的仿真预测装置与结构表面磨损的仿真预测方法实施例属于同一构思，其具体实现过程详见方法实施例，这里不再赘述。

上述本公开实施例序号仅仅为了描述，不代表实施例的优劣。

本领域普通技术人员可以理解实现上述实施例的全部或部分步骤可以通过硬件来完成，也可以通过程序来指令相关的硬件完成，所述的程序可以存储于一种计算机可读存储介质中，上述提到的存储介质可以是只读存储器，磁盘或光盘等。

以上所述仅为本公开的较佳实施例，并不用以限制本公开，凡在本公开的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本公开的保护范围之内。

Claims (10)

1.一种结构表面磨损的仿真预测方法，其特征在于，所述仿真预测方法包括：

形成连续排列的多个球形单元，所述多个球形单元的包络曲面用于表示结构表面；

获取各个所述球形单元中所述结构表面的表示区域的法向量和面积；

根据各个所述球形单元中所述结构表面的表示区域的法向量和面积，对各个球形颗粒与所述多个球形单元之间的碰撞进行离散元仿真，得到各个所述球形单元吸收的来自各个所述球形颗粒的碰撞能量；

根据各个所述球形单元吸收的来自各个所述球形颗粒的碰撞能量，改变各个所述球形单元的球心位置，所述多个球形单元的球心位置的改变用于表示所述结构表面的磨损。

2.根据权利要求1所述的仿真预测方法，其特征在于，所述形成连续排列的多个球形单元，所述多个球形单元的包络曲面用于表示结构表面，包括：

绘制所述结构表面；

将所述结构表面沿所述结构表面的内法线方向偏移设定距离，生成内偏移曲面；

将所述内偏移曲面划分为多个三角形单元；

以各个所述三角形单元的顶点为球心，设定距离为半径，形成所述球形单元。

3.根据权利要求2所述的仿真预测方法，其特征在于，所述获取各个所述球形单元中所述结构表面的表示区域的法向量和面积，包括：

根据以所述球形单元的球心为顶点的所有三角形单元的法向量之和，得到所述球形单元中所述结构表面的表示区域的法向量；

根据以所述球形单元的球心为顶点的所有三角形单元中所述结构表面的表示区域的面积之和，得到所述球形单元中所述结构表面的表示区域的面积。

4.根据权利要求2或3所述的仿真预测方法，其特征在于，所述根据各个所述球形单元中所述结构表面的表示区域的法向量和面积，对各个球形颗粒与所述多个球形单元之间的碰撞进行离散元仿真，得到各个所述球形单元吸收的来自各个所述球形颗粒的碰撞能量，包括：

在离散元仿真中每隔设定时长判定各个所述球形颗粒与各个所述球形单元是否接触；

当本次判定所述球形颗粒与所述球形单元发生接触，且在上一次判定中与所述球形颗粒发生接触的球形单元的数量等于0时，确定所述多个球形单元吸收的来自所述球形颗粒的碰撞能量；

当本次判定所述球形颗粒与所述球形单元发生接触时，记录与所述球形颗粒发生接触的球形单元；

当上一次判定所述球形颗粒与所述球形单元发生接触，且在本次判定中与所述球形颗粒发生接触的球形单元的数量等于0时，将所述多个球形单元吸收的来自所述球形颗粒的碰撞能量分配给与所述球形颗粒发生接触的球形单元。

5.根据权利要求4所述的仿真预测方法，其特征在于，所述根据各个所述球形单元吸收的来自各个所述球形颗粒的碰撞能量，改变各个所述球形单元的球心位置，所述多个球形单元的球心位置的改变用于表示所述结构表面的磨损，包括：

根据各个所述球形单元吸收的来自各个所述球形颗粒的碰撞能量，预测各个所述球形单元的磨损体积；

根据各个所述球形单元的磨损体积，确定各个所述球形单元的球心位置；

按照确定的球心位置移动所述球形单元。

6.根据权利要求1～3任一项所述的仿真预测方法，其特征在于，所述仿真预测方法还包括：

当相邻两个所述球形单元之间的缝隙长度达到阈值时，在相邻两个所述球形单元之间增设一个所述球形单元，增设的所述球形单元的球心与相邻两个所述球形单元的球心连线的中点重合。

7.根据权利要求6所述的仿真预测方法，其特征在于，所述仿真预测方法还包括：

确定所述阈值。

8.根据权利要求1～3任一项所述的仿真预测方法，其特征在于，所述仿真预测方法还包括：

在所述球形单元的球面上选择采样点；

将所述多个球形单元的采样点进行拟合，形成多个球形单元的包络曲面。

9.一种结构表面磨损的仿真预测装置，其特征在于，所述仿真预测装置包括：

形成模块，用于形成连续排列的多个球形单元，所述多个球形单元的包络曲面用于表示结构表面；

获取模块，用于获取各个所述球形单元中所述结构表面的表示区域的法向量和面积；

仿真模块，用于根据各个所述球形单元中所述结构表面的表示区域的法向量和面积，对各个球形颗粒与所述多个球形单元之间的碰撞进行离散元仿真，得到各个所述球形单元吸收的来自各个所述球形颗粒的碰撞能量；

改变模块，用于根据各个所述球形单元吸收的来自各个所述球形颗粒的碰撞能量，改变各个所述球形单元的球心位置，所述多个球形单元的球心位置的改变用于表示所述结构表面的磨损。

10.根据权利要求9所述的仿真预测装置，其特征在于，所述仿真预测装置还包括：

采样模块，用于在所述球形单元的球面上选择采样点；

表示模块，用于将所述多个球形单元的采样点进行拟合，形成多个球形单元的包络曲面。

Images