CN102981454B - 一种基于势能场及能量泛函优化的刀具轨迹路径产生方法 - Google Patents

Abstract

本发明公开了一种多轴刀具建立基于能量泛函优化的刀路轨迹计算模型，包括能量泛函、边界条件和行宽约束条件，并采用有限元方法获取刀路轨迹优化计算模型基于三角形离散网格的数值近似计算公式，根据刀路轨迹计算模型得到势能场Ф在每个节点I上的近似计算值φ_I，截取火山口的等高线，并向加工区域Ω投影，形成等势轮廓线近似轨迹，对得到的等势轮廓线近似轨迹进行平滑和采样，得到从外向内偏置的刀具轨迹轮廓基准线，在刀具轨迹轮廓基准线之间进行线性插值，形成一段螺旋刀轨，然后将所有的螺旋刀轨连成一条螺旋轨迹，并对螺旋轨迹进行曲线拟合，以形成几何二阶连续的B样条曲线刀轨。本发明可实现高速加工、优化切削力分布，并提高表面加工质量。

Description

一种基于势能场及能量泛函优化的刀具轨迹路径产生方法

技术领域

本发明属于铣削加工领域，更具体地，涉及一种基于势能场及能量泛函优化的刀具轨迹路径产生方法。

背景技术

高质量、高效的数控加工要求计算机辅助制造系统（Computer AidedManufacturing，简称CAM）产生的刀具轨迹具有光滑和行距均匀（即稳定的材料去除率）的特点。利用传统的参数曲面等参线或者轮廓等距偏置方法产生的曲线刀轨在一定程度上能有效控制行距（包括行切，环切等），但其形状严重受制于加工区域的参数构成或边界形状，光滑性不好；新发展的曲线轨迹方法如螺旋线加工在加工区域产生类似于电势场的等势轮廓线族，然后通过离散、修正的方式构造尽可能均匀分布的偏置轮廓线，在此基础上形成螺旋线轨迹。这种方法在应用到轮廓外形或者加工区域几何形状特别复杂时，各偏置轮廓线之间的距离变化剧烈，难以产生只需一次进退刀并且行距变化稳定的刀轨。

美国专利US6591158B1在加工区域的内外边界设置电荷，并利用计算的电场线作为辅助线来构建偏置轮廓线，进而构造螺旋状刀具轨迹。该方法中对电场线的求解依赖于毛坯与成型零件轮廓上的起始电荷分布。在电荷分布不均的情况下会严重影响所生成的电场辅助线的形状，最终影响刀具轨迹质量，例如造成螺旋刀轨的不均匀。另外由于需要根据电荷分布求电场线，涉及到对一个平面向量，即2个标量的求解，计算相对复杂，且不能对行距进行有效控制。该专利只适用于产生叶片粗加工刀轨，其计算模型要求必须具有内外边界，所以不能运用到无岛屿的型腔如狭长型的空腔加工中。

美国专利US6591158利用求解稳恒电磁场本征值的方法得到偏置轮廓线，进一步产生光滑的螺旋刀具轨迹。但是该方法只适用于无岛屿型腔加工，且同样存在不能控制切削行距的缺点。

申请号为200810207221.0的中国专利采用美国加州大学伯克利分校J.A.SETHIAN提出的波前传播方法构造水平集函数。该发明根据加工要求给定偏置距离，利用求解水平集方程得到型腔边界的偏置线，当遇到岛屿时型腔边界的推进停止。这种方法只能从外边界或者内边界进行偏置，缺乏全局范围内的偏置线分布优化（光滑性）和行距控制；同时在型腔的推进中会自动产生区域分割，形成多次进退刀路径；对于复杂加工曲面来说，由于偏置轮廓线之间的任意可能形状，在使用通用的渐变（morphing）算法构建螺旋线时，也不能完全避免实际切削轨迹的重叠。

中国专利“一种平面螺旋状或环形铣削刀具路径生成方法”发展了前述方法。该发明假设加工区域存在一个由调和方程（拉普拉斯方程）确定的边值问题界，通过有限元方求解该问题，形成光滑的等值线分布。然后按照边界特征进行分区，分区边界为该标量场中的梯度曲线（GradientCurve），通过对最长梯度曲线按照行距要求进行优化布置得到一系列代表等势轮廓基准线的势能值，依据这些势能值形成偏置轮廓基准线并以此形成螺旋刀轨。该方法的优点在于满足了刀轨光滑性需求并进行了行距控制，但行距控制未优化，且光滑性需求和行距控制是先后满足的，没有在一个统一的数学模型中一起优化。

叶片加工、三轴型腔加工等加工区域可能具有任意形状，或具有任意岛屿分布，前述的已有方法都不能产生全局光滑和稳定行距分布的刀轨，其原因在于前述算法都是将偏置轮廓线的生成转化为偏微分方程求解，而无法考虑除能量分布之外的其他约束如行距控制等，即使专利“一种平面螺旋状或环形铣削刀具路径生成方法”进行了有效行距控制，但未将刀轨形状与行距控制一起进行优化。

发明内容

针对现有技术的缺陷，本发明的目的在于提供一种基于势能场及能量泛函优化的刀具轨迹路径产生方法，其能够针对具体加工需求，将约束条件整合到问题求解中，进行统一的优化求解，从而实现刀轨形状与行距控制的同时优化，并可针对任意形状的加工区域产生光滑和行距全局优化的轮廓轨迹，在此基础上形成具有更精准和优化的行距控制且只有一次进退刀的螺旋轨迹曲线，以实现高速加工、优化切削力分布，并提高表面加工质量。

为实现上述目的，本发明提供了一种基于势能场及能量泛函优化的刀具轨迹路径产生方法，包括以下步骤：

（1）为多轴刀具建立基于能量泛函优化的刀路轨迹计算模型，包括能量泛函、边界条件和行宽约束条件，并采用有限元方法获取刀路轨迹优化计算模型基于三角形离散网格的数值近似计算公式，具体包括以下子步骤：

（a）确定多轴刀具的加工区域及其外边界轮廓和内边界轮廓，并用三角形网格对整个加工区域进行网格划分，以构建势能场Ф，其中加工区域及其三角网格划分共具有k个节点，多轴刀具的加工区域上的所有节点集合为Ω，加工区域外边界上的节点集合为Г₁，加工区域内边界上的节点集合为Г₂，Ф_i(i＝1,2,...,k)表示加工区域中任意节点i所具有的势能值，整个加工区域形成一势能场Ф；

（b）为势能场Ф建立能量泛函F(φ)，具体公式为：

$F\left(\mathrm{\φ}\right)\≡\underset{D}{\∫}{\mathrm{du}}^1{\mathrm{du}}^2\sqrt{G\left(u\right)}(\frac{1}{2}{g}^{\mathrm{\α\β}}\left(u\right)\frac{\mathrm{d\Φ}\left(u\right)}{{\mathrm{du}}^{\mathrm{\α}}}\·\frac{\mathrm{d\Φ}\left(u\right)}{{\mathrm{du}}^{\mathrm{\β}}})$

其中α,β为按爱因斯坦求和约定的标号，取值分别为1和2，u¹,u²为曲面参数，u＝(u¹,u²)^T， $g^{\mathrm{\α\β}}\left(u\right)\frac{\mathrm{d\Φ}\left(u\right)}{{\mathrm{du}}^{\mathrm{\α}}}\·\frac{\mathrm{d\Φ}\left(u\right)}{{\mathrm{du}}^{\mathrm{\β}}}=\underset{\mathrm{\α}=\mathrm{1,2}}{\mathrm{\Σ}}\underset{\mathrm{\β}=\mathrm{1,2}}{\mathrm{\Σ}}{g}^{\mathrm{\α\β}}\left(u\right)\frac{\mathrm{d\Φ}\left(u\right)}{{\mathrm{du}}^{\mathrm{\α}}}\·\frac{\mathrm{d\Φ}\left(u\right)}{{\mathrm{du}}^{\mathrm{\β}}},$ $\frac{\mathrm{d\φ}}{{\mathrm{du}}^{\mathrm{\α}}}$ 表示向量G^αβ为曲面度量张量的矩阵，且 $G_{\mathrm{\α\β}}=\left[\begin{array}{cc}{G}_{11}& G_{12}\\ G_{21}& G_{22}\end{array}\right],$ G(u)＝det(G_αβ)表示曲面度量张量G_αβ的行列式计算值，det（）表示求行列式，g^αβ(u)为曲面度量张量的逆，u^α和u^β满足爱因斯坦求和约定规则，表示曲面参数（u¹,u²）；

（c）获取能量泛函F(φ)的有限元数值近似计算公式，采用以下公式

$F\left(\mathrm{\φ}\right)=\underset{(\mathrm{IJK}\∈\mathrm{\Π})}{\mathrm{\Σ}}\frac{1}{2A_{\left(\mathrm{IJK}\right)}}[{S_{\mathrm{IK}}}^2{({\mathrm{\φ}}_J-{\mathrm{\φ}}_I)}^2{S_{\mathrm{IJ}}}^2{({\mathrm{\φ}}_K-{\mathrm{\φ}}_I)}^2-2(S_{\mathrm{IK}}\·S_{\mathrm{IJ}})({\mathrm{\φ}}_K-{\mathrm{\φ}}_I)\·({\mathrm{\φ}}_J-{\mathrm{\φ}}_I)]$

其中（I，J，K）表示加工区域中的一个三角形，三角形网格的顶点为u_I，u_J，u_K，I，J，K分别表示三个顶点的节点序号，u_I，u_J，u_K对应物理空间中的三点S_I=S（u_I），S_J=S（u_J），S_K=S（u_K），S表示从参数空间到物理空间的坐标转换，S_IJ，S_JK，S_KI表示物理空间中三角形的三边可用矢量，分别为S_IJ=S_J-S_I，S_JK=S_K-S_J，S_KI=S_I-S_K， $A_{\left(\mathrm{IJK}\right)}=\sqrt{\det \left(G_{\mathrm{\α\β}}\right)}=\sqrt{{S_{\mathrm{IK}}}^2{S_{\mathrm{IJ}}}^2-{(S_{\mathrm{IK}}\·S_{\mathrm{IJ}})}^2},$ det（）表示求度量张量G_αβ的行列式，φ_I＝φ(u_I)代表势能场Ф在节点I上的值，φ_K＝φ(u_K)代表势能场Ф在节点K上的值,φ_J＝φ(u_J)代表势能场Ф在节点J上的值,п表示全体三角形网格的集合；

（d）利用行距条件，建立行距控制约束泛函以及其数值近似计算公式，采用以下公式：

$\frac{1}{{h_s}^2}-\frac{1}{{A_{\left(\mathrm{IJK}\right)}}^2}({S_{\mathrm{IK}}}^2{({\mathrm{\φ}}_J-{\mathrm{\φ}}_I)}^2+{S_{\mathrm{IJ}}}^2{({\mathrm{\φ}}_K-{\mathrm{\φ}}_I)}^2-2(S_{\mathrm{IK}}\·S_{\mathrm{IJ}})({\mathrm{\φ}}_K-{\mathrm{\φ}}_I)\·({\mathrm{\φ}}_J-{\mathrm{\φ}}_I))\≤0$

(IJK∈П)

φ_I＝0           在边界Г₁上

φ_I＝const＝n    在边界Г₂上

其中h_s表示加工行距；

（e）建立基于能量泛函优化的刀路轨迹计算边界条件；

（2）根据刀路轨迹计算模型得到势能场Ф在每个节点I上的近似计算值φ_I；

（3）按Ф=0，1，…，6，截取势能场的等高线，并向加工区域Ω投影，形成等势轮廓线近似轨迹；

（4）对得到的等势轮廓线近似轨迹进行平滑和采样，得到从外向内偏置的刀具轨迹轮廓基准线；

（5）在刀具轨迹轮廓基准线之间进行线性插值，形成一段螺旋刀轨，将所有的螺旋刀轨连成一条螺旋轨迹，并对螺旋轨迹进行曲线拟合，以形成几何二阶连续的B样条曲线刀轨。

步骤（2）具体为，首先在加工区域的三角网格上建立刀路轨迹计算模型如下：

min:

$F\left(\mathrm{\φ}\right)=\underset{(\mathrm{IJK}\∈\mathrm{\Π})}{\mathrm{\Σ}}\frac{1}{2A_{\left(\mathrm{IJK}\right)}}[{S_{\mathrm{IK}}}^2{({\mathrm{\φ}}_J-{\mathrm{\φ}}_I)}^2{S_{\mathrm{IJ}}}^2{({\mathrm{\φ}}_K-{\mathrm{\φ}}_I)}^2-2(S_{\mathrm{IK}}\·S_{\mathrm{IJ}})({\mathrm{\φ}}_K-{\mathrm{\φ}}_I)\·({\mathrm{\φ}}_J-{\mathrm{\φ}}_I)]$

s.t. $\frac{1}{{h_s}^2}-\frac{1}{{A_{\left(\mathrm{IJK}\right)}}^2}({S_{\mathrm{IK}}}^2{({\mathrm{\φ}}_J-{\mathrm{\φ}}_I)}^2+{S_{\mathrm{IJ}}}^2{({\mathrm{\φ}}_K-{\mathrm{\φ}}_I)}^2-2(S_{\mathrm{IK}}\·S_{\mathrm{IJ}})({\mathrm{\φ}}_K-{\mathrm{\φ}}_I)\·({\mathrm{\φ}}_J-{\mathrm{\φ}}_I))\≤0$ (IJK∈П)

φ_I＝0        在边界Г₁上

φ_I＝const＝n 在边界Г₂上

然后，利用序列二次规划方法，求解该刀路轨迹计算模型，得到势能场在每个节点上的近似计算值φ_I，将各节点的势能值作为高度值绘制出来，得到势能场分布。

通过本发明所构思的以上技术方案，与现有技术相比，本发明具有以下的有益效果：

1、将刀具路径光滑和行距控制优化纳入统一且具开放性的数学框架：由于采用了步骤（1）和步骤（2）基于能量泛函的约束优化模型以及有限元近似求解算法，计算出的加工区域上的势能场能用来产生光滑和行距优化的刀具轨迹轮廓线基准。该数学框架是开放性的，指的是其他加工需求，比如等残余高度、加工时间最短等等，可以以新的约束泛函的形式，纳入该求解框架进行统一优化求解；

2、产生只有一次进退刀的螺旋轨迹曲线：由于采用了步骤（3）的等势线投影和步骤（4）、（5）的轮廓线之间线性插值（Morphing）方法，能在具有任意形状的加工区域内产生只有一次进退刀的螺旋刀具轨迹；

3、产生适用于高速加工的光滑刀具轨迹：由于采用了步骤（5）的样条拟合方法，能产生出具有几何二阶（C2）连续的光滑刀具轨迹。

附图说明

图1是本发明基于势能场及能量泛函优化的刀具轨迹路径产生方法的流程图。

图2（a）示出加工区域的边界形状以及三角网格划分。

图2（b）示出加工区域内一连续变化的势能场示意图。

图2（c）示出势能场的等势线向加工区域的投影所形成的刀具轨迹基准轮廓线。

图3示出针对方钢叶片加工的势能场优化结果。

图4（a）示出根据势能场得到的势能场分布。

图4（b）示出根据势能场得到的等势轮廓线近似轨迹。

图5示出根据等势轮廓线得到的刀具轨迹轮廓基准。

图6示出通过在相邻轮廓线之间进行线性插值得到的螺旋轨迹。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。

本发明将加工区域视为有界曲面域，并用离散的三角网格来近似，如图2（a）所示，其网格区域内每个节点都具有确定的势能值，从而形成一标量势能场。若为加工区域边界赋予等常势能值，而不同的边界具有不同等常势能值，由于该标量场为连续场，则在整个加工区域内形成一连续变化的势能曲面，如图2（b）所示，其等势轮廓线代表了具有相同势能值的点的集合，场中一点的梯度向量代表了该点在势能曲面上“坡度”最陡的方向以及“坡度”到底有多“陡”。理想的刀具轨迹应该提供从外边界轮廓到内边界轮廓的光滑以及均匀过渡，可由设计良好（光滑且坡度均匀）的势能曲面上的等势轮廓线向加工曲面区域投影形成。将优化分布的等势轮廓线作为偏置轮廓线基准，在此基础上形成轮廓偏置或者螺旋线形状的刀轨。图2（c）示出等势轮廓线向加工曲面区域的投影，其中的行距未做优化控制。为了形成光滑和行距尽可能均匀的等势轮廓线，在势能场上构造能代表拉伸能量的泛函表达式，并利用其梯度场构造行距控制约束条件，将寻求优化分布的等势轮廓线问题转化成带约束的泛函极值问题。

前述带约束的泛函极值问题无法取得解析解，必须采用数值方法（有限差分法或者有限元方法）近似求解。有限元网格提供了对任意求解区域形状更好的适应性，所以本发明采用有限元方法，将泛函极值问题转化成多元非线性函数的约束优化问题，然后利用约束非线性优化数值方法求出每个有限元节点上的近似势能值，并以此构造等势轮廓线。由于优化后的势能场趋于光滑，相邻的等势轮廓线形状是渐进变化（Morphing）的，基于此种等势轮廓线并采用线性插值建立的螺旋轨迹行进平稳，不易产生相交和重叠。有限元数值计算是一种近似方法，基于有限元网格产生的轮廓线虽然提供了很好的刀轨分布趋势，但属于分段线性的，只有G0连续性。在此基础上可以采用曲线拟合的方式，产生光滑的刀轨。

本发明不仅能应用于产生行距优化的光滑刀轨，同时通过适当修改边界条件，也能将等距偏置、行切、环切等刀轨产生方法纳入统一的数学解决框架中。以本发明为基础，不仅能产生螺旋曲线刀轨，还能产生等距偏置的行切或环切刀轨。

同时本发明所建立数学求解框架是开放性的，意味着除了行距控制约束条件外，与其他数控加工需求相关的约束条件也可通过适当的数学建模，并纳入到本发明建立的统一优化框架中。

如图1所示，本发明基于势能场及能量泛函优化的刀具轨迹路径产生方法包括以下步骤：

（1）为多轴刀具建立基于能量泛函优化的刀路轨迹计算模型，包括能量泛函、边界条件和行宽约束条件，并采用有限元方法获取刀路轨迹优化计算模型基于三角形离散网格的数值近似计算公式，具体包括以下子步骤

（a）确定多轴刀具的加工区域及其外边界轮廓和内边界轮廓，并用三角形网格对整个加工区域进行网格划分，以构建势能场Ф；

如图2所示，图2（a）示出加工区域及其三角网格划分，共具有k个节点。图中Ω表示多轴刀具的加工区域（不包括边界轮廓）上的所有节点集合，Г₁表示加工区域外边界上的节点集合，Г₂表示加工区域内边界（即岛屿轮廓）上的节点集合；Ф_i(i＝1,2,...,k)表示加工区域中任意节点i所具有的势能值，整个加工区域形成一势能场Ф。

（b）为势能场Ф建立能量泛函F(φ)；

具体而言，本发明采用能量泛函优化的方法，求出一个类似于图2（b）中所示的势能场分布，其中外边界Г₁上Ф＝0，内边界Г₂上Ф＝n，n为由加工行距所决定的等势轮廓线的圈数；利用数值方法求解该优化问题后得到势能场Ф的值分布，然后求与Ф=1，2，…n对应的等势轮廓线，如图2（c）所示，作为产生刀具轨迹的基准轮廓线，并依此形成螺旋轨迹刀轨。

所采用的能量泛函目标函数为：

$\min :F\left(\mathrm{\Φ}\right)={\∫}_{{\mathrm{\Γ}}_1{\mathrm{\Γ}}_2\∪\mathrm{\Ω}}{(\▿\mathrm{\Φ})}^2\mathrm{dxdy}$

其中：min表示求最小值；F(Ф)表示以势能场Ф函数为变量的能量泛函，该泛函的物理含义类似于图2（b）所示的势能曲面中的拉伸能；表示势能场Ф中的梯度场（也称矢量场），其计算式为dx和dy表示平面坐标的微分符号。

设Ω对应于二维参数曲面S（u¹，u²），则该能量泛函计算式为：

$F\left(\mathrm{\φ}\right)\≡\underset{D}{\∫}{\mathrm{du}}^1{\mathrm{du}}^2\sqrt{G\left(u\right)}(\frac{1}{2}{g}^{\mathrm{\α\β}}\left(u\right)\frac{\mathrm{d\Φ}\left(u\right)}{{\mathrm{du}}^{\mathrm{\α}}}\·\frac{\mathrm{d\Φ}\left(u\right)}{{\mathrm{du}}^{\mathrm{\β}}})---\left(1\right)$

其中，α,β为按爱因斯坦求和约定的标号，取值分别为1和2。u¹,u²为曲面参数，u＝(u¹,u²)^T。当标号α或者β在方程的一个单独项目内重复出现时，表示求和，如

$g^{\mathrm{\α\β}}\left(u\right)\frac{\mathrm{d\Φ}\left(u\right)}{{\mathrm{du}}^{\mathrm{\α}}}\·\frac{\mathrm{d\Φ}\left(u\right)}{{\mathrm{du}}^{\mathrm{\β}}}=\underset{\mathrm{\α}=\mathrm{1,2}}{\mathrm{\Σ}}\underset{\mathrm{\β}=\mathrm{1,2}}{\mathrm{\Σ}}{g}^{\mathrm{\α\β}}\left(u\right)\frac{\mathrm{d\Φ}\left(u\right)}{{\mathrm{du}}^{\mathrm{\α}}}\·\frac{\mathrm{d\Φ}\left(u\right)}{{\mathrm{du}}^{\mathrm{\β}}};$

而单独出现时，表示一个向量，其元素按照标号的取值选取，如表示向量G^αβ为曲面度量张量（曲面高斯第一基本形）的矩阵， $G_{\mathrm{\α\β}}=\left[\begin{array}{cc}{G}_{11}& G_{12}\\ G_{21}& G_{22}\end{array}\right],$ G(u)＝det(G_αβ)，为曲面度量张量G_αβ的行列式计算值，det（）表示求行列式。g^αβ(u)为曲面度量张量的逆。

（c）获取能量泛函F(φ)的有限元数值近似计算公式；

图3示出加工区域中的一个三角形网格在参数空间和其对应的物理空间中的定义，用一个三元数（I，J，K）表示加工区域中的一个三角形，三角形网格的顶点为u_I，u_J，u_K，I，J，K分别表示三个顶点的节点序号，u_I，u_J，u_K对应物理空间中的三点S_I=S（u_I），S_J＝S（u_J），S_K=S（u_K），S表示从参数空间到物理空间的坐标转换。物理空间中三角形的三边可用矢量S_IJ，S_JK，S_KI表示，分别计算为S_IJ=S_J-S_I，S_JK=S_K-S_J，S_KI=S_I-S_K。

如图3所示，在每个三角形网格的参数空间建立以u_I为原点的非正交坐标系，其坐标向量分别为du¹＝u_J-u_I，du²=u_K-u_I。利用非正交坐标系以及微分几何知识，该坐标系下的曲面度量张量 $G_{\mathrm{\α\β}}=\left[\begin{array}{cc}{G}_{11}& G_{12}\\ G_{21}& G_{22}\end{array}\right]$ 的四个元素可计算为：

G₁₁＝S_IJ ²

G₂₂＝S_IK ²

G₁₂＝G₂₁＝S_IJ·S_IK

度量张量G_αβ的逆g^αβ可计算为：

$g^{\mathrm{\α\β}}=\left[\begin{array}{cc}{g}^{11}& g^{12}\\ g^{21}& g^{22}\end{array}\right]$

其中：

$g^{11}=\frac{{S_{\mathrm{IK}}}^2}{{A_{\left(\mathrm{IJK}\right)}}^2}$

$g^{22}=\frac{{S_{\mathrm{IJ}}}^2}{{A_{\left(\mathrm{IJK}\right)}}^2}$

$g^{12}=g^{21}=-\frac{S_{\mathrm{IJ}}\·S_{\mathrm{IK}}}{{A_{\left(\mathrm{IJK}\right)}}^2}$

其中 $A_{\left(\mathrm{IJK}\right)}=\sqrt{\det \left(G_{\mathrm{\α\β}}\right)}=\sqrt{{S_{\mathrm{IK}}}^2{S_{\mathrm{IJ}}}^2-{(S_{\mathrm{IK}}\·S_{\mathrm{IJ}})}^2},$ det（）表示求度量张量G_αβ的行列式。

通过以上等式容易得出式（1）中各项在有限元网格（I,J,K）上的计算式：

$\sqrt{G\left(u\right)}=A_{\left(\mathrm{IJK}\right)}$

$g^{\mathrm{\α\β}}\left(u\right)=\left[\begin{array}{cc}\frac{{S_{\mathrm{IK}}}^2}{{A_{\left(\mathrm{IJK}\right)}}^2}& \text{-}\frac{S_{\mathrm{IJ}}\·S_{\mathrm{IK}}}{{A_{\left(\mathrm{IJK}\right)}}^2}\\ -\frac{S_{\mathrm{IJ}}\·S_{\mathrm{IK}}}{{A_{\left(\mathrm{IJK}\right)}}^2}& \frac{{S_{\mathrm{IJ}}}^2}{{A_{\left(\mathrm{IJK}\right)}}^2}\end{array}\right]$

设势能场Ф在节点I上的值为φ_I＝φ(u_I)，φ_K＝φ(u_K)代表势能场Ф在节点K上的值，φ_J＝φ(u_J)代表势能场Ф在节点J上的值；

该势能场的各种导数在三角网格上可近似计算为:

$\frac{\mathrm{d\φ}}{{\mathrm{du}}^1}={\mathrm{\φ}}_J-{\mathrm{\φ}}_I$

$\frac{\mathrm{d\φ}}{{\mathrm{du}}^2}={\mathrm{\φ}}_K-{\mathrm{\φ}}_I$

其中u^α和u^β满足爱因斯坦求和约定规则，表示向量（u¹,u²）。

将前面列出的有限元网格上的G(u)，g^αβ(u),以及的计算式代入式（1），可得到与能量泛函式（1）等价的有限元近似计算公式为：

$F\left(\mathrm{\φ}\right)=\underset{(\mathrm{IJK}\∈\mathrm{\Π})}{\mathrm{\Σ}}\frac{1}{2A_{\left(\mathrm{IJK}\right)}}[{s_{\mathrm{IK}}}^2{({\mathrm{\φ}}_J-{\mathrm{\φ}}_I)}^2+{s_{\mathrm{IJ}}}^2{({\mathrm{\φ}}_K-{\mathrm{\φ}}_I)}^2-2(s_{\mathrm{IK}}\·s_{\mathrm{IJ}})({\mathrm{\φ}}_K-{\mathrm{\φ}}_I)\·({\mathrm{\φ}}_J-{\mathrm{\φ}}_I)]---\left(2\right)$

其中П表示全体三角形网格的集合。

从式（2）可以看出，能量泛函已经变成了一个以各节点I的势能值φ_I为自变量，并包含三角形网格几何参数信息的多元目标函数。

（d）利用行距条件，建立行距控制约束泛函以及其数值近似计算公式；

设等势轮廓线j上的势能值为φ_j，与之无限接近的等势轮廓线（j+1）的势能值为φ_j+Δφ。则轮廓线j和j+1之间的势能差为Δφ；同理，轮廓线j和j+1向参数曲面（加工区域Ω）投影而形成的两条平面轮廓线之间的距离用Δx表示。则势能场Ф的梯度可计算为

$\▿\mathrm{\φ}=\frac{\mathrm{d\φ}}{\mathrm{dx}}\≈\frac{\mathrm{\Δ\φ}}{\mathrm{\Δx}}$

其中dφ表示等势轮廓线之间的势能值微分，dx表示等势轮廓线在加工区域Ω的投影轮廓线之间的位置微分。人为设定两相邻刀轨轮廓基准线之间的Δφ＝1，则其对应的加工区域投影轮廓之间的距离应满足：

$\left|\mathrm{\Δx}\right|=\frac{\left|\mathrm{\Δ\φ}\right|}{|\▿\mathrm{\φ}|}=\frac{1}{|\▿\mathrm{\φ}|}\≤h_s$

其中h_s表示加工行距，通常取刀具直径的70~80%。

从而所采用的行距控制约束泛函为：

$|\▿\mathrm{\φ}|\≥\frac{1}{h_s}\⇔\frac{1}{{h_s}^2}-{(\▿\mathrm{\φ})}^2\≤o$ 在区域Ω上

由于：

$\▿\mathrm{\φ}={(\frac{\mathrm{d\φ}}{{\mathrm{du}}^1},\frac{\mathrm{d\φ}}{{\mathrm{du}}^2})}^T$

${(\▿\mathrm{\φ})}^2=\frac{\mathrm{d\φ}}{{\mathrm{du}}^{\mathrm{\α}}}\·\frac{\mathrm{d\φ}}{{\mathrm{du}}^{\mathrm{\β}}}$

则行距控制约束泛函对应的参数空间计算式为

$\frac{1}{{h_s}^2}-g^{\mathrm{\α\β}}\left(u\right)\frac{\mathrm{d\φ}\left(u\right)}{{\mathrm{du}}^{\mathrm{\α}}}\·\frac{\mathrm{d\φ}\left(u\right)}{{\mathrm{du}}^{\mathrm{\β}}}\≤o---\left(3\right)$

利用类似于针对能量泛函的有限元数值近似方法，将前面列出的有限元网格上的G(u)，g^αβ(u),以及的计算式代入式（3），可以获取行距控制约束泛函的有限元近似计算公式为：

$\frac{1}{{h_s}^2}-\frac{1}{{A_{\left(\mathrm{IJK}\right)}}^2}({S_{\mathrm{IK}}}^2{({\mathrm{\φ}}_J-{\mathrm{\φ}}_I)}^2+{S_{\mathrm{IJ}}}^2{({\mathrm{\φ}}_K-{\mathrm{\φ}}_I)}^2-2(S_{\mathrm{IK}}\·S_{\mathrm{IJ}})({\mathrm{\φ}}_K-{\mathrm{\φ}}_I)\·({\mathrm{\φ}}_J-{\mathrm{\φ}}_I))\≤0$

(IJK∈П)       （4）

φ_I＝0           在边界Г₁上

φ_I＝const＝n    在边界Г₂上

从式（4）可以看出，行距控制约束泛函已经变成了一系列以三角网格节点I的势能值φ_I为自变量，并包含该三角网格几何参数信息的多元不等式函数。

（e）建立基于能量泛函优化的刀路轨迹计算边界条件；

该优化问题需要设定边界条件才可获得具体解。不同的边界条件组合可以产生不同的优化结果。对于图2所示的带岛屿的平面区域螺旋加工采用的边界条件为

φ＝0    在边界Г₁上

φ＝n    在边界Г₂上

（2）根据刀路轨迹计算模型得到势能场Ф在每个节点I上的近似计算值φ_I；

具体而言，如前所述，在加工区域的三角网格上为了求得光滑并满足行宽约束的刀具轨迹基准线而建立的刀路轨迹计算模型如下：

min:

$F\left(\mathrm{\φ}\right)=\underset{(\mathrm{IJK}\∈\mathrm{\Π})}{\mathrm{\Σ}}\frac{1}{2A_{\left(\mathrm{IJK}\right)}}[{S_{\mathrm{IK}}}^2{({\mathrm{\φ}}_J-{\mathrm{\φ}}_I)}^2{S_{\mathrm{IJ}}}^2{({\mathrm{\φ}}_K-{\mathrm{\φ}}_I)}^2-2(S_{\mathrm{IK}}\·S_{\mathrm{IJ}})({\mathrm{\φ}}_K-{\mathrm{\φ}}_I)\·({\mathrm{\φ}}_J-{\mathrm{\φ}}_I)]$

s.t. $\frac{1}{{h_s}^2}-\frac{1}{{A_{\left(\mathrm{IJK}\right)}}^2}({S_{\mathrm{IK}}}^2{({\mathrm{\φ}}_J-{\mathrm{\φ}}_I)}^2+{S_{\mathrm{IJ}}}^2{({\mathrm{\φ}}_K-{\mathrm{\φ}}_I)}^2-2(S_{\mathrm{IK}}\·S_{\mathrm{IJ}})({\mathrm{\φ}}_K-{\mathrm{\φ}}_I)\·({\mathrm{\φ}}_J-{\mathrm{\φ}}_I))\≤0$ (IJK∈П)

φ_I＝0         在边界Г₁上

φ_I＝const＝n  在边界Г₂上

然后，本步骤利用序列二次规划（Sequential Quadratic Programming，简称SQP）方法，求解该刀路轨迹计算模型，得到势能场在每个节点上的近似计算值φ_I。将各节点的势能值作为高度值绘制出来，得到的势能场分布如图4（a）所示，其形状类似一个火山口。火山口底部为外边界，其势能值为0；火山口顶部为内边界，其势能值为6。该非线性约束优化问题的求解应不局限于SQP方法，也可采用其他适用于非线性约束优化问题的求解方法。

（3）按Ф=0，1，…，6，截取火山口的等高线，并向加工区域Ω投影，形成等势轮廓线近似轨迹；如图4（b）所示，由于采用有限元近似计算技术，得到的近似轨迹是分段线性的；

（4）对得到的等势轮廓线近似轨迹进行平滑和采样，得到从外向内偏置的刀具轨迹轮廓基准线；如图5所示，相邻轮廓基准线之间的距离满足加工行距条件，且其变化是光滑的；

（5）在刀具轨迹轮廓基准线之间进行线性插值，形成一段螺旋刀轨，然后将所有的螺旋刀轨连成一条螺旋轨迹，并对螺旋轨迹进行曲线拟合，以形成几何二阶连续（C2）的B样条曲线刀轨，如图6所示。

本领域的技术人员容易理解，以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。

Claims (1)

1.一种基于势能场及能量泛函优化的刀具轨迹路径产生方法，其特征在于，包括以下步骤：

(1)为多轴刀具建立基于能量泛函优化的刀路轨迹计算模型，包括能量泛函、边界条件和行宽约束条件，并采用有限元方法获取刀路轨迹优化计算模型基于三角形离散网格的数值近似计算公式，具体包括以下子步骤：

(a)确定多轴刀具的加工区域及其外边界轮廓和内边界轮廓，并用三角形网格对整个加工区域进行网格划分，以构建势能场Φ，其中加工区域及其三角网格划分共具有k个节点，多轴刀具的加工区域上的所有节点集合为Ω，加工区域外边界上的节点集合为Г₁，加工区域内边界上的节点集合为Г₂，Φ_i(i＝1,2,...,k)表示加工区域中任意节点i所具有的势能值，整个加工区域形成一势能场Φ；

(b)为势能场Φ建立能量泛函F(φ)，具体公式为：

$F\left(\mathrm{\φ}\right)\≡\underset{D}{\∫}{\mathrm{du}}^1{\mathrm{du}}^2\sqrt{G\left(u\right)}(\frac{1}{2}{g}^{\mathrm{\α\β}}\left(u\right)\frac{\mathrm{d\Φ}\left(u\right)}{{\mathrm{du}}^{\mathrm{\α}}}\·\frac{\mathrm{d\Φ}\left(u\right)}{{\mathrm{du}}^{\mathrm{\β}}})$

其中α,β为按爱因斯坦求和约定的标号，取值分别为1和2，u¹,u²为曲面参数， $u={(u^1,u^2)}^T,g^{\mathrm{\α\β}}\left(u\right)\frac{\mathrm{d\Φ}\left(u\right)}{{\mathrm{du}}^{\mathrm{\α}}}\·\frac{\mathrm{d\Φ}\left(u\right)}{{\mathrm{du}}^{\mathrm{\β}}}=\underset{\mathrm{\α}=\mathrm{1,2}}{\mathrm{\Σ}}\underset{\mathrm{\β}=1.2}{\mathrm{\Σ}}{g}^{\mathrm{\α\β}}\left(u\right)\frac{\mathrm{d\Φ}\left(u\right)}{{\mathrm{du}}^{\mathrm{\α}}}\·\frac{\mathrm{d\Φ}\left(u\right)}{{\mathrm{du}}^{\mathrm{\β}}},\frac{\mathrm{d\φ}}{{\mathrm{du}}^{\mathrm{\α}}}$ 表示向量G^αβ为曲面度量张量的矩阵，且 $G_{\mathrm{\α\β}}=\left[\begin{array}{cc}{G}_{11}& G_{12}\\ G_{21}& G_{22}\end{array}\right],$ G(u)＝det(G_αβ)表示曲面度量张量G_αβ的行列式计算值，det()表示求行列式，g^αβ(u)为曲面度量张量的逆，u^α和u^β满足爱因斯坦求和约定规则，表示曲面参数(u¹,u²)；

(c)获取能量泛函F(φ)的有限元数值近似计算公式，采用以下公式

$F\left(\mathrm{\φ}\right)=\underset{(\mathrm{IJK}\∈\mathrm{\Π})}{\mathrm{\Σ}}\frac{1}{{2A}_{\left(\mathrm{IJK}\right)}}[{S_{\mathrm{IK}}}^2{({\mathrm{\φ}}_J-{\mathrm{\φ}}_I)}^2+{S_{\mathrm{IJ}}}^2{({\mathrm{\φ}}_K-{\mathrm{\φ}}_I)}^2-2(S_{\mathrm{IK}}\·S_{\mathrm{IJ}})({\mathrm{\φ}}_K-{\mathrm{\φ}}_I)\·({\mathrm{\φ}}_J-{\mathrm{\φ}}_I)]$

其中(I，J，K)表示加工区域中的一个三角形，三角形网格的顶点为u_I，u_J，u_K，I，J，K分别表示三个顶点的节点序号，u_I，u_J，u_K对应物理空间中的三点S_I＝S(u_I)，S_J＝S(u_J)，S_K＝S(u_K)，S表示从参数空间到物理空间的坐标转换，S_IJ，S_JK，S_KI表示物理空间中三角形的三边可用矢量，分别为S_IJ＝S_J-S_I，S_JK＝S_K-S_J，S_KI＝S_I-S_K， $A_{\left(\mathrm{IJK}\right)}=\sqrt{\det \left(G_{\mathrm{\α\β}}\right)}=\sqrt{{S_{\mathrm{IK}}}^2{S_{\mathrm{IJ}}}^2-{(S_{\mathrm{IK}}\·S_{\mathrm{IJ}})}^2},$ det()表示求度量张量G_αβ的行列式，φ_I＝φ(u_I)代表势能场Φ在节点I上的值，φ_K＝φ(u_K)代表势能场Φ在节点K上的值,φ_J＝φ(u_J)代表势能场Φ在节点J上的值,Π表示全体三角形网格的集合；

(d)利用行距条件，建立行距控制约束泛函以及其数值近似计算公式，采用以下公式：

$\frac{1}{{h_s}^2}-\frac{1}{{A_{\left(\mathrm{IJK}\right)}}^2}({{S_{\mathrm{IK}}}^2({\mathrm{\φ}}_J-{\mathrm{\φ}}_I)}^2+{S_{\mathrm{IJ}}}^2{({\mathrm{\φ}}_K-{\mathrm{\φ}}_I)}^2-2(S_{\mathrm{IK}}\·S_{\mathrm{IJ}})({\mathrm{\φ}}_k-{\mathrm{\φ}}_I)\·({\mathrm{\φ}}_J-{\mathrm{\φ}}_I))\≤0$

(IJK∈Π)

φ_I＝0            在边界Γ₁上

φ_I＝const＝n     在边界Γ₂上

其中h_s表示加工行距；

(e)建立基于能量泛函优化的刀路轨迹计算边界条件；

(2)根据刀路轨迹计算模型得到势能场Φ在每个节点I上的近似计算值φ_I；具体为，首先在加工区域的三角网格上建立刀路轨迹计算模型如下：

min:

$F\left(\mathrm{\φ}\right)=\underset{(\mathrm{IJK}\∈\mathrm{\Π})}{\mathrm{\Σ}}\frac{1}{{2A}_{\left(\mathrm{IJK}\right)}}[{S_{\mathrm{IK}}}^2{({\mathrm{\φ}}_J-{\mathrm{\φ}}_I)}^2+{S_{\mathrm{IJ}}}^2{({\mathrm{\φ}}_K-{\mathrm{\φ}}_I)}^2-2(S_{\mathrm{IK}}\·S_{\mathrm{IJ}})({\mathrm{\φ}}_K-{\mathrm{\φ}}_I)\·({\mathrm{\φ}}_J-{\mathrm{\φ}}_I)]$

s.t. $\frac{1}{{h_s}^2}-\frac{1}{{A_{\left(\mathrm{IJK}\right)}}^2}({{S_{\mathrm{IK}}}^2({\mathrm{\φ}}_J-{\mathrm{\φ}}_I)}^2+{S_{\mathrm{IJ}}}^2{({\mathrm{\φ}}_K-{\mathrm{\φ}}_I)}^2-2(S_{\mathrm{IK}}\·S_{\mathrm{IJ}})({\mathrm{\φ}}_k-{\mathrm{\φ}}_I)\·({\mathrm{\φ}}_J-{\mathrm{\φ}}_I))\≤0(\mathrm{IJK}\∈\mathrm{\Π})$

φ_I＝0            在边界Γ₁上

φ_I＝const＝n     在边界Γ₂上

然后，利用序列二次规划方法，求解该刀路轨迹计算模型，得到势能场在每个节点上的近似计算值φ_I，将各节点的势能值作为高度值绘制出来，得到势能场分布；

(3)按Φ＝0，1，…，6，截取势能场的等高线，并向加工区域Ω投影，形成等势轮廓线近似轨迹；

(4)对得到的等势轮廓线近似轨迹进行平滑和采样，得到从外向内偏置的刀具轨迹轮廓基准线；

(5)在刀具轨迹轮廓基准线之间进行线性插值，形成一段螺旋刀轨，将所有的螺旋刀轨连成一条螺旋轨迹，并对螺旋轨迹进行曲线拟合，以形成几何二阶连续的B样条曲线刀轨。