张现东，卜 昆，刘连喜，窦杨柳

（1.西北工业大学现代设计与集成制造技术教育部重点实验室，西安 710072；2.北京星航机电装备有限公司，北京 100074）

张现东

西北工业大学航空宇航制造工程专业博士研究生，主要研究方向为CAD/CAM、精密测量、熔模精密铸造模具设计与优化。

空心涡轮叶片成形精度低，一直是困扰我国航空制造业的难题。与国外高达80%的合格率相比，国内的空心涡轮叶片精密铸造合格率不到40%[1]，差距很大。统计结果表明：在涡轮叶片无余量熔模精密铸造工艺过程中，同批次不合格空心叶片铸件中由于壁厚尺寸漂移引起的形状超差占50%左右[2]，是造成我国空心涡轮叶片合格率偏低的主要原因之一。涡轮叶片的壁厚是由外型与内腔共同保证的，内腔是精铸过程中由陶芯作为转接件形成的，因此陶芯的弯扭变形直接关系到空心涡轮叶片的壁厚分布情况。精确计算陶芯的弯扭变形，对于指导陶芯在蜡模型腔中的定位、控制蜡型壁厚、提高空心涡轮叶片成形精度都有重要的作用。

叶片类截面形状的弯扭变形一般认为是相对于理论截面的刚体位移，可以分解为相对于理论截面的平移与扭转。目前，国内外对叶片弯扭分析有大量的研究，王仲奇等[3]研究了弯扭叶片对降低能量损失的作用，并对弯扭叶片级和级组用于压气机做了大量理论研究；Houbolt等[4]将转子叶片的弯扭变形分为翼面向弯曲、翼弦向弯曲和扭曲变形；Mohaghegh等[5]通过定义叶片积叠轴，给出了基于积叠轴的叶片弯曲变形的定义；Omer[6]研究了叶片轮廓度误差、弯曲度误差和扭曲度误差定义及测量方法；陈福兴等[7]研究了基于三坐标测量数据的叶片弯扭变形误差计算方法；王红霞[8]研究了基于位移场的叶片弯扭变形分析；彭志光[9]研究了基于凸包算法的叶片弯扭变形计算方法。上述所有方法都是针对叶片类零件提出的，并不适用于结构更加复杂的陶芯。陶芯的制备及成形工艺在各国是保密的，对陶芯成形精度的研究更是很少见报道，陶芯弯扭变形的研究仅有黄胜利[10]、DOU等[11]提出的基于主动轮廓模型的陶芯弯扭变形计算方法，该方法通过Snake模型逼近陶芯外轮廓来获得陶芯截面外轮廓线，但Snake模型在陶芯纵向肋间易收敛成直线，导致提取的轮廓线光顺度不好，从而影响弯扭变形的计算精度。本文根据弯曲变形和扭曲变形的定义，研究和比较了几种弯扭变形算法，并通过几组截面轮廓线计算分析了陶芯的弯曲变形和扭曲变形。

弯曲变形和扭转变形的定义

二维平面内的弯曲变形和扭转变形分别表示为测量截面相对于理论截面的平移位移和扭转角位移。如图1所示，实线表示理论CAD截面轮廓线，虚线表示测量截面线。图1（a）中点Mc为测量截面的形心，点Dc为CAD截面的形心即为截面弯曲变形矢量，用Tc来表示。Tc可以分解为沿X轴向的位移Tx和沿Y轴向的位移Ty。图1（b）中，点A和点B分别为CAD截面的后缘点和测量截面后缘点（中弧线延长线与后缘轮廓线的交点），直线McB与DcA的夹角θ表示测量截面相对于CAD截面的扭转变形，θ逆时针为正，顺时针为负。

图1 弯扭变形定义Fig.1 Definition of bending and twist deformation

弯扭变形的计算方法

空心涡轮叶片陶芯包含众多细小结构：数个纵向肋槽，数十个横向肋槽，近百个扰流柱孔、排气边、定位窗等，截面轮廓复杂，增加了弯扭变形的计算难度。本文根据陶芯的结构特点，将陶芯弯扭变形的计算方法分为两类：基于轮廓线的算法和基于数据点的算法。基于轮廓线的算法是通过构造实测点云截面轮廓线，并与CAD截面轮廓线对比，根据弯扭变形定义计算弯扭变形，主要有基于主动轮廓模型的算法和基于样条拟合的算法等。基于数据点的算法直接通过实测点云与CAD截面轮廓线计算弯扭变形，主要有凸包算法、二维配准算法、距离权值法等。

1 基于轮廓线的算法

1.1 主动轮廓模型提取轮廓线

主动轮廓模型也称为Snake模型,是由Kass等[12]在1987年提出的,广泛应用于数字图像分析和计算机视觉领域。点云数据中的Snake可以用离散点pi和连接邻点的直线段li表示：

其中，li= {pi，pi+1}（i=1，2，…，n-1）。

基于Snake的二维点云数据轮廓线提取主要包括3部分：初始轮廓点的选取、能量方程的表示、轮廓点运动方向的确定。Snake提取陶芯截面轮廓线的步骤为：

（1）提取STL格式涡轮叶片外型截面轮廓点，并采用Kochanek曲线插值加密，以加密后的点作为Snake初始轮廓点。

（2）建立带权值的Snake能量方程，

其中，Etension和Ebend分别表示由内力（包括拉伸力和弯曲力）作用引起的拉伸能和弯曲能，Edist表示点云自身的特征能。

其中，θi是li-1与li的夹角，其中pi是初始轮廓上的第i点，mj是点云数据中的第j点，N1、N2分别代表初始轮廓点个数和点云数据点个数，|pimj|代表点pi到mj的距离。

（3）以涡轮叶片截面线在初始轮廓点处内法矢方向作为该轮廓点的运动方向。

（4）如果所有轮廓点的弯曲能都小于给定阈值，即max{Edist(i)}＜ε，则输出Snake；否则，增加初始轮廓点数，重复步骤（1） ～（4），直至收敛。

1.2 样条拟合轮廓线

样条拟合陶芯轮廓线包含4步：截面数据的获取；截面轮廓数据的提取；型面轮廓线拟合；后缘拟合。

（1） 截面数据提取：通过ATOS或ICT获得陶芯表面三维点云数据，将三维点云数据与陶芯CAD模型配准，采用切片法获得陶芯截面数据（图 2（a））；

（2） 截面轮廓数据提取：计算截面数据点到CAD截面轮廓线段（图2（b））的距离，给定距离阈值，提取截面轮廓数据（图2（c））；

（3）截面轮廓线拟合：对截面数据排序，然后以NUBRS拟合陶芯型面轮廓线；

（4）后缘拟合：由于排气边的影响，陶芯后缘测量数据较少（图3（a）），直接采用样条拟合偏差较大（图3（b）），理论上将陶芯后缘近似成圆弧，可采用最小二乘圆的拟合陶芯后缘轮廓线（图3（c））。

如图4所示，AO'和BO'分别为陶芯叶背延拓线和叶盆延拓线，内切圆与AO'和BO'分别切于点C和点D，O'C的长度为a，e1和e2分别为O'A和O'B的单位方向法矢。则内

图2 截面外轮廓数据点提取Fig.2 Abstraction of measured contour points for cross-section

图3 最小二乘法拟合陶芯后缘轮廓线Fig.3 Trailing edge fitting by least-squares method

图4 叶盆叶背延拓线内切圆Fig.4 Inscribed circle of blade basin and back extension cord

切圆方程可表示为：

将后缘点坐标代入内切圆方程，就可得到后缘点到内切圆的代数距离pi（a）：

后缘数据点的坐标变换到坐标系X'O'Y'下，带入式（7），记方程的>值为p1（a），p2（a），p3（a），…，以min{p1（a）+p2（a）+p3（a）+…}为目标函数，求取满足最优解的a，即得到距所有后缘点距离和最小的圆。最后，将所得圆心坐标变换到陶芯模型坐标系中，得到陶芯后缘的拟合轮廓线。

之后，采用等半径法计算中弧线及后缘点，采用等弧长或等参数的方式采样轮廓线，采样点坐标均值作为形心。采用同样的方法得到CAD截面线的后缘点及形心，按照弯曲和扭曲变形的定义计算陶芯的弯扭变形。

2 基于数据点的计算方法

2.1 凸包算法

文献[9]采用凸包算法求解叶片测量数据形心，首先通过凸包计算将叶片测量数据点分为凸包点和凹弧点，然后采用分割三角形法分别计算凸包点和凹弧点的形心，最后通过面积加权计算出叶片截面测量数据点的形心。由于陶芯纵向肋结构会产生肋面数据点，可采用1.2节中提取截面轮廓点的方法去除内部数据点，然后计算凸包。凸包算法计算陶芯测量数据形心的步骤为：（1）提取截面轮廓数据（图 2（c））；（2）计算凸包（图 5（a））；（3）采用三角形分割分别计算凸包点（图5（b））和凹弧点形心（图 5（c））；（4）计算陶芯截面形心。

设（xi，yi），（i=1，2，…，n）为凸包点，凸包第i个三角形顶点分别为（x1，y1）、（xi+1，yi+1）和（xi+2，yi+2），面积为Ai=0.5|x1（yi+1-yi+2）+xi+1（yi+2-y1）+xi+2（y1-yi+1）|；（x'j，y'j），（j=1，2，…，m）为凹弧点，凹弧第j个三角形顶点为（x'1，y'1）、（x'j+1，y'j+1）和（x'j+2，y'j+2），面 积为Bj=0.5|x'1（y'j+1-y'j+2）+x'j+1（y'j+2-y'1）+x'j+2（y1-y'j+1）|。其中n、m分别为凸包点和凹弧点个数，点（x1，y1）和点（x'1，y'1）相同，为凸包（凹弧）起点；点（xn，yn）和点（x'm，y'm）相同，为凸包（凹弧）终点。设陶芯截面形心坐标（xc，yc）为：

图5 陶芯测量数据凸包计算Fig.5 Convex hull of measured contour points

凸包算法的扭转变形，以CAD截面弦线与测量截面弦线（凸包多边形的最大边长）夹角来表示。

2.2 二维配准算法

二维配准法直接采用SVD-ICP算法对实测截面点云与CAD截面配准，在配准的过程中记录平移及旋转量，累加平移量为弯曲变形量，累加旋转量为扭转变形，其流程如图6所示。

图6 二维配准算法流程Fig.6 Flowchart of 2-D registration method

2.3 距离权值算法

距离权值法与凸包算法类似，不通过拟合截面线直接计算陶芯截面测量数据形心；不同点是凸包算法将陶芯截面看作片体，用三角分割法计算形心，距离权值法将陶芯截面轮廓看作线体，以相邻点间的距离作为权值计算形心。算法步骤如下：

（1）提取陶芯截面轮廓数据（图2（c））；

（2）对截面轮廓数据排序；

（3）计算距离加权形心。

设排序后截面轮廓数据点为p1，p2，p3，…pn，令pn+1=p1，对点序列p1，p2，p3，…，pn，pn+1计算距离加权形心pc：

式中，di为点pi和点pi+1的距离。

以测量数据点中距离形心最远的点近似后缘点计算扭转变形。

仿真计算

为验证上述算法的准确性，截取某型号陶芯CAD模型Z=50截面，均匀离散数据点，将点集沿X轴和Y轴各偏移0.5mm，绕Z轴旋转2°，如图7，红色实线为CAD模型，黑色点为变换后的数据点。

变换后的数据点视为仿真测量点，采用Snake模型算法、样条拟合轮廓线法、凸包算法、二维配准算法、距离权值算法，进行弯扭变形计算分析。

算例1：对仿真数据不添加随机误差，计算结果如表1所示。

算例2：对数据点集添加服从正态分布N（0，0.01）的随机法向误差，计算结果如表2所示。

算例3：对数据点集添加服从正态分布N（0，0.1）的随机法向误差，计算结果如表3所示。

图7 弯扭变形仿真数据Fig.7 Simulated points with known bending and twist

通过表1、2、3看出，当测量数据不存在误差时，基于数据点的弯扭变形算法的计算精度明显高于基于轮廓线的算法精度；当测量数据误差较小时，5种算法的误差均不大，样条拟合算法和二维配准算法计算精度略高于其他3种算法；当测量数据误差较大时，Snake模型、二维配准、凸包算法、距离权值算法的最大相对误差均超过或接近10%，而样条拟合算法的相对误差低于这4种算法，对测量数据点误差具有更好的计算稳定性。

比较Snake模型算法与样条拟合算法，样条拟合算法对弯曲和扭曲变形的计算精度均高于Snake模型算法，因为Snake模型易在陶芯纵向肋间收敛成直线，导致提取轮廓线的光顺度不是很好；凸包算法和距离权值算法的弯曲变形计算精度要高于扭曲变形计算精度，由于这两种算法是通过所有（或大部分）点计算形心的，而扭曲变形是通过前后缘与弦线切点或后缘点计算的，所以计算扭曲变形对点的误差更敏感。

表1 算例1计算结果

表2 算例2计算结果

表3 算例3计算结果

复杂陶蕊实测弯扭分析

对某型号涡轮叶片进行CT扫描点云提取7组水平截面，采用样条拟合轮廓线的方法计算弯扭变形。水平截面位置如图8（a）所示，截面线如图8（b）所示，图8（b）中实线为CAD轮廓线，虚线为拟合的测量截面线。弯扭变形计算结果见表4。

将表4中的弯扭变形计算结果绘制成弯曲变形、扭转变形曲线，如图9所示。

图8 截面位置及拟合轮廓线Fig.8 Location and contour lines of cross-sections for ceramic core

表4 截面平移与扭转结果

综合表4及图9可以得到陶芯的弯曲和扭曲变形的整体趋势，发现存在弯扭不连续截面，说明陶芯在制坯、烧结和强化处理工艺过程中有局部变形，陶芯的模具设计和生产工艺还有待改进。下一步工作是研究陶芯的变形机理和基于变形调控的陶芯模具优化设计及工艺改进技术。

结论

本文基于刚体位移的假设，给出了陶芯弯扭变形的定义；根据陶芯的结构特点，研究了基于轮廓线和数据点的弯扭变形计算方法，通过仿真分析和试例验证可得出：（1）当测量数据点与CAD模型完全对应不包含误差时，基于数据点的凸包算法、二维配准与距离权值算法对弯扭变形的计算精度要高于基于轮廓线的算法；（2）样条拟合轮廓线法对数据点误差的计算稳定性要高于Snake模型、凸包算法、二维配准和距离权值算法；（3）样条拟合轮廓线法对弯扭变形的计算精度要高于Snake模型算法，这是由于样条拟合轮廓线的光顺度要高于Snake模型；（4）凸包算法和距离权值算法对扭曲变形的计算精度要低于弯曲变形，这是由于两种算法计算扭曲变形对测量数据点的误差更敏感。

图9 实测陶芯弯曲变形曲线Fig.9 Bending and torsion curves of measured ceramic core

样条拟合轮廓线法能准确地计算陶芯的弯曲、扭曲变形，而且对后续陶芯变形位移场建模、轮廓度误差分析及烧成收缩率计算都很方便，具有较大的工程实用价值。

[1]李世峰,卜昆,董一巍,等.空心涡轮叶片近净形过程变形分析与控制[J].中国电机工程学报,2011,31(14):109-112.

LI Shifeng,BU Kun,DONG Yiwei,et al.Analyzing and controlling deformation of hollow turbine blade in near net-shape process[J].Proceedings of the CSEE,2011,31(14):109-112.

[2]WINENS K.Airfoil thickness as a life limiting factor of gas turbine blades[C]//20th Symposium of the Industrial Application of Gas Turbines Committee,Banff,Alberta,Canada,2013: 1-14.

[3]王仲奇,苏杰先,钟兢军.弯扭叶片栅内减少能量损失机理研究的新进展[J].工程热物理学报,1994,15(2):147-152.

WANG Zhongqi,SU Jiexian,ZHONG Jinjun.New process of investigation into mechanism of reducing energy loss in cascades with curved and twisted blades[J].Journal of Engineering Thermophysics,1994,15(2):147-152.

[4]HOUBOLT J C,BROOKS G W.Differential equations of motion for combined flapwise bending,chordwise bending,and torsion of twisted nonuniform rotor blades[M].Washington: NACA ,1957.

[5]MOHAGHEGH K,SADEGH M H,ABDULLAH A,et al.Improvement of reverseengineered turbine blades using construction geometry[J].The International Journal of Advanced Manufacturing Technology,2010,49:675–687.

[6]OMER L H.Recent advances in laser triangulation-based measurement of airfoil surfaces[R].Industrial Optical Sensors for Metrology and Inspection.Proceedings of SPIE,1995.

[7]陈福兴,张秋菊.基于UG的叶片型面质量评估[J].工具技术,2006,40(8):66-69.

CHEN Fuxing,ZHANG Qiuju.Evaluation of blade surface quality based on UG[J].Tool Engineering,2006,40(8):66-69.

[8]王红霞.基于位移场的涡轮叶片变形特征分离技术研究[D].西安: 西北工业大学,2009:31-42.

WANG Hongxia.Separation of the deformation characteristics based on the displacement field of the turbine blade[D].Xi’an: Northwestern Polytechnical University,2009:31-42.

[9]彭志光.基于改进凸包算法的叶片型面检测[D].武汉: 华中科技大学,2012:51-54.

PENG Zhiguang.The Blade surface inspection based on the improved convex hull algorithm[D].Wuhan: Huazhong University of Science and Technology,2012:51-54.

[10]黄胜利,卜昆,徐岳林,等.基于主动轮廓模型的复杂陶芯弯扭变形分析[J].计算机集成制造系统,2010,16(8):1643-1648.

HUANG Shengli,BU Kun,XU Yuelin.Analysis of bending and torsion deformation for complex ceramic core based on active contour model[J].Computer Integrated Manufacturing Systems,2010,16(8):1643-1648.

[11]DOU Y L,XU Q Y,BU K.Reversing design methodology of ceramic core for hollow turbine blade based on measured data[C]//ASME.Proceedings of ASME 2014 International Mechanical Engineering Congress and Exposition,Quebec,Canada,2014.

[12]KASS M,WITKIN A,TERZOPOLOS D.Snake: active contour model[J].International Journal of Computer Vision,1988,1(4):321-331.