一种虚拟激光加工过程几何仿真方法

摘要

本发明公开了一种虚拟激光柔性加工过程几何仿真方法。该方法包括：1)虚拟加工平台的建立；2)基于CAD数据的轨迹规划；3)轨迹优化；4)基于加工头效应的三角简化算法；5)碰撞检测模型。利用本发明可在计算机可视化环境中观察到加工运动过程，直观显示激光加工机器人与其它设备及工件的碰撞干涉情形及不合理加工轨迹，及早预见加工缺陷，为工艺参数调整及工艺优化提供有效依据，节省人力、物力、时间与成本。

说明书

技术领域

本发明涉及一种仿真方法，更具体地，本发明涉及一种激光柔性加工过程几何仿真系统及实现方法。

背景技术

虚拟制造是集现代制造工艺、计算机图形学、并行工程、人工智能、人工现实等多种高新技术为一体，对所有的制造环境及制造活动进行建模与仿真，是对产品全生命周期的仿真。其中的加工过程仿真是虚拟制造的底层与核心。由于加工过程中涉及因素多，过程复杂，而且对加工过程的研究难以集成到虚拟制造的大仿真系统下，这一直是虚拟制造的瓶颈问题之一。而加工过程几何仿真可为激光制造系统提供完备的仿真环境、加工过程轨迹验证等重要功能环节，对提高系统运行效率、消除工艺设计缺陷提供支持。

集成化激光制造系统是由计算机控制系统在同一五轴框架式机器人平台下实现对不同种类的加工方式(切割、焊接、表面处理等)的一种自动化系统。在现有方法中，几何仿真均是针对某种特定类型的加工方式，如华中科技大学机械学院基于分形扫描的选择性激光烧结过程动态几何仿真但这种方式存在仿真应用范围局限性大、可维护性差等缺点；也有与物理模型集成对激光加工过程进行完整仿真，如中国科学院力学所张桃红博士提出虚拟激光柔性加工的完整过程仿真进行了探索。然而，其中几何仿真使用商业化IGRIP软件，具有开放性、可扩展性差，而且软件独立性差(软件只能在配有相关硬件的条件下使用)，难于在不同的系统之间移植。

如果以激光制造系统中框架式机器人为载体，采用模块化设计并使用VC编程环境和OpenGL图形函数库构建虚拟激光加工环境，数据接口、轨迹规划、图元排序以及碰撞检测等模块嵌入到该可扩展的虚拟环境中实现仿真。该方法可以克服上述难题，但是目前尚未见相关的研究报道或者是公布的专利技术。

发明内容

本发明的目的是在虚拟制造的仿真大系统下，建立包括虚拟环境与虚拟设备的虚拟加工平台，分析加工轨迹规划优化方法及算法实现，建立基于三角网格简化的碰撞检测模型，嵌入到虚拟加工平台中。

为了达到上述目的，本发明采取如下技术方案：

一种嵌入式激光柔性加工过程仿真方法，其步骤包括：

(1)建立虚拟加工平台：在虚拟加工环境中加入虚拟设备即虚拟激光加工机器人；虚拟激光加工机器人的建立包括几何模型的建立和控制设备几何动作的运动学模型的建立。几何模型指加工机器人的CAD实体模型；运动学模型包括其正解与运动学逆解，正解指由机器人各轴的运动量计算得加工头处的位姿，即坐标值与法向量值，逆解指由加工头处的位姿反算出机器人各轴的运动量值，反映的是机器人各关节轴值与加工头位姿间的关系。

(2)基于CAD数据的轨迹规划：在虚拟激光加工环境中，激光加工机器人加工轨迹驱动仿真加工过程，它是虚拟环境中加工动作的根本和优化对象。而虚拟仿真系统中待加工工件的面型数据主要来源于除CAD模型，针对大型汽车覆盖件冲压模具表面形貌复杂的特点，从工艺力学分析的角度出发，结合大量的工程实践经验，从中抽取出一些有共同几何特征的典型型面——棱脊，并基于STL数据格式实现加工轨迹规划，便于驱动仿真进一步进行。

(3)轨迹优化：加工点1，2，...，n，采用C空间法计算任两点i，j(i＝1，2，...，n；j＝1，2，...，n)的无碰撞最短轨迹；如图2所示，加工点i用O₁表示加工点j用O₂表示，碰撞点为p，建立通过p点和j点的连线(图中虚线表示)的w个平面与模具型面相交得到w条交线，图中用实线表示出其中的两条，取w条交线中的最短交线，即为加工点i，j之间的无碰撞最短轨迹；w的取值大小决定了计算量大小，w的取值越大，计算量越大，取得的最短轨迹也越准确；w的取值至少为8；两点轨迹算出后采用蚂蚁算法计算加工这n个加工点的最优加工轨迹；

(4)三角简化算法：针对加工过程仿真干涉检测的特点，设计了一种基于网格覆盖，误差控制及加工头尺寸效应的三角网格简化算法。该算法首先对三角网格寻找一种初始的简化网格覆盖T，给出相应的初始覆盖网格点集V。进一步在T的各网格区域内寻找与相应网格距离最大的网格点，并根据设定的最大及最小误差范围及加工头的尺寸决定该点是否保留。该算法能够在控制误差的前提下快速有效的简化三角网格，特别适用于加工过程仿真干涉检测应用的三角网格简化。

(5)碰撞检测模型：对于虚拟加工环境，可视为由n个形体(每个形体由m个三角网格组成)构成的集合，若对每个形体直接实施加工头(加工头由k个三角网格组成)与周围环境形体的碰撞检测算法，则进行一次碰撞检测需n×m×k次判断，显然算法的效率十分低下。多数碰撞检测算法均存在时间步长问题，即碰撞检测的精度与频率问题，频率过高，则计算量大，增加了系统负担，频率过低，则可能错过碰撞检测。为此，在虚拟加工环境中建立了加工头包围球，利用包容球在整个加工轨迹形成的扫略体与加工件实体三角网格进行求交计算。该算法可避免时间步长问题，并提高碰撞检测算法的效率。

在上述方法中，使用VC编程开发环境和OpenGL图形函数库建立了具有通用性、易维护性和可扩展性的计算机模拟系统。

与现有技术相比，本发明的有益效果是：

在虚拟激光加工环境中直观显示激光加工机器人与其它设备及工件的碰撞干涉情形及不合理加工轨迹；及早预见加工缺陷，为工艺参数调整及工艺优化提供有效依据，节省人力、物力、时间与成本。

附图说明

图1(a)是加工机器人的几何模型整体结构图；

图1(b)是图1(a)中虚线圈中的腕部加工头的放大图；

图2是棱脊示意图；

图3是系统总体结构流程图

图4是本发明加工点最优加工路径的蚂蚁算法流程图；

图5是三角网格简化算法流程示意图；

图6是顶点P及其所在三角网格的投影示意图；

图7是包容球扫略体碰撞检测示意图。

具体实施方式

下面结合附图和具体实施方式对本发明作进一步详细描述

本发明的方法主要是仿真过程中的建模、轨迹规划优化及碰撞干涉检测模型，其具体步骤如下：

1.建立虚拟激光加工平台：该虚拟平台是为所有激光加工提供的加工平台，所要虚拟的是激光加工设备和虚拟环境。本实施例针对由中国科学院力学研究所研制的“一种具有柔性传输和多轴联动的激光加工装置”，专利号为98101217.5；其虚拟环境主要是加工间灯光、辅助设施及布局等。虚拟加工设备主要是虚拟激光加工机器人，包括机器人几何模型和可执行任何轨迹的机器人运动学控制模型。如图1(a)所示，加工机器人为5轴框架式机器人，有x移动轴1、y移动轴2、z移动轴3和A转动轴5、C转动轴6，几何模型如图1所示，其框架4尺寸为5.77m×3.63m×2.0m，硬限位为x：4.45m；y：2.755m；z：1.085m，运动学关系矩阵为：

$T_5=\left(\begin{array}{cccc}1& 0& 0& x\\ 0& 1& 0& y\\ 0& 0& 1& z\\ 0& 0& 0& 1\end{array}\right)·\left(\begin{array}{cccc}{C}_4& -S_4& 0& 0\\ S_4& C_4& 0& 0\\ 0& 0& 1& 0\\ 0& 0& 0& 1\end{array}\right)·\left(\begin{array}{cccc}1& 0& 0& 0\\ 0& 1& 0& {-\mathrm{ll}}_1\\ 0& 0& 1& 0\\ 0& 0& 0& 1\end{array}\right)·\left(\begin{array}{cccc}{C}_5& 0& S_5& 0\\ 0& 1& 0& 0\\ -S_5& 0& C_5& 0\\ 0& 0& 0& 1\end{array}\right)$

$·\left(\begin{array}{cccc}-1& 0& 0& 0\\ 0& 1& 0& 0\\ 0& 0& -1& 0\\ 0& 0& 0& 1\end{array}\right)$

其逆解为：

P_x＝x+S₄ll₁+a_x·T_l   P_y＝y-C₄ll₁+a_y·T_l  P_z＝z+a_z·T_l

A＝arctan(a_y/a_x)     C＝arccos(-a_z)

式中P_x，P_y，P_z，A，C为各轴的运动量值；x，y，z为腕部加工工具头7的坐标值；法向矢量(a_x，a_y，a_z)是机器人腕部加工工具头7的姿态，如图1(b)所示，腕部加工工具头7的三个方向余弦为姿态)；T_l是腕部加工工具头7的长度；ll₁是腕部手臂8的长度；θ_A代表A转动轴转动角度；θ_C代表C转动轴转动角度；S₄＝sinθ_A代表A转动轴转动角度的正弦计算；C₄＝cosθ_A代表A转动轴转动角度的余弦计算；S₅＝sinθ_C代表C转动轴转动角度的正弦计算；C₅＝cosθ_C代表C转动轴转动角度的余弦计算。

上述运动学关系通过软件IGRIP的Shared Library二次开发功能导入给虚拟机器人形成机器人运动学控制模型。

2.基于STL格式的棱脊轨迹规划：棱脊是加工表面中常见的一种面型，但因为应用环境，应用目的不同，加之棱脊的种类很多，所以还没有一个通用的定义。从特征线的角度以及激光强化加工的目的出发，可以构造性地给出如下的棱脊定义：给定一条光滑的空间曲线L，则在L上任一点可做相应的法平面。在法平面上做过该点的曲线R，令R沿L运动，R扫掠所构成的曲面即称为棱脊。曲线L称为棱脊导线，R称为棱脊母线(如图2所示)。

由上述定义可知，导线与母线就是棱脊的两条特征线，知道一条棱脊导线以及导线上某一点处的母线，则整个棱脊就可以确定下来。基于STL棱脊轨迹规划的流程如下：

(1)STL文件中，实体表面的所有三角网格被无序地列出，不存在任何拓扑信息，然而在基于STL轨迹规划中，建立三角网格拓扑信息是首要的。

(2)对于复杂表面STL文件，实体表面以大量甚至是海量的三角网格数据记录的，对于零件表面棱脊这类特殊形貌提取特征点。由曲面论可知，零件的棱脊线处曲面的曲率较大。为计算每一顶点的曲率，可在顶点P处建立曲面S(u，v)＝(u，v，h(u，v))，其中h(u，v)＝au²+buv+cv²，曲面再P点得局部形状可由Darboux结构D(p)＝(p，m₁，m₂，N，k₁，k₂)完全描述。V_j(1≤j≤m)在局部坐标系(P_huv)下的坐标值为(u_j，v_j，h_j)，由m个邻点得到的线性方程组为式，用最小二乘法解此方程组，即求得曲面S(u，v)。在求得S(u，v)的基础上，利用曲面的第一、第二基本公式可以得到曲面上P处的法曲率k：

(3)特征点确定后，对特征点集范围内搜索构造样条曲线生成棱脊导线和棱脊母线。

(4)按照一定的步长，分别在棱脊导线和已知母线上取等距点。并可在此基础上进一步计算已知母线上各等距点的法向量。

(5)在导线等距点上都有相应的母线，母线上的等距点(包括坐标及对应法向量)可由上面已知母线上的对应等距点经一定的坐标变换得到。

(6)将上面的不同母线上的对应等距点按照一定顺序排列，即得到所需的加工轨迹。

3.轨迹优化：激光加工的轨迹优化与加工仿真结合起来才能保证正确加工及加工效率。激光加工点的轨迹优化模型为求单次遍访各加工点的最短回路公式表述如下：

$\min \underset{i\ne j}{\Sigma }{d}_{\mathrm{ij}}{x}_{\mathrm{ij}}$

$\underset{i\ne j}{\Sigma }{x}_{\mathrm{ij}}=1,$i，j＝0，1，…，n

$\underset{i,j\in S}{\Sigma }{x}_{\mathrm{ij}}\le \left|S\right|-1,S⋐\{\mathrm{0,1},...,n\}$

x_ij∈{0，1}，

其中，|S|表示集合S中所含元素的个数，n为激光硬化加工的总点数，d_ij为加工点i到加工点j间的无碰撞最短轨迹长度。

加工点i(图4中用O₁点表示)到加工点j(图4中用O₂点表示)间的无碰撞轨迹长度在C空间中求得。如图4所示，取工件上O点为坐标参考原点，加工点O₁相对于工件基础坐标原点O点的坐标为(272.2，24.3，14.7)，加工点O₂的相对坐标为(283.4，44.6，29.8)，从O₁到O₂点发生碰撞的点为P(279.4，34.4，19.8)，以PO₂为连线建立不同的8个面与工件型面相交得到8个交线，其中长度最短的交线即为无碰撞最短轨迹；本实施例建立不同的8个面，还可以取更多个面以获得更准确的无碰撞最短轨迹。

两两点无碰撞轨迹得知后，用蚂蚁算法求单次遍访各加工点的最短回路。设有n个加工点，有m个搜索器，本实施例中m＝20。搜索器在走过的路径上留下搜索标记量，这个量随着时间而更新，搜索器k在t次迭代中从加工点i选择下一加工点j时，根据转移概率来选择，见公式(1)中P^k_ij(t)；其中τ_ij(t)为第t次迭代搜索路径地图上加工点i到加工点j轨迹的搜索标记量，迭代总次数为T，例如T＝10，η_ij为搜索器搜寻法则，在这里η_ij＝1/d_ij，G_k(i)为加工头尚未到达过的加工点集合，α、β表示相对重要性。公式(2)中Δτ_ij^k(t)为第k个搜索器于第t次搜索中加工点i到加工点j的单位轨迹长度搜索标记数，Q为单个搜索器的搜索标记总量，L_k为第k个搜索器搜索路径的总距离，公式(3)中为搜索标记量的更新，ρ为轨迹衰减度，公式(4)中Δτ_ij(t)表示轨迹上所有搜索器的搜索标记量。算法公式为：

$P_{\mathrm{ij}}^k\left(t\right)=\left(\begin{array}{cc}\frac{{\left[{\tau }_{\mathrm{ij}}\left(t\right)\right]}^{\alpha }{\left[{\eta }_{\mathrm{ij}}\right]}^{\beta }}{\underset{h\in G_k\left(i\right)}{\Sigma }{\left[{\tau }_{\mathrm{ih}}\left(t\right)\right]}^{\alpha }{\left[{\eta }_{\mathrm{ij}}\right]}^{\beta }}& j\in G_k\left(i\right)\\ 0& j\notin G_k\left(i\right)\end{array}\right)---\left(1\right)$

τ_ij(t+1)＝(1-ρ)τ_ij(t)+Δτ_ij(t)(3)

${\mathrm{\Delta \tau }}_{\mathrm{ij}}\left(t\right)=\underset{k=1}{\overset{m}{\Sigma }}{\mathrm{\Delta \tau }}_{\mathrm{ij}}^k\left(t\right)---\left(4\right)$

i，j＝0，1，...，n；

k＝1，2，...，m；

算法中α、β、ρ等参数对算法性能有很大的影响。α值的大小表明留在每个加工点上的标记量受重视的程度，α值越大，搜索器选择以前经过的路线的可能性越大，但过大会使搜索过早陷于局部最小解；β的大小表明启发式信息受重视的程度，β值越大，搜索器选择离它近的加工点的可能性也越大；ρ表示轨迹标记量的保留率，如果它的值取得不恰当，得到的结果会很差，计算中采用最佳参数：α＝1，β＝5，ρ＝0.5，Q＝100。如图3所示为上述算法流程图。

4.基于加工头效应的三角网格简化算法：设加工头的特征尺寸为d_t，三角网格三边的边长分别为L₁，L₂，L₃，如果满足：L_i＜d_t，i＝1，2，3，且三角网格投影域内没有高出三角网格面的顶点，则从干涉检测的角度而言，该三角网格已没有必要再加以细分。这就是所谓的加工头尺寸效应。即当三角网格的尺寸相对于加工头的特征尺寸很小时，往往就会出现网格冗余，特别是在一些凹陷区域。所以可针对不同的加工头形状，选择合适的加工头特征尺寸，用以过滤一些满足上述条件的三角网格及其顶点。

基于上面的数据结构及所定义的运算和判断法则，算法流程如图5，具体描述如下：

第一步：设定所需的三个指标：点到三角网格面的最大距离指标D_max；点到三角网格面的最小距离指标D_min；加工头特征尺寸D_t。

第二步：在原三角网格的基础上构造初始的三角面片覆盖。首先将三角网格投影到水平面上(假定为XY面)，搜索得到投影多边形的各个顶点，并确定它们所对应的三角网格顶点，记这些网格顶点为v₁，v₂，...v_n；继而得到三角网格顶点中投影点位于上述投影域中部的某点v_c，连接该点与v₁，v₂，...v_n，然后顺次(顺时针或逆时针)连接v₁，v₂，...v_n，从而构成一初始的三角网格覆盖T_c。T_c中所有的顶点构成初始顶点集V。判断某顶点的投影是否为投影区域边缘点的方法如下：给定顶点p，找到所有以点p为顶点的三角网格，将这些网格投影到水平面上，p′为点p对应的投影点。如图6所示，如果此时投影中所有不以点p′为端点的边构成一个封闭曲线，那么点p′就是投影域的一个内点(见图6(a))；反之，如果这些边不能构成一个封闭曲线，则p′就位于投影域的边缘上(见图6(b))。

第三步：对三角网格覆盖T_c中的某一三角网格t∈T_c，利用前面判断一点是否在三角形区域内方法判断是否有某些顶点的投影落入t的投影域内。如果没有，直接返回，执行对T_c中其它三角网格的判断；如果有，找到这些顶点中与t距离最大的点，记为p_max，相应的最大距离记为d_max。类似的找到与t距离最小的点p_min，相应的最小距离记为d_min。另外以L₁，L₂，L₃分别标记t的三条边的边长。然后做如下的判断：

(1)当|d_max|＜D_max，且|d_min|＜D_min时，返回，执行对T_c中其它三角网格的判断。

(2)当|d_max|＜D_max，且L_i＜D_t，i＝1，2，3时，返回，执行对T_c中其它三角网格的判断。

(3)对于其他情形，将该顶点添加到顶点集V；

第四步：对初始三角面片覆盖中的每一三角面片进行第三步，从而得到一新的顶点集V，对V进行三角剖分，得新的三角网格覆盖T_c。

第五步：对顶点集V和三角面片覆盖T_c重复第三步和第四步。

第六步：当顶点集V不再发生变化时，终止。此时经三角剖分所得的三角面片集T_c即为所求的三角网格。用N_v表示此时V中的顶点个数，N_t表示T_c的三角网格个数。

5.碰撞检测模型：在激光加工中，轨迹规划源于测量数据和CAD数据这两种方式。对于测量数据的规划是首先通过选取代加工区域，再利用智能测量特征点集，最终对离散点构造样条曲线生成加工轨迹；而对于CAD源数据，轨迹规划是先从实体中选取特征点，在此基础上构造样条曲线生成加工轨迹因此，这里加工轨迹可以使用累加弦长三次参数样条曲线描述如下：

P(s)＝[x(s)  y(s)  z(s)]

s为累加弦长，可表示为：

s₀＝0

$s_k=\underset{j=1}{\overset{k}{\Sigma }}{l}_k=\underset{j=1}{\overset{k}{\Sigma }}|P_{j-1}-P_j|$

$=\underset{j=1}{\overset{k}{\Sigma }}\sqrt{{(x_j-x_{j-1})}^2+{(y_j-y_{j-1})}^2+{(z_j-z_{j-1})}^2},$

k＝0，1，...，n

其中，

$\left(\begin{array}{c}x\left(s\right)=F_{j-1}\left(s\right)x_{j-1}+F_j\left(s\right)x_j+G_{j-1}\left(s\right)m_{x,j-1}+G_j\left(s\right)m_{x,j}\\ y\left(s\right)=F_{j-1}\left(s\right)y_{j-1}+F_j\left(s\right)y_j+G_{j-1}\left(s\right)m_{y,j-1}+G_j\left(s\right)m_{y,j}\\ z\left(s\right)=F_{j-1}\left(s\right)z_{j-1}+F_j\left(s\right)z_j+G_{j-1}\left(s\right)m_{z,j-1}+G_j\left(s\right)m_{z,j}\end{array}\right),$

$\left(\begin{array}{c}{F}_{j-1}\left(s\right)=\frac{{(s_j-s)}^2[2(s-s_{j-1})+h_j]}{h_j^3}\\ F_j\left(s\right)=\frac{{(s-s_{j-1})}^2[2(s_j-s)+h_j]}{h_j^3}\\ G_{j-1}\left(s\right)=\frac{{(s_j-s)}^2(s-s_{j-1})}{h_j^2}\\ G_j\left(s\right)=\frac{{(s-s_{j-1})}^2(s_j-s)}{h_j^2}\end{array}\right)\left(\begin{array}{c}{h}_j=s_j-s_{j-1}\\ j=\mathrm{1,2},...,n\end{array}\right)$

如图7所示，当包围球沿加工轨迹运动时，三角形所在平面最先与球发生碰撞的点称为接触点，球与三角形发生碰撞的点位碰撞点。那么球与三角形是否碰撞问题可转化为以下两个子问题：判断接触点是否在三角形内与求出碰撞点。对于接触点在三角形内的情况，碰撞点即为接触点利用球正下方点代入式可以求解得到三角形所在平面碰撞点。对于接触点在三角形外的情况，为了简化计算，可以等分s用等分点作为端点的多线段来近似样条曲线，并求各线段到三角网格边的距离，如果大于包容球半径则不发生碰撞；如果小于等于半径则可判定发生碰撞。