基于视觉与边界元法的复杂曲面砂带磨削接触状态快速获取

任绪凯，余焕伟，陈仙凤，杜锡勇，王国彪，陈小奇

（1. 绍兴市特种设备检测院，绍兴 312071；2. 绍兴市特种设备智能检测与评价重点实验室，绍兴 312071 3. 天津大学，天津 300072；4. 上海交通大学，上海 200241）

机器人砂带磨削是一种高效、高质量的加工方式，与砂轮或切削加工方式不同，砂带与弹性接触轮形成的加工工具是柔性的，这赋予其加工复杂曲面工件的独特优势[1–2]，但也给砂带磨削过程的控形、控性带来困难[3]。一方面，砂带磨削中的滑擦作用明显[4]，大多数磨削能在磨削区转变为热量[5]，柔性的接触状态使热量的分布不均匀，给磨削温度预测及控制带来困难，进而导致接触区受力大的部位可能出现高温损伤。另一方面，磨削去除量受局部接触状态的影响，然而目前缺少快速获取局部接触状态的方法。

砂带磨削中磨具与工件间的接触一般认为是弹性接触方式[6–8]，解决弹性接触的方法可分为解析法、数值法与数据驱动的方法。解析法是将接触问题归纳为Hertz接触[9]问题或Signorini[10]、Lurie等[11]接触问题。在砂带的抛光加工中，砂带与工件表面之间具有小的接触压力与摩擦力，有学者将其看作是Hertz接触问题并求解接触区的压力分布[12–14]，即

式中，P0为椭球体压力分布中心的压力；a为椭球体长半轴；b为椭球体短半轴。考虑到磨削过程中接触区材料不断被去除，导致接触压力的变化，Wang等[14]进一步提出了停留时间变化下的接触区压力分布，计算公式为

式中，xi为平行长于半轴a的x轴上的点；Pt为接触压力；t为停留时间；vw为工件移动速度。Wang等[15]同样将接触轮的变形简化为Hertz接触问题，平板磨削试验结果表明准确率不低于96.9%。

Signorini接触问题致力于解决弹性–刚性体之间的接触状态问题，与砂带磨削中的磨具–工件接触情形类似。然而，使用积分法对以上两种模型求解仍然是比较困难的，同时将接触问题简化为简单的二维接触问题，对于工件表面轮廓复杂的砂带磨削加工来说，难以满足实际生产中精确控形的需求。

有限元法（Finite element method，FEM）是将连续的方程离散化进而解决复杂的工程问题[16]。在解决多体接触问题上，商业有限元软件多以节点到表面（Node-to-Surface）形式处理物体间的接触问题。对FEM来说，离散化的程度一般决定了最终的计算精度，同时也是计算耗时多少的主要因素。在忽略时间成本的前提下，FEM在解决静态、准静态接触问题上具有比较高的准确度。Schröder等[17]通过改进FEM解决了具有Tresca摩擦现象的Signorini接触问题，砂带磨削水龙头试验结果证明考虑摩擦现象计算得到的接触状态更加准确，虽然FEM能够比较准确地计算三维、任意形状工件的接触状态，但是时间成本较高，一次接触状态仿真计算短则数分钟，长则数小时乃至数天。为解决FEM计算效率问题，Weinert等[18]引入误差自主控制的有限元思想，即在接触区自主划分高密度网格，而非接触区划分粗网格，以此提高计算效率。Blum[19]和Suttmeier[20]等同样引入自适应网格划分技术来提高FEM计算接触状态的效率。但优化后的FEM仍然需要较长的时间来获取接触状态计算结果。考虑到接触状态分析中已有研究只对工件表面的受力状态感兴趣，而内部的应力状态则不是计算目标，采用边界元法（Boundary element method，BEM）则能体现出这种优势[21]，只需要划分、求解工件表面网格节点，极大缩减计算量。BEM需要利用格林公式与基本解建立变形网格与外部力的积分方程，边界条件为在边界处的网格受外力为0，但是接触边界仍然需要迭代的方式计算得到，因此许多研究工作致力于提高迭代的效率或缩减迭代步骤。

考虑到FEM和BEM在快速获得接触状态方面的不足，有学者提出采用数据驱动的技术将繁琐耗时的计算过程采用机器学习模型进行代替[22–23]，仿真得到的接触状态结果用来训练机器学习模型。Zhang等[24–25]将获取接触问题看作是由磨具与工件间的距离信息到接触状态的回归问题，采用SVR建立预测模型，模型的输入为特别选择的接触体间的距离信息。模型的预测误差不超过5%，一次接触状态的预测时间为1 s左右，远高于数值法计算效率（15 min左右）。Lipiński等[26]同样提出了一个基于改进的神经网络模型进行接触状态预测的方法，测试结果表明该方法在预测磨削接触状态上表现优异。

综上，目前国内外在求解砂带磨削中接触状态时，广泛采用简化接触模型并结合数值法通过迭代计算得到结果，此种方法效率低，过度的简化难以满足高精度磨削加工要求；引入机器学习算法建立工件轮廓与接触状态映射模型是一种高效率的得到接触状态的方法，然而可能的接触范围难以界定，容易陷入局部最优。获取大量的、可靠的训练数据集也是建立机器学习模型预测接触状态的难点。此外，未考虑到磨具弹性参数的方法在计算接触状态时，实用性难以保证。因此，为了高效准确地计算砂带磨削接触状态，本研究引入视觉获得接触区域宽度，结合BEM快速计算得到接触状态。

1 基于数值法的接触状态计算

1.1 基于FEM的接触状态计算

要想获得磨削复杂曲面工件的高精度接触状态，需要建立磨具与工件的三维模型。考虑到工件轮廓的复杂性，无法将模型简化为二维进行数值求解，在ABAQUS有限元分析软件平台建立三维模型并求解。模型如图1所示，工件与磨具的接触姿态通过调整工件的角度确定。磨具由多孔PU轮与砂带组成，可以看作为具有确定力学参数的弹性体，力学参数通过简单试验对比FEM模型得到。最终得到的FEM设置参数如表1所示。

表1 有限元模型参数Table 1 Parameters of finite element model

图1 有限元模型Fig.1 Finite element model

1.2 基于BEM的接触状态计算

磨削过程中的接触是动态的、时变的，对计算的速度要求很高，FEM难以满足加工过程中接触状态快速获取的要求。在FEM的基础上，BEM被发展起来。与FEM不同的是，BEM仅需要将物体表面离散化，内部无网格、无节点，因此极大缩减了计算量。忽略摩擦力与切向应力，如图2所示，Love等[27]分析了均匀法向压力下，范围为–a≤x≤a，–b≤y≤b平面内任意一点p（x，y，z）在z方向上的位移与点A、B、C、D的关系为

图2 均匀分布法向压力下矩形内点法向位移Fig.2 Normal displacement of point in rectangular under uniform normal pressure

当上述结论推广到三维曲面时，如图3所示，在点（i，j）处的变形受到点1～4的影响，同时点1～4又受到其周围点的影响。因此对于砂带表面接触区域，变形的表面是一个整体，某点的变形服从于整体曲面变形的平衡。基于此，将曲面划分为宽度为c，大小为m×n的网格区域。由式（3）可得，应变与应力的矩阵表达式[28]为

图3 离散曲面示意图Fig.3 Diagram of discreted surface

式中，为点（i，j）处的法向位移；pij为点（i，j）的应力；为参考点为（i，j）时的全局影响系数，由式（3）可知

式中，A为（m×n）×（m×n）影响系数矩阵。

在实际接触问题中，砂带表面接触区域是未知的，即m、n未知，接触问题的解需要通过不断迭代来满足边界条件。在接触区域两物体接触间距为0（即在接触区，弹性体变形是已知的，服从于刚体表面轮廓）。在两物体接触区以外，压应力为0，而弹性体的位移，一般来说不为0。在开始迭代前，需要假设一个可能的砂带表面接触区域，应力与位移可以分为接触区内应力pi与接触区内位移ui，接触区外应力po与接触区外位移uo。若此状态成立，则ui已知，且po= 0。依据式（6），则有

进一步的，

接触状态pi最终可由式（8）得到。同时，砂带表面接触区外的弹性体变形可通过式（9）计算得到。通常，在第1次迭代时可能在接触区域内产生一个负的压应力（即拉应力），同时在接触区域外产生一个不合理的负分离距离。第2次迭代选择新的接触区域，使所有产生拉应力的接触点被消除，同时产生负分离的点被考虑为接触点。重复第2步直到迭代至一个合理的接触状态（没有受拉应力的点，没有负距离的点）。若没有给定接触区域范围的参考，上述迭代将是一个计算量很大且耗时的过程。

2 基于视觉的接触状态快速获取

与FEM相比，尽管BEM缩减了网格数量，但在精度要求高的计算过程中仍要求解大量的非线性方程组，迭代过程将成倍地增加计算时间。为进一步提高BEM计算接触状态效率，本研究提出采用基于视觉的方式得到接触区宽度b作为BEM的边界条件，即b+BEM方法，减少迭代求解大量非线性方程组计算量，一次计算得到整个接触区的接触状态。

2.1 试验设置

接触区域与磨削中的其他特征不同，其大小难以用声音、磨削力来映射。工件在与砂带接触过程中的接触区域也难以使用相机直接拍摄。但是参与磨削的磨粒顶部会裸露出光亮面，基于此特点，可采用磨削前对砂带进行着色+磨削后对砂带磨痕抓拍的方式得到实时的接触区宽度。

工件为3D激光选区熔覆成形的复杂曲面Inconel 718工件，如图4（a）所示，采用KEYENCE，LJ–7200线激光扫描仪获取扫描面点云数据，将点云数据作为曲面工件离散数据，扫描密度为0.2 mm。实时获取接触区宽度试验设置如图4（b）所示，工业相机+辅助超亮LED光源保证相机在50 μs的曝光时间下能够采集到运行中砂带的清晰图像。工业相机型号为Basler acA4112–8gc，曝光时间最小为30 μs，传感器为CMOS，分辨率不低于4112×3008，可采集彩色或黑白图像，采样率设置为4帧/s。其中砂带被黑色涂料粉刷，涂料为水溶性颜料，不会造成砂带空隙堵塞。粉刷位置为整条砂带的表面，可使磨削后磨痕清晰。需要特别说明的是，实际磨削过程中磨削宽度是时变的、没有规律的，新的磨削痕迹可能被前面磨削道次痕迹所覆盖，需要在相机与磨削区之间增加喷涂装置实时涂敷掉前面磨削痕迹。本试验为了简化试验装置，通过试验参数的特别设置，保证随着磨削的进行试验中接触区域是逐渐增大的，因此在本试验中识别到的最大砂带表面接触区宽度即为当前工件接触区宽度，不受前面磨削过程产生的痕迹影响。

图4 获得接触区宽度试验设置Fig.4 Experimental setup to obtain contact width

2.2 基于图像的接触区域获取

基于图像在线获取接触区宽度流程如图5所示。图像尺寸标定是为了获得相机视场下的实际尺寸，参照物为已知宽度的砂带，如图6（a）所示。图6（b）为50 μs曝光时间下采集到的磨削过程中砂带图像，可以看出参与磨削的磨粒在图像中呈明显的光亮色，而未参与磨削的区域在黑色涂料染色下为黑色。图6（c）为对图6（b）进行二值化处理，二值化阈值设置为0.3。砂带表面接触区域外的少量噪点是由反光现象造成的，噪点随机分布且噪点之间相隔较远，因此对采集到的原始图像进行聚类去噪，k参数设置为4。为方便得出磨削区域宽度b，对去噪后的图片进行线性膨胀形态学处理，线性膨胀方向为砂带运行方向。经过上述步骤处理后结果如图6（d）所示，其中的磨削区域更加清晰，在选取宽度内测量白色区域平均宽度，得到实际的工件接触区宽度。

图5 实时获取接触区流程Fig.5 Procedures to obtain contact area online

图6 基于图像处理获取接触区宽度Fig.6 Procedures to obtain contact area based on image processing

因为砂带在与接触轮配合形成磨具时受到一定的张紧力，同时自身强度较高，因此忽略因转动和切向摩擦力带来的切向应力导致的接触状态变化。如图7所示，在已知工件与磨具接触姿态的条件下，由轮廓数据可得二者在点pij处的法向距离信息为

图7 磨具与工件接触姿态示意图（mm）Fig.7 Diagram of contact posture between grinding tool and workpiece (mm)

式中，R为磨具半径。

如图8（a）所示，总的可能接触区域内磨具与工件的法向距离集合为

若接触区域宽度b已知，则可通过等高点数据寻优与最近邻搜索获得具体接触区域，寻优算法采用BADS（Bayesian adaptive direct search）。由接触区宽度获得接触区域的伪代码如下。

第1步：由工件轮廓、磨具轮廓和接触姿态确定可能的接触区域间距，确定可能的接触区间g。

第2步：由图像信息已知接触区宽度。

第3步：寻优。

BADS（@targetfun，hmin）：

targetfun：

搜索集合g内距离范围为hmin–0.1～hmin+0.1的点集合g1；

搜索g1集合内y轴数值最小的点ymin；

搜索g1集合内y轴数值最大的点ymax；

可能的接触区宽度Dis=|ymax–ymin|；

返回与真实接触区宽度b的差方D=（b–Dis）2

结束；

判断D是否满足优化准则；

是则返回hmin；

否则重复第3步；

结束。

第4步：输出接触区域集合g1。

其中，hmin为工件接触区域边缘在初始状态下的接触体法向距离；等高线宽度范围设置为0.1 mm；由于BADS在快速寻优问题中的优势，将其作为获得hmin的优化算法。假设b=20 mm，最终得到的工件接触区域如图8（b）所示。

图8 接触区域示意图Fig.8 Diagram of contact area

3 计算与验证

在已知工件接触区与轮廓离散点云的基础上，由式（8）可计算得到接触状态，将所得结果与磨削试验结果和FEM结果对比，验证本研究提出的接触状态获取方法的准确性与效率。在进行接触磨削试验前对工件表面涂色，以便清晰地看出工件接触区域。接触磨削试验参数如表2所示，为保证工件表面接触区域上的涂色在被去除的同时不会磨削过多的母材导致接触区域的变化，设置一个非常低的砂带速度。为模仿实际磨削的复杂工况，接触姿态在x、y、z轴方向分别进行了–5°、–15°和–20°的旋转。

表2 获取磨削接触状态试验设置Table 2 Experimental settings to obtain grinding contact status

磨削试验结果、FEM计算结果与b+BEM获得的接触状态结果如图9所示。图9（a）为接触磨削试验结果，被磨掉的涂色部分为接触区域，其中沿y轴的磨削区宽度b可由垂直划痕的磨削区域最宽处测得，约为12 mm。图9（b）为FEM模型获得的接触状态，接触区宽度为11.23 mm，与接触磨削试验结果基本吻合，说明FEM模型合理，接触区最大压力为3.05 MPa。图9（c）为b+BEM获得的接触状态示意图，其中接触区宽度由接触磨削试验获得，b=12 mm，接触区最大压力为3.22 MPa。对比接触区域形状，FEM与b+BEM获得的接触区域和磨削试验获得的接触区域相似。

图9 不同方法获得接触状态Fig.9 Contact status obtained by different methods

对比数据如表3所示，FEM获得的接触状态作为对照组，b+BEM获得接触状态为试验组。所有数据处理与仿真试验在个人电脑（CPU i5 7400，RAM 16G）上完成。FEM获得接触状态在ABAQUS软件中完成，接触区域附近网格尺寸精度最高为0.5 mm，接触状态的一次计算时间约为25 min。由于b+BEM的离散基于点云数据，因此网格尺寸精度与点云数据密度相关，都为0.2 mm。在采用b+BEM计算接触状态时，每次接触状态的计算约需要8 s，优于FEM计算效率，但是还达不到实时要求。造成b+BEM计算效率过低的原因是工件点云中与磨具接触的点需要一个寻优的步骤获得，当点云范围很大时，寻优的步骤增多，增加了计算成本。选取可能的接触区域点云集合，缩减感兴趣点云数据量，可以达到提高计算效率的目的。选取规则为以最邻近磨具的点为中心，宽度±（b+5），垂直于y轴的矩形带内的点云数据，选取结果如9（c）所示。在进行数据量缩减后，b+BEM计算一次接触状态的时间缩短到1 s内。

表3 不同方法获得接触状态对比Table 3 Comparisons of different methods to obtain contact status

在最大接触压力与接触面积方面，b+BEM获得的接触状态与FEM结果吻合度较好，接触区最大压力与接触面积偏差都没有超过5.6%，一方面是因为接触区宽度由实际试验得到，保证了接触区域的准确性；另一方面是因为BEM在计算接触状态时将接触区看作一个整体，每个点的压力计算都会将其他接触点的影响计算在内，最终得到比较准确的接触状态。

综上，在已知材料参数的条件下，FEM获得的接触状态准确度较高，但计算成本也很高，难以满足实际生产活动中对于后续局部材料去除量与热流分布快速获取的需求，但可将其作为参照组来验证其他可行的获取接触状态的方法。b+BEM获取接触状态的方法需要额外的辅助得到接触区宽度，但得到的结果可信度高，采用选取可能接触区域，缩减点云数据量的方式提高计算效率，能够快速计算局部接触状态。未来，可通过升级计算机硬件、优化网格密度、改进寻优算法等手段进一步提高计算效率。因此，采用b+BEM获取接触状态的方法在磨削应用中具有明显的优势。

4 结论

本研究首先分析了FEM与BEM在接触状态获取方面的应用，提出了基于视觉获取接触区宽度并结合BEM，即b+BEM方法，实现了基于点云轮廓的局部接触状态快速计算。与数值法相比，该方法可计算复杂曲面工件接触状态，效率与适用性更高；与机器学习模型相比，该方法具有明确的物理意义，不会陷入局部最优且不需要额外的训练样本。计算结果与FEM仿真结果对比，b+BEM计算获得接触状态的误差不超过5.6%，而且在计算效率上远高于FEM。结合点云数据量缩减规则，该方法能够在1 s内获得18 mm×18 mm，精度0.2 mm区域内的接触状态，为后续计算局部材料去除量与热流分布奠定基础。