流线型加工路径

基本介绍

流线是表征流体运动情况的空间曲线，指在某一瞬间与流场中流体质点流速方向相切的光滑曲线，流线上各点的切线方向就是该点速度方向.

型腔加工中的刀具运动，类似于流体质点在平面流场中的流动，因此将平面流场中的流线做适当的行距控制后，即可作为型腔加工的刀具轨迹.

优势: 流线既不相交，也不突然转折，并且是光滑曲线，非常适合型腔的高速铣削.

流函数的数学表达

在平面流动中，不可压缩流体的连续性方程为
$$\nabla v=\frac{\partial v_x}{\partial x}+\frac{\partial v_y}{\partial y}=0.$$
此式也是微分表达式$-v_{y}dx+v_{x}dy$成为某一数量函数$\Psi$ 全微分的充要条件，即
$$\begin{array}{c}\{\begin{array}{l}\frac{\partial \Psi }{\partial y}=v_x,\\ \frac{\partial \Psi }{\partial x}=-v_y.\end{array}\end{array}$$
符合条件的数量函数$\Psi$ 被称为二维不可压缩流体平面运动的流函数.

基本性质

（1）平面有势流动的流函数为调和函数.

不可压缩流体的平面势流满足
$${\omega }_z=\frac{1}{2}(\frac{\partial v_y}{\partial x}-\frac{\partial v_x}{\partial y})=0,$$
则有
$$\Delta \Psi =\frac{{\partial }^2\Psi }{\partial x^2}+\frac{{\partial }^2\Psi }{\partial y^2}=0.$$
说明平面有势流动的流函数满足$Laplace$方程，是调和函数.

（2）平面流函数的等值线是流线.

$$d\Psi =\frac{\partial \Psi }{\partial x}dx+\frac{\partial \Psi }{\partial y}dy=(-v_y)dx+v_{x}dy,$$
在流函数的任意等值线$\Psi =c_i$上，$d\Psi =0$.代入上式可得
$$\frac{dx}{v_x}=\frac{dy}{v_y}.$$
该微分关系式与流线关于流线上每一点的切线与速度矢量重合的定义相符，由此可得流函数的等值线就是流线.

平面速度矢量场的重构

基本步骤

• 为加工区域中的关键点指定速度矢量，或在给定的边界条件下通过有限元方法获得加工区域离散网格点上的速度分布；
• 利用插值或逼近技术，构造该速度矢量场的流函数；
• 对平面流线微分方程进行数值积分，获得用于构造加工路径的流线场.

基于B样条基函数多项式的矢量场重构

为增强重构矢量场的整体平滑效果，提高对矢量场的局部调控能力，采用B样条基函数多项式来构造速度矢量场的流函数，实现矢量场的整体重建.

假定加工区域网格节点$v_s(x_s,y_s)$上的速度矢量分布为${\tau }_s=[{\tau }_s^u,{\tau }_s^v{]}^T(s=0,...,t)$,节点对应参数值为$(u_s,v_s)$,并得到节点矢量U,V.

设由B样条基函数表示的流函数
$$\Psi (u,v)=\sum _{i=0}^m\sum _{j=0}^n{d}_{i,j}{N}_{i,k}(u)N_{j,l}(v),$$
其梯度
$$\nabla \Psi (u,v)=[\frac{\partial \Psi (u,v)}{\partial u},\frac{\partial \Psi (u,v)}{\partial v}{]}^T,$$
由于流函数的等值线是流线，而梯度是等值线法线，可得
$${\tau }_s\cdot \nabla \Psi =0.$$
即
$$\begin{gathered} \begin{array}{c}\frac{\partial \Psi (u,v)}{\partial u}=-{\tau }_s^v, \\ \frac{\partial \Psi (u,v)}{\partial v}={\tau }_s^u.\end{array} \end{gathered}$$
由B样条曲面的偏导矢公式，在每一点处均有
$$AX=b,$$
其中
$$A=\left[\begin{array}{ccccccc}{a}_{0,0}^u& a_{0,1}^u& ...& a_{0,n}^u& a_{1,1}^u& ...& a_{m,n}^u\\ a_{0,0}^v& a_{0,1}^v& ...& a_{0,n}^v& a_{1,1}^v& ...& a_{m,n}^v\end{array}\right].$$

$$\begin{array}{c}{a}_{i,j}^u=\{\begin{array}{l}-k{N}_{j,l}(v_s)\frac{N_{1,k-1}(u_s)}{u_{k+1}-u_1},i=0,\\ k{N}_{j,l}(v_s)(\frac{N_{i,k-1}(u_s)}{u_{i+k}-u_i}-\frac{N_{i+1,k-1}(u_s)}{u_{i+k+1}-u_{i+1}}),i=1,...,m-1,\\ k{N}_{j,l}(v_s)\frac{N_{m,k-1}(u_s)}{u_{m+k}-u_m},i=m.\end{array}\end{array}$$

$$\begin{array}{c}{a}_{i,j}^v=\{\begin{array}{l}-l{N}_{i,k}(u_s)\frac{N_{1,l-1}(v_s)}{v_{l+1}-v_1},j=0,\\ l{N}_{i,k}(u_s)(\frac{N_{j,l-1}(v_s)}{v_{j+l}-v_j}-\frac{N_{j+1,l-1}(v_s)}{v_{j+l+1}-v_{j+1}}),j=1,...,n-1,\\ l{N}_{i,k}(u_s)\frac{N_{n,l-1}(v_s)}{v_{n+k}-v_n},j=n.\end{array}\end{array}$$

将所有采样点处的方程汇总
$$HX=B.$$
其中$H=[A_0,A_1,...,A_t],B=[b_0,b_1,...,b_t],X=[d_{0,0},d_{0,1},...,d_{0,n},d_{1,0},...,d_{m,n}{]}^T$.

利用最小二乘法可确定上述方程的解，进而得到流函数$\Psi (u,v)$的具体形式.

加工路径生成

平面流线的微分方程可改写为
$$\frac{dy}{dx}=\frac{v_x}{v_y}=f(x,y)$$
式中$(x,y)$为加工区域上的点，与流函数的参数$(u,v)$对应，对其进行求解时可采用较高求解精度的4阶$Runge-Kutta$方法：
$$\begin{array}{c}\{\begin{array}{l}{y}_{n+1}=y_n+\frac{h}{6}(k_1+2k_2+2k_3+k_4),\\ k_1=f(x_n,y_n),\\ k_2=f(x_n+\frac{1}{2}h,y_n+\frac{1}{2}h{k}_1),\\ k_3=f(x_n+\frac{1}{2}h,y_n+\frac{1}{2}h{k}_2),\\ k_4=f(x_n+h,y_n+h{k}_3).\end{array}\end{array}$$
这样，就可从加工区域任取一点$(x_0,y_0)$作为初始点生成一条初始加工路径，然后在可行走刀行距L下，计算初始加工路径采样点 { p i } \{p_i\} {pi​}在相邻路径的采样点 { p i 0 } \{p_i^0\} {pi0​}，根据采样点间最小增量$(\Delta x,\Delta y)$得到下一条流线的开始点，重复上述步骤，直至流线型加工路径充满整个加工区域.

代码

Remark:在流线生成环节，直接用MATLAB自带函数contour代替Runge-Kutta法自动生成等值线(即流线).

%主程序

load('PFD_My_Test_by_Stream_Function.mat');%载入矢量场数据

U(15)=U(15)-0.001;%对参数值1进行处理
V(15)=V(15)-0.001;
[x,y]=meshgrid(U,V);

%3*3次B样条基函数
k=3;l=3;
%控制点阵为7*7
m=6;n=6;
%KTP确定节点矢量
knotU=KTP(m,U);
knotV=KTP(n,V);

%构造线性方程组HX=B
H=[];B=[];
for i=1:length(U)
    for j=1:length(V);
    u0=U(i);v0=V(j);
    
    b=[0,1;-1,0]*PFD{i,j}';%注意修改
    A=coefficient_matrix(k,l,m,n,knotU,knotV,u0,v0);
    H=[H;A];
    B=[B;b];
    end
end

X=pinv(H)*B;%未知系数的列向量

%绘制B样条基函数表示的流函数
figure();
D=reshape(X,m+1,n+1);
z=zeros(length(U),length(V));
for i=1:length(U)
    for j=1:length(V)
        z(i,j)=streamfunction_matrix(U(i),V(j),k,l,m,n,knotU,knotV,D);
    end
end

surf(x,y,z)
title('流函数拟合')
figure();
du=zeros(length(U),length(V));dv=du;
for i=1:length(U)
    for j=1:length(V)
        du(i,j)=-B((i-1)*2*length(U)+2*j);
        dv(i,j)=B((i-1)*2*length(V)+2*j-1);
    end
end
du=du';dv=dv';


quiver(x,y,du,dv);%绘制矢量场图
hold on

title('矢量场图')
figure();
quiver(x,y,du,dv);%绘制矢量场图
hold on
[c,h]=contour(x,y,z,20);%绘制拟合函数流线图
clabel(c,h);
title('拟合函数流线图')

function U = KTP(n,T)

%KTP方法计算3次B样条曲线节点向量U
%n=控制顶点-1
%m=数据点个数-1
p=3;
m=length(T)-1;
U=zeros(1,n+p+2);
if m==n
    for j=1:n-p
        U(j+p+1)=1/p*sum(T(j+1:j+p-1+1));
    end
else  
    c=m/(n-p+1);
    for j=1:n-p
        i=fix(j*c);
        alpha=j*c-i;
        U(p+j+1)=(1-alpha)*T(i-1+1)+alpha*T(i+1);%参数标号从0开始，matlab从1开始记
    end
end
U(n+2:n+p+2)=1;

function A=coefficient_matrix(k,l,m,n,knotU,knotV,u0,v0)

%生成每个采样点处的系数矩阵A与右端项b

A1=zeros(m+1,n+1);A2=zeros(m+1,n+1);

for j=0:n
   
    A1(1,j+1)=-k*BaseFunction(j,l,v0,knotV)*BaseFunction(1,k-1,u0,knotU)/(knotU(k+2)-knotU(2));
    A1(m+1,j+1)=k*BaseFunction(j,l,v0,knotV)*BaseFunction(m,k-1,u0,knotU)/(knotU(m+k+1)-knotU(m+1));
    for i=1:m-1
        A1(i+1,j+1)=k*BaseFunction(j,l,v0,knotV)*(BaseFunction(i,k-1,u0,knotU)/(knotU(k+i+1)-knotU(i+1))-BaseFunction(i+1,k-1,u0,knotU)/(knotU(i+k+2)-knotU(i+2)));
    end
end
Au=[];
for i=1:m+1
    Au=[Au,A1(i,:)];
end

for i=0:m
    A2(i+1,1)=-l*BaseFunction(i,k,u0,knotU)*BaseFunction(1,l-1,v0,knotV)/(knotV(l+2)-knotV(2));
    A2(i+1,n+1)=l*BaseFunction(i,k,u0,knotU)*BaseFunction(n,l-1,v0,knotV)/(knotV(n+k+1)-knotV(n+1));
    for j=1:n-1
        A2(i+1,j+1)=l*BaseFunction(i,k,u0,knotU)*(BaseFunction(j,l-1,v0,knotV)/(knotV(j+l+1)-knotV(j+1))-BaseFunction(j+1,l-1,v0,knotV)/(knotV(j+l+2)-knotV(j+2)));
    end
end
Av=[];
for i=1:m+1
    Av=[Av,A2(i,:)];
end

A=[Au;Av];

end

function Psi=streamfunction_matrix(u,v,k,l,m,n,knotU,knotV,D)

%流函数表达式

M=zeros(1,m+1);N=zeros(1,n+1);
for i=0:m
    M(i+1)=BaseFunction(i,k,u,knotU);
end
for j=0:n
    N(j+1)=BaseFunction(j,l,v,knotV);
end
Psi=M*D*N';
end
    

function [Nip_u]=BaseFunction(i,p,u,NodeVector)
%利用de Boor-Cox 公式计算基函数Ni_p(u),i是节点序号,p是阶数,NodeVector为节点向量
%采用递归方式实现
if p == 0
    if (u >= NodeVector(i+1)) && (u < NodeVector(i+2)) %节点序号从0开始，但matlab从1开始，所以这里用i+1
        Nip_u = 1;
    else
        Nip_u = 0;
    end
else
    length1 = NodeVector(i+p+1) - NodeVector(i+1);
    length2 = NodeVector(i+p+2) - NodeVector(i+2); %支撑区间长度
    if length1 == 0  %规定0/0=0
        length1 = 1;
    end
    if length2 == 0
        length2 = 1;
    end
    Nip_u=(u-NodeVector(i+1))/length1*BaseFunction(i,p-1,u,NodeVector)+...
        +(NodeVector(i+p+2)-u)/length2*BaseFunction(i+1,p-1,u,NodeVector);
end    
end

运行结果

Reference

[1] 孙玉文，徐金亭，任斐，郭强著. 复杂曲面高性能加工技术与方法[M]. 北京：科学出版社, 2014:84-88.
[2] Sun, Yuwen,Sun, Shuoxue,Xu, Jinting,Guo, Dongming.A unified method of generating tool path based on multiple vector fields for CNC machining of compound NURBS surfaces[J],COMPUTER-AIDED DESIGN,2017,91:14-26.