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

摘要

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

说明书

技术领域

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

背景技术

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

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

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

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

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

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

发明内容

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

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

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

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

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

$F\left(\phi \right)\equiv \underset{D}{\int }{\mathrm{du}}^1{\mathrm{du}}^2\sqrt{G\left(u\right)}(\frac{1}{2}{g}^{\mathrm{\alpha \beta }}\left(u\right)\frac{\mathrm{d\Phi }\left(u\right)}{{\mathrm{du}}^{\alpha }}·\frac{\mathrm{d\Phi }\left(u\right)}{{\mathrm{du}}^{\beta }})$

其中α,β为按爱因斯坦求和约定的标号，取值分别为1和2，u¹,u²为曲 面参数，u＝(u¹,u²)^T，$g^{\mathrm{\alpha \beta }}\left(u\right)\frac{\mathrm{d\Phi }\left(u\right)}{{\mathrm{du}}^{\alpha }}·\frac{\mathrm{d\Phi }\left(u\right)}{{\mathrm{du}}^{\beta }}=\underset{\alpha =\mathrm{1,2}}{\Sigma }\underset{\beta =\mathrm{1,2}}{\Sigma }{g}^{\mathrm{\alpha \beta }}\left(u\right)\frac{\mathrm{d\Phi }\left(u\right)}{{\mathrm{du}}^{\alpha }}·\frac{\mathrm{d\Phi }\left(u\right)}{{\mathrm{du}}^{\beta }},$$\frac{\mathrm{d\phi }}{{\mathrm{du}}^{\alpha }}$表 示向量G^αβ为曲面度量张量的矩阵，且$G_{\mathrm{\alpha \beta }}=\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(\phi \right)=\underset{(\mathrm{IJK}\in \Pi )}{\Sigma }\frac{1}{2A_{\left(\mathrm{IJK}\right)}}[{S_{\mathrm{IK}}}^2{({\phi }_J-{\phi }_I)}^2{S_{\mathrm{IJ}}}^2{({\phi }_K-{\phi }_I)}^2-2(S_{\mathrm{IK}}·S_{\mathrm{IJ}})({\phi }_K-{\phi }_I)·({\phi }_J-{\phi }_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{\alpha \beta }}\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{({\phi }_J-{\phi }_I)}^2+{S_{\mathrm{IJ}}}^2{({\phi }_K-{\phi }_I)}^2-2(S_{\mathrm{IK}}·S_{\mathrm{IJ}})({\phi }_K-{\phi }_I)·({\phi }_J-{\phi }_I))\le 0$

(IJK∈П)

φ_I＝0           在边界Г₁上

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

其中h_s表示加工行距；

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

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

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

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

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

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

min:

$F\left(\phi \right)=\underset{(\mathrm{IJK}\in \Pi )}{\Sigma }\frac{1}{2A_{\left(\mathrm{IJK}\right)}}[{S_{\mathrm{IK}}}^2{({\phi }_J-{\phi }_I)}^2{S_{\mathrm{IJ}}}^2{({\phi }_K-{\phi }_I)}^2-2(S_{\mathrm{IK}}·S_{\mathrm{IJ}})({\phi }_K-{\phi }_I)·({\phi }_J-{\phi }_I)]$

s.t.$\frac{1}{{h_s}^2}-\frac{1}{{A_{\left(\mathrm{IJK}\right)}}^2}({S_{\mathrm{IK}}}^2{({\phi }_J-{\phi }_I)}^2+{S_{\mathrm{IJ}}}^2{({\phi }_K-{\phi }_I)}^2-2(S_{\mathrm{IK}}·S_{\mathrm{IJ}})({\phi }_K-{\phi }_I)·({\phi }_J-{\phi }_I))\le 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(\Phi \right)={\int }_{{\Gamma }_1{\Gamma }_2\cup \Omega }{(▿\Phi )}^2\mathrm{dxdy}$

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

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

$F\left(\phi \right)\equiv \underset{D}{\int }{\mathrm{du}}^1{\mathrm{du}}^2\sqrt{G\left(u\right)}(\frac{1}{2}{g}^{\mathrm{\alpha \beta }}\left(u\right)\frac{\mathrm{d\Phi }\left(u\right)}{{\mathrm{du}}^{\alpha }}·\frac{\mathrm{d\Phi }\left(u\right)}{{\mathrm{du}}^{\beta }})---\left(1\right)$

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

$g^{\mathrm{\alpha \beta }}\left(u\right)\frac{\mathrm{d\Phi }\left(u\right)}{{\mathrm{du}}^{\alpha }}·\frac{\mathrm{d\Phi }\left(u\right)}{{\mathrm{du}}^{\beta }}=\underset{\alpha =\mathrm{1,2}}{\Sigma }\underset{\beta =\mathrm{1,2}}{\Sigma }{g}^{\mathrm{\alpha \beta }}\left(u\right)\frac{\mathrm{d\Phi }\left(u\right)}{{\mathrm{du}}^{\alpha }}·\frac{\mathrm{d\Phi }\left(u\right)}{{\mathrm{du}}^{\beta }};$

而单独出现时，表示一个向量，其元素按照标号的取值选取，如表示 向量G^αβ为曲面度量张量（曲面高斯第一基本形）的矩阵， $G_{\mathrm{\alpha \beta }}=\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{\alpha \beta }}=\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{\alpha \beta }}=\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{\alpha \beta }}\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{\alpha \beta }}\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\phi }}{{\mathrm{du}}^1}={\phi }_J-{\phi }_I$

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

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

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

$F\left(\phi \right)=\underset{(\mathrm{IJK}\in \Pi )}{\Sigma }\frac{1}{2A_{\left(\mathrm{IJK}\right)}}[{s_{\mathrm{IK}}}^2{({\phi }_J-{\phi }_I)}^2+{s_{\mathrm{IJ}}}^2{({\phi }_K-{\phi }_I)}^2-2(s_{\mathrm{IK}}·s_{\mathrm{IJ}})({\phi }_K-{\phi }_I)·({\phi }_J-{\phi }_I)]---\left(2\right)$

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

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

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

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

$▿\phi =\frac{\mathrm{d\phi }}{\mathrm{dx}}\approx \frac{\mathrm{\Delta \phi }}{\mathrm{\Delta x}}$

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

$\left|\mathrm{\Delta x}\right|=\frac{\left|\mathrm{\Delta \phi }\right|}{|▿\phi |}=\frac{1}{|▿\phi |}\le h_s$

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

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

$|▿\phi |\ge \frac{1}{h_s}\iff \frac{1}{{h_s}^2}-{(▿\phi )}^2\le o$在区域Ω上

由于：

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

${(▿\phi )}^2=\frac{\mathrm{d\phi }}{{\mathrm{du}}^{\alpha }}·\frac{\mathrm{d\phi }}{{\mathrm{du}}^{\beta }}$

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

$\frac{1}{{h_s}^2}-g^{\mathrm{\alpha \beta }}\left(u\right)\frac{\mathrm{d\phi }\left(u\right)}{{\mathrm{du}}^{\alpha }}·\frac{\mathrm{d\phi }\left(u\right)}{{\mathrm{du}}^{\beta }}\le 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{({\phi }_J-{\phi }_I)}^2+{S_{\mathrm{IJ}}}^2{({\phi }_K-{\phi }_I)}^2-2(S_{\mathrm{IK}}·S_{\mathrm{IJ}})({\phi }_K-{\phi }_I)·({\phi }_J-{\phi }_I))\le 0$

(IJK∈П)       （4）

φ_I＝0           在边界Г₁上

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

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

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

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

φ＝0    在边界Г₁上

φ＝n    在边界Г₂上

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

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

min:

$F\left(\phi \right)=\underset{(\mathrm{IJK}\in \Pi )}{\Sigma }\frac{1}{2A_{\left(\mathrm{IJK}\right)}}[{S_{\mathrm{IK}}}^2{({\phi }_J-{\phi }_I)}^2{S_{\mathrm{IJ}}}^2{({\phi }_K-{\phi }_I)}^2-2(S_{\mathrm{IK}}·S_{\mathrm{IJ}})({\phi }_K-{\phi }_I)·({\phi }_J-{\phi }_I)]$

s.t.$\frac{1}{{h_s}^2}-\frac{1}{{A_{\left(\mathrm{IJK}\right)}}^2}({S_{\mathrm{IK}}}^2{({\phi }_J-{\phi }_I)}^2+{S_{\mathrm{IJ}}}^2{({\phi }_K-{\phi }_I)}^2-2(S_{\mathrm{IK}}·S_{\mathrm{IJ}})({\phi }_K-{\phi }_I)·({\phi }_J-{\phi }_I))\le 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所示。

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