翼型叶片类曲面B样条曲线拟合算法研究

李传军,王立萍

(1.天津中德应用技术大学 机械工程学院,天津 300350;2.天津中德应用技术大学 基础实验实训中心,天津 300350)

0 引言

进行高性能的机翼和翼型设计是现代飞行器面临的两个重要任务,如何在满足机翼几何参数要求条件下,对翼型进行高性能设计是重要的技术指标。飞行器机翼沿着来流的方向切下来的剖面称为翼型,翼型是飞行器产生升力的基本单元,为了取得较大升力,其流线形外形一般制作成头部圆滑,尾巴尖瘦。翼型轮廓线的内切圆的圆心连线称为中弧线,不同的翼型沿翼弦的厚度分布规律不同。翼型叶片类曲面零件其表达形式不能用有理多项式曲线段全部表示,当约束条件较多时,需要利用高次曲线,高次曲线在拟合局部复杂形状时有较高精度,但是高次曲线计算机处理效率低。翼型叶片类曲面零件在基于B样条曲线的轨迹生成过程中通过修改控制顶点具有全局性,为了实现高效率、高稳定和局部控制,采用基于分段多项式的B样条基函数实现。B样条曲线是控制点P在B样条基函数约束下的控制多边形。按照变量[0,1]范围内计算B样条基函数并将非零的基函数值与对应的控制点相乘求和得到B样条曲线[1]。

叶片类曲面的数控加工通常用普通、传统CAM刀具路径规划,用五轴机床进行加工。样条曲线[2]本身具有G1或更高阶的几何连续,消除了微小直线段之间的段间转接,因此加工表面会更加工光滑,加工质量会得到提高。为了使加工轨迹光顺,现在主要的光顺方法包括最小二乘法、能量法和小波方法等。能量法[3]是应用最为广泛的一种曲面光顺处理方法,其基本思想是使曲面的整体能量在一定约束条件下达到最小,由于复杂曲面加工时计算速度慢,效率不高,并且只从型面上进行要求,没有考虑受力合理性。

本文从流体力学性能角度分析叶片类曲线曲面,研究其外形参数与流体力学性能的关系;提出叶片类曲面以流体力学为约束,去光顺这些不光滑的轨迹;将力学约束应用到加工轨迹上,建立一种B样条曲线等效升力模型,参考流体力学中的升力原理,通过翼型叶片类曲面,既满足光顺的特性,又符合力学特性。

1 流线与迹线

流线是流场中某一瞬时的空间曲线,按照跟踪质点的描述法,记录质点的位移、速度和加速度和其他物理参数的变化,即流线的切线方向和速度流团的变化情况,分为定常流动和非定常流动两种[4]。迹线是以流场中的流体微团为研究对象,用拉格朗日描述法确定其轨迹的运动方向,如图1所示。

图1 流线与速度矢量

设ds表示流线上的微段,按照流线的研究对象,其速度v与轨迹ds平行,

ds×v=0。

(1)

按照右手坐标系在(i,j,k)三个方向上展开得:

ds=dxi+dyj+dzk;

(2)

v=ui+vj+wk;

(3)

i(wdy-vdz)+j(udz-wdx)+k(vdx-udy)=0。

(4)

得出:

(5)

即为流线微分方程,对于任意二维情况,流线微分方程为:

vdx-udy=0;

(6)

(7)

式(7)表示其速度方向与流线方向相切。

对于二维定常不可压缩的无旋流动,连续性方程可写为:

(8)

根据偏微分方程可知,式(8)表明udy-vdx是任意一个函数ψ的全微分,即:

(9)

函数ψ(x,y)即为流函数。流函数ψ(x,y)为常数的曲线求导即得到式(9),因此曲线ψ(x,y)=C即为一条流线,ψ(x,y)=C即为流线方程,流函数满足拉普拉斯方程:

(10)

根据推导过程,确定在流线位置与质点运动的位置变量是一一对应的连续函数,便于离散刀具轨迹与流体力学特性的升力约束建立联系。

2 B样条曲线的等效升力模型

2.1 B样条曲线

B样条曲线是为了解决描述复杂形状的连接问题而提出,其整体表示时具有局部性质和Bezier曲线的优点,使其有统一、通用和有效的标准算法和解决方案,是目前最广泛和流行的形状数学方法,并且B样条方法也成为工业产品几何定义的国际标准。用最常用的递推公式方法对B样条基函数Ni,p(u)进行定义。设U为有m个实数的节点矢量,其中ui为矢量节点,Ni,p(u)根据式(11)递推,其表达形式为[5]:

(11)

利用式(11)的B样条基函数,根据参考文献[6-10]提出的B样条曲线说明,定义p次B样条曲线C(u),其表达形式为:

(12)

其中:U表示节点矢量;Pi为B样条曲线的控制点;p为曲线次数;(n+1)为控制顶点数;(m+1)为节点数。

根据B样条基函数的求导公式,得到其B样条曲线的k阶导数为:

(13)

通过以上分析可得,一个p次与p-1次B样条曲线的导数在同一个节点向量上,通过递归方法可以用计算机得到更高阶导数。

2.2 二维翼型的升力模型

叶片类曲面功能曲线的型线形状,决定其截面的流体力学性能要求,型值点组成的形状是相对固定的,受到流体力学约束。以发动机动子叶片的受力情况为例,其截面近似于一条表示翼型等叶片类曲面的B样条曲线,其气动升力是由曲线几何、流体条件等因素决定的,本文参考流体力学升力原理[11-16],建立一种B样条曲线等效升力模型。

二维机翼工作原理如图2所示,当气流以均匀速度吹向翼型曲面时,在各个截面处会对机翼产生向上的升力Fa和向后的阻力Fu,由于机翼升力的产生是一个非常复杂的纳维—斯托克斯方程,其力学建模和求解非常复杂,本文仅考虑二维翼型的情况并对其进行线性化。

图2 二维翼型工作原理

根据库塔—茹科夫斯基定理,对于低速无黏均匀来流中的二维机翼,其单位展长上的升力为:

Fa=ρv0Γ。

(13)

式中:ρ为流体密度,v0为相对流体速度,Γ为绕机翼的环量。

假设翼型为薄翼型,在低速定常的不可压缩流体中相对运动,迎角为0,翼型截面可以看为是由对称翼型和中弧线叠加而成,如图3所示。

图3 对称翼型上半部分

根据薄翼型理论,其向上的升力可近似为有迎角中弧线柱面升力和无迎角对称翼型曲面升力的线性叠加,而无迎角时对称翼型曲面的升力为0,则翼型的总升力可看作是由带迎角的中弧线柱面产生的,如图4所示为绕机翼的环量表示,根据中弧线定义和表示方法,等效为如图5所示。

图4 绕机翼环量表示

图5 中弧线柱面线性化表示

中弧线的速度环量为

(14)

其中γ(s)为涡流强度分布,s为弧长参数。利用三角级数的形式进行展开,涡流强度分布为:

(15)

(16)

其中b为弦长长度。

升力为:

(17)

2.3 基于B样条曲线的等效升力

利用B样条曲线来表示中弧线:

(18)

并对节点矢量U和控制顶点Pi(Xi,Yi)进行约束,使其满足x关于参数u为线性关系,即

(19)

则有

(20)

可得:

(21)

(22)

故翼型的升力用B样条曲线可表示为

(23)

可以得出利用B样条曲线描述一条薄翼型的中弧线,当B样条的节点矢量不变时,其二维翼型单位展长的升力Fa和B样条的控制顶点{Yi}是线性相关的。定义B样条等效升力为:

(24)

将其作为优化目标,可使曲线趋向于流线型,即可获得更加光顺的B样条曲线。

3 B样条刀轨离散刀位点识别和重构

3.1 离散刀轨B样条曲线重构流程

传统的样条曲线重构仅考虑了几何约束,利用最小二乘法使型值点和曲线的误差平方和最小,存在局部曲率波动、刀轨曲线光顺性差的情况,插补时会出现局部的加减速,从而降低了加工效率和加工精度。本文考虑B样条曲线的流体力学特性,对叶片类曲面离散刀轨的B样条曲线重构引入升力最优优化目标,以获得具有更加光顺的B样条刀轨曲线。曲线重构的流程图如图6所示。首先对离散刀位点进行分段识别,获取无拐点的曲线段;其次对每一段的离散点进行坐标变换,将起末点变换至翼型体轴坐标系的X轴上,并对其进行参数化,确定B样条曲线的节点矢量;最后建立B样条曲线的约束模型和优化模型,在体轴坐标系中对B样条曲线进行优化重构,然后通过坐标变换变换至工件坐标系,即完成离散刀位点的B样条重构。

图6 离散刀轨B样条曲线重构流程图

3.2 离散刀位点识别

根据离散刀位点间的几何特征,可将离散刀位点集分为长直线段和无拐点微小直线段集,其中无拐点微小直线段集可以利用B样条曲线进行重构。长直线段和无拐点微小直线段集之间的转接关系可分为长直线—长直线、长直线—微段集和微段集—微段集3种形式,如图7所示。

图7 叶片类曲面离散刀位点分段转接形式

利用最小转角限制和最大弓高误差限制对离散刀位点进行分段。如图8所示为曲线分段点示意图,实线为指令轨迹,虚线为实际轨迹,则最小转角限制为:

图8 长直线段轨迹表示

θcr=2×arccos(Kpε/F/(1-Kf))。

(25)

式中:F为指令速度,ε为预设最大转角误差,Kf和Kp分别为伺服前馈增益和位置增益。

然后计算两个最大弓高误差限制为:

δ1=R(1-cosα1),

(26)

δ2=R(1-cos(π-θ-α1))。

(27)

其中:R=2×l1/sinα1,α1=arctan(l1sin(π-θ)/(l2+l1cos(π-θ))),l1和l2分别为前后曲线段的长度,θ为前后两曲线段的切矢夹角。

再根据曲率的极小值为分段点,作为无拐点微小直线段集的充分条件。最终获得了一组由长直线段点集和无拐点微小直线段点集组成的轨迹段。

3.3 型值点的坐标系变换

翼型叶片类曲面,其上的型值点在工件坐标系下的坐标已知,只要获得体轴坐标系到工件坐标系的转换矩阵,就可以得到型值点在体轴坐标系下的坐标值。其转换矩阵与叶片的前后缘点及弦线b长度有关。

根据合同变换和活动标架[16]的相关原理,给出空间坐标系R3中一个直角坐标系{O;e1,e2,e3},坐标系中任一点Q对应的坐标为{x1,x2,x3},矢量表示为:

(28)

(29)

根据合同变换矩阵定义得出:

(30)

其中aij是合同变换的正交矩阵。

其中正交矩阵表示为:

(31)

3.4 无拐点微小直线段集的B样条曲线重构

(32)

变换矩阵

(33)

根据中弧线X坐标关于参数u的线性化假设,根据体轴坐标系点集的X坐标进行参数化和归一化。

(34)

在体轴坐标系中,已知微小刀位点集{Qk}={Qk|k=0,1,…,l},采用X坐标参数化。用p次B样条曲线进行重构

(35)

其中Pi=(Xi,Yi),节点矢量U采用准均匀形式,Pi·X可根据参数化的Qk·X选取,以保证重构曲线的X坐标关于参数u的线性关系。

为保证重构的各段曲线之间的连续性,除了端点插值外,还应在首末端点处保证切矢连续,即:

(36)

曲线C(u)在型值点处的误差可由Y轴坐标表示:

(37)

将B样条曲线的等效升力作为优化目标,可以得到体轴坐标系中微小刀位点集{Qk}的B样条曲线重构模型

s.t.

C(0)=Q0;

C(1)=Ql;

k=0,…,l。

(38)

式(38)以求等效升力最大值作为优化目标,求解的最大值含有u值,要确定最大值,需根据u值范围进行分别代入求解,u∈(0,1),任意间隔代入求解判断。最小二乘分段曲线升力:{1096.24,637.073,755.709,172.428,12.5507,-263.627,-799.027,-432.151,-1052.46};本文算法分段曲线升力:{1096.56,640.764,755.835,238.701,13.0012,-257.671,-798.802,-590.045,-1052.13}。通过两者对比,本文算法计算的分段曲线等效升力优于最小二乘分段曲线的等效升力。

(39)

4 海星型轨迹的B样条等效升力拟合

节点矢量采用准均匀节点矢量方法。权因子全为1,即等权因子,构造的B样条曲线为均匀B样条曲线。曲线次数一般取固定值,通常3次就满足了数控加工中的速度、加速度和加加速度连续。取高次亦可,在仿真过程中采用3次进行仿真。

控制顶点的求解方法。采用顶点个数递增的方案,依次计算由给定初始个数条件下的控制顶点,计算完成后判断曲线是否满足精度要求,满足精度要求输出曲线,不满足精度要求则增加控制顶点,进行循环迭代,直至满足要求为止。

根据控制顶点,在已知次数、控制顶点个数,给定的节点矢量采用均匀节点矢量方法。总的节点矢量为p+n+1。p为曲线次数,n为控制顶点个数;则具体计算方法为:前p+1个节点矢量为0,最后的p+1节点矢量为1,中间节点均匀分布。

海星型轨迹是典型轨迹拟合参考曲线[16]。通过模拟基于流体力学特性的海星型离散刀轨的拟合过程,以离散海星型B样条曲线重构刀具轨迹为验证对象,根据B样条等效升力模型,通过控制弓形高度间的最大误差范围为0.001,获得200个基于B样条等效升力的离散点,如图9所示,通过离散点的B样条光顺拟合作为约束验证条件,对本文算法进行验证。

图9 基于流体力学特性的海星型B样条离散点

4.1 拟合效率对比

基于流体力学特性的B样条曲线采用等效升力进行光顺平滑拟合。为保证重构的各段曲线之间的连续性,给出了端点插值约束,采用准连续点和准均匀节点矢量,只要保证首尾控制顶点重合即可,即将分段拐点作为B样条曲线的控制顶点的始末点,保证整段曲线的G0连续;同时,分段函数的切矢连续,保证分段B样条曲线的G1连续;通过以上约束条件,保证海星分段曲线各段曲线首尾相连,且连接点处一阶导矢相同。如图10所示,本文根据离散点的曲率估计拐点个数为9个,两端点作为B样条曲线的拐点,曲线上其他点的切削矢量计算结果及拟合轨迹方向如图10a所示,根据拐点及切矢计算结果,本仿真算法是以其他9个离散点为转折进行光顺拟合,其拟合结果如图10b所示。

图10 海星型轨迹离散点拐点、切矢计算和拟合结果

最小二乘B样条拟合法和基于流体力学特性的等效升力最大B样条拟合法,对比两种方法的拟合计算效率。以Windows 10家庭中文版的测试操作系统,其配置为处理器为Intel(R) Core(TM)I7-8550U CPU@1.99 GHz E75,对两种方法的拟合方式分别进行500次拟合测试,用MATLAB 2017作为测试环境,用曲线轨迹数据压缩率[15]表示效率。曲线轨迹数据压缩率=(初始离散点数-重构的控制顶点数)/初始离散点数。曲线控制顶点为201个,分为9段进行拟合。最小二乘法与等效升力最大拟合法控制顶点如表1所示。

表1 最小二乘法与等效升力最大拟合法控制顶点

最小二乘法拟合控制顶点数133个,曲线轨迹数据压缩率为33.83%;等效升力最大拟合法拟合控制顶点数128个,曲线轨迹数据压缩率为36.32%,测试结果如表2所示。可见,本文的基于流体力学特性的升力最大拟合法在控制顶点个数以及曲线轨迹数据压缩率方面均优于最小二乘拟合法。

表2 拟合效率比较

4.2 加工效率对比

以离散微直线段刀具轨迹和基于流体力学特性的B样条曲线轨迹,对比加工效率,说明基于流体力学特性的B样条曲线算法在数控加工中的高效性。按照文献[19-25]中采用分段S型加减速,小直线段采用S型加减速,B样条采用直接插补方法进行速度规划。速度规划是以最小时间为目标函数,虽然速度前瞻可以减少进给速度波动,但不能避免微小线段段内加减速引起的微小波动。微小波动是由于段内引起的,但由于微小直线段间会有速度方向的突变,段内会有加减速。结合常规数控系统采用的进给速度规划,规划约束项如表3所示。

表3 基于流体力学特性的B样条曲线进给速度规划仿真参数

基于流体力学特性的B样条曲线和脉冲增量法离散微直线段的进给速度曲线仿真如图11所示。图11a采用脉冲增量法离散微直线段刀具轨迹加工时间为7.911 s,由于加速度和加加速度变化明显,速度不断变化,加工效率不高。而B样条光顺压缩刀具轨迹仅需4.894 s,基于流体力学特性的B样条曲线在样条曲线拟合的段间和段内速度切向矢量连续,进给速度曲线更光顺,如图11b所示,因此加工效率高。

图11 基于流体力学特性的B样条曲线和脉冲增量法离散微直线段的进给速度曲线仿真

4.3 拟合精度对比

当已知控制顶点时,可以计算出每个型值点处的误差值,即给定的201个离散点的误差值(首点、尾点重复,保证连续),则平均误差和最大误差均可以计算,并通过图形显示出来,如图12所示。本文的优化目标是升力最大,一般曲线上的最大误差值为整条曲线的精度。基于等效升力最大算法拟合误差其中最大误差为0.001,平均误差为8.09×10-4。最小二乘法的最大误差0.001,最小二乘法的平均误差是2.12×10-4。

图12 两种拟合算法轨迹误差对比

通过海星型轨迹拟合和计算效率分析,本文提出基于流体力学特性的B样条曲线光顺拟合算法在离散刀具轨迹拟合重构和轨迹的光顺应用是可行有效的。在满足离散轨迹最小误差要求的约束条件下,重构的刀具轨迹根据仿真实验验证可以满足误差要求。通过仿真光顺拟合后的刀具轨迹具有C2连续,速度曲线更加光顺平滑,提高了数控机床的加工效率。

5 平面叶栅叶片曲面加工试验对比

5.1 试验加工

以亚音速压气机叶片翼型的平面叶栅叶片为加工对象,分别采用直线插补和B样条曲线直接插补两种不同插补方式进行加工,以验证三次B样条曲线直接插补加工的加工效果。

(1)选用最大厚度相对位置0.4的NACA67-010叶型如图13和图14所示,叶型选用零攻角设计,中弧线选用双圆弧,加密前后缘数据点。

图13 叶片型值点分布图

图14 平面叶栅叶片模型

(2)叶栅的精加工分为三轴叶背加工与三轴叶盆加工,采用单向走刀顺铣加工,以获得较好的加工表面,如图15所示为叶栅叶片叶背精加工刀轨图,如图16所示为叶栅叶片叶盆精加工刀具轨迹图。

图15 叶栅叶片叶背精加工刀轨

图16 叶栅叶片叶盆精加工刀轨

(3)加工切削参数如表4所示。由于切深和切宽较小,对刀具每齿进给量进行切削半径补偿,在图17所示自主开发的数控机床中设置G54刀具半径偏置数值,如图18所示。

表4 叶片加工切削参数

图17 工作台平动的AC双转台五轴联动机床

图18 自主开发的数控机床刀具偏置界面

(4)将刀具轨迹导出生成刀位文件,提取刀具轨迹离散点,利用本文提出的等效升力拟合方法进行离线拟合,生成样条曲线格式的NC程序,再以鲁南三轴数控加工中心XK7132A为基础,再工作台上加装AC旋转台实现五轴联动结构,课题组系统自主开发的数控机床整体如图17所示实现。叶片加工的结果如图19和图20所示。

图19 叶栅叶片三轴精加工

图20 平面叶栅叶片加工结果

5.2 试验加工结果分析

通过三坐标测量仪对两种方法加工的叶栅叶片的加工截面轮廓度进行检测,并通过光学显微镜检测叶栅叶片的表面三维形貌。利用海克斯康三坐标测量仪分别对两种叶栅进行轮廓度检测,获得每组叶栅叶片的检测结果,如图21所示。叶栅的检测结果如表5所示。

表5 三坐标检测结果 mm

图21 第1号叶栅叶片三坐标轮廓检测结果图

通过奥林巴斯DSX500光学数码3D显微镜检测两种加工方法加工的平面叶栅叶片的加工表面三维形貌,检测结果统计信息如表6所示。3组叶栅叶片均沿流线方向走刀,设计行距0.523 mm,理论残留高度为0.023 mm,实际加工中,直线插补的残留高度大于理论残留高度,B样条曲线直接插补的残留高度小于理论残留高度。另外,B样条曲线直接插补方法加工的表面波峰的平均高度差1 μm,而直线插补方式加工的表面波峰平均高度差为2.8 μm。

表6 第2、6、7号平面叶栅叶片加工表面三维形貌检测结果统计

通过三维表面形貌对3种方式加工的叶片进行检测,结果如图22所示。从图中可以看出,直线插补的加工表面有明显分段的凹痕,不同分段表面内部也不平滑,这是因为微小直线段段间转接只能达到G0连续,且在转接点出会出现加减速,造成切削过程的不平稳;而B 样条曲线直接插补没有该现象,在进给方向上的刀痕均匀,曲线直接插补方式的加工表面质量比直线插补有一定的提高。

图22 叶片加工表面三维形貌检测立体图

由上述检测结果对比可知,本文提出的B 样条曲线直接插补加工的叶片表面三维形貌比直线插补方式加工的更光滑、均匀,能够获得更低的残留高度和更高的表面质量。

6 结束语

本文首先根据B样条曲线性质,将翼形等叶片类曲面的流体力学特性结合等效升力约束条件,描述一条薄翼形的中弧线。利用B样条曲线的等效升力,对基于流体力学特性的离散刀轨进行B样条曲线重构,利用曲率的极小值分段,建立满足端点插值、端点切矢连续和拟合最大误差下的等效升力刀轨模型。通过海星型轨迹的B样条等效升力拟合验证和平面叶栅叶片曲面加工试验对比,其计算效率和加工效率均有明显提高。基于B样条曲线对基于流体力学特性的叶片类曲面进行刀具轨迹拟合,由于其数学模型存在权因子、迭代计算和误差判断等内容,对实时性要求高的数控加工还要进一步研究。