空心渦輪葉片復雜陶芯彎扭變形分析方法比較 *

張現(xiàn)東，卜 昆，劉連喜，竇楊柳

（1.西北工業(yè)大學現(xiàn)代設(shè)計與集成制造技術(shù)教育部重點實驗室，西安 710072；2.北京星航機電裝備有限公司，北京 100074）

張現(xiàn)東

西北工業(yè)大學航空宇航制造工程專業(yè)博士研究生，主要研究方向為CAD/CAM、精密測量、熔模精密鑄造模具設(shè)計與優(yōu)化。

空心渦輪葉片成形精度低，一直是困擾我國航空制造業(yè)的難題。與國外高達80%的合格率相比，國內(nèi)的空心渦輪葉片精密鑄造合格率不到40%[1]，差距很大。統(tǒng)計結(jié)果表明：在渦輪葉片無余量熔模精密鑄造工藝過程中，同批次不合格空心葉片鑄件中由于壁厚尺寸漂移引起的形狀超差占50%左右[2]，是造成我國空心渦輪葉片合格率偏低的主要原因之一。渦輪葉片的壁厚是由外型與內(nèi)腔共同保證的，內(nèi)腔是精鑄過程中由陶芯作為轉(zhuǎn)接件形成的，因此陶芯的彎扭變形直接關(guān)系到空心渦輪葉片的壁厚分布情況。精確計算陶芯的彎扭變形，對于指導陶芯在蠟?zāi)Ｐ颓恢械亩ㄎ?、控制蠟型壁厚、提高空心渦輪葉片成形精度都有重要的作用。

葉片類截面形狀的彎扭變形一般認為是相對于理論截面的剛體位移，可以分解為相對于理論截面的平移與扭轉(zhuǎn)。目前，國內(nèi)外對葉片彎扭分析有大量的研究，王仲奇等[3]研究了彎扭葉片對降低能量損失的作用，并對彎扭葉片級和級組用于壓氣機做了大量理論研究；Houbolt等[4]將轉(zhuǎn)子葉片的彎扭變形分為翼面向彎曲、翼弦向彎曲和扭曲變形；Mohaghegh等[5]通過定義葉片積疊軸，給出了基于積疊軸的葉片彎曲變形的定義；Omer[6]研究了葉片輪廓度誤差、彎曲度誤差和扭曲度誤差定義及測量方法；陳福興等[7]研究了基于三坐標測量數(shù)據(jù)的葉片彎扭變形誤差計算方法；王紅霞[8]研究了基于位移場的葉片彎扭變形分析；彭志光[9]研究了基于凸包算法的葉片彎扭變形計算方法。上述所有方法都是針對葉片類零件提出的，并不適用于結(jié)構(gòu)更加復雜的陶芯。陶芯的制備及成形工藝在各國是保密的，對陶芯成形精度的研究更是很少見報道，陶芯彎扭變形的研究僅有黃勝利[10]、DOU等[11]提出的基于主動輪廓模型的陶芯彎扭變形計算方法，該方法通過Snake模型逼近陶芯外輪廓來獲得陶芯截面外輪廓線，但Snake模型在陶芯縱向肋間易收斂成直線，導致提取的輪廓線光順度不好，從而影響彎扭變形的計算精度。本文根據(jù)彎曲變形和扭曲變形的定義，研究和比較了幾種彎扭變形算法，并通過幾組截面輪廓線計算分析了陶芯的彎曲變形和扭曲變形。

彎曲變形和扭轉(zhuǎn)變形的定義

二維平面內(nèi)的彎曲變形和扭轉(zhuǎn)變形分別表示為測量截面相對于理論截面的平移位移和扭轉(zhuǎn)角位移。如圖1所示，實線表示理論CAD截面輪廓線，虛線表示測量截面線。圖1（a）中點Mc為測量截面的形心，點Dc為CAD截面的形心即為截面彎曲變形矢量，用Tc來表示。Tc可以分解為沿X軸向的位移Tx和沿Y軸向的位移Ty。圖1（b）中，點A和點B分別為CAD截面的后緣點和測量截面后緣點（中弧線延長線與后緣輪廓線的交點），直線McB與DcA的夾角θ表示測量截面相對于CAD截面的扭轉(zhuǎn)變形，θ逆時針為正，順時針為負。

圖1 彎扭變形定義Fig.1 Definition of bending and twist deformation

彎扭變形的計算方法

空心渦輪葉片陶芯包含眾多細小結(jié)構(gòu)：數(shù)個縱向肋槽，數(shù)十個橫向肋槽，近百個擾流柱孔、排氣邊、定位窗等，截面輪廓復雜，增加了彎扭變形的計算難度。本文根據(jù)陶芯的結(jié)構(gòu)特點，將陶芯彎扭變形的計算方法分為兩類：基于輪廓線的算法和基于數(shù)據(jù)點的算法?；谳喞€的算法是通過構(gòu)造實測點云截面輪廓線，并與CAD截面輪廓線對比，根據(jù)彎扭變形定義計算彎扭變形，主要有基于主動輪廓模型的算法和基于樣條擬合的算法等。基于數(shù)據(jù)點的算法直接通過實測點云與CAD截面輪廓線計算彎扭變形，主要有凸包算法、二維配準算法、距離權(quán)值法等。

1 基于輪廓線的算法

1.1 主動輪廓模型提取輪廓線

主動輪廓模型也稱為Snake模型,是由Kass等[12]在1987年提出的,廣泛應(yīng)用于數(shù)字圖像分析和計算機視覺領(lǐng)域。點云數(shù)據(jù)中的Snake可以用離散點pi和連接鄰點的直線段li表示：

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

基于Snake的二維點云數(shù)據(jù)輪廓線提取主要包括3部分：初始輪廓點的選取、能量方程的表示、輪廓點運動方向的確定。Snake提取陶芯截面輪廓線的步驟為：

（1）提取STL格式渦輪葉片外型截面輪廓點，并采用Kochanek曲線插值加密，以加密后的點作為Snake初始輪廓點。

（2）建立帶權(quán)值的Snake能量方程，

其中，Etension和Ebend分別表示由內(nèi)力（包括拉伸力和彎曲力）作用引起的拉伸能和彎曲能，Edist表示點云自身的特征能。

其中，θi是li-1與li的夾角，其中pi是初始輪廓上的第i點，mj是點云數(shù)據(jù)中的第j點，N1、N2分別代表初始輪廓點個數(shù)和點云數(shù)據(jù)點個數(shù)，|pimj|代表點pi到mj的距離。

（3）以渦輪葉片截面線在初始輪廓點處內(nèi)法矢方向作為該輪廓點的運動方向。

（4）如果所有輪廓點的彎曲能都小于給定閾值，即max{Edist(i)}＜ε，則輸出Snake；否則，增加初始輪廓點數(shù)，重復步驟（1） ～（4），直至收斂。

1.2 樣條擬合輪廓線

樣條擬合陶芯輪廓線包含4步：截面數(shù)據(jù)的獲取；截面輪廓數(shù)據(jù)的提?。恍兔孑喞€擬合；后緣擬合。

（1） 截面數(shù)據(jù)提?。和ㄟ^ATOS或ICT獲得陶芯表面三維點云數(shù)據(jù)，將三維點云數(shù)據(jù)與陶芯CAD模型配準，采用切片法獲得陶芯截面數(shù)據(jù)（圖 2（a））；

（2） 截面輪廓數(shù)據(jù)提?。河嬎憬孛鏀?shù)據(jù)點到CAD截面輪廓線段（圖2（b））的距離，給定距離閾值，提取截面輪廓數(shù)據(jù)（圖2（c））；

（3）截面輪廓線擬合：對截面數(shù)據(jù)排序，然后以NUBRS擬合陶芯型面輪廓線；

（4）后緣擬合：由于排氣邊的影響，陶芯后緣測量數(shù)據(jù)較少（圖3（a）），直接采用樣條擬合偏差較大（圖3（b）），理論上將陶芯后緣近似成圓弧，可采用最小二乘圓的擬合陶芯后緣輪廓線（圖3（c））。

如圖4所示，AO'和BO'分別為陶芯葉背延拓線和葉盆延拓線，內(nèi)切圓與AO'和BO'分別切于點C和點D，O'C的長度為a，e1和e2分別為O'A和O'B的單位方向法矢。則內(nèi)

圖2 截面外輪廓數(shù)據(jù)點提取Fig.2 Abstraction of measured contour points for cross-section

圖3 最小二乘法擬合陶芯后緣輪廓線Fig.3 Trailing edge fitting by least-squares method

圖4 葉盆葉背延拓線內(nèi)切圓Fig.4 Inscribed circle of blade basin and back extension cord

切圓方程可表示為：

將后緣點坐標代入內(nèi)切圓方程，就可得到后緣點到內(nèi)切圓的代數(shù)距離pi（a）：

后緣數(shù)據(jù)點的坐標變換到坐標系X'O'Y'下，帶入式（7），記方程的>值為p1（a），p2（a），p3（a），…，以min{p1（a）+p2（a）+p3（a）+…}為目標函數(shù)，求取滿足最優(yōu)解的a，即得到距所有后緣點距離和最小的圓。最后，將所得圓心坐標變換到陶芯模型坐標系中，得到陶芯后緣的擬合輪廓線。

之后，采用等半徑法計算中弧線及后緣點，采用等弧長或等參數(shù)的方式采樣輪廓線，采樣點坐標均值作為形心。采用同樣的方法得到CAD截面線的后緣點及形心，按照彎曲和扭曲變形的定義計算陶芯的彎扭變形。

2 基于數(shù)據(jù)點的計算方法

2.1 凸包算法

文獻[9]采用凸包算法求解葉片測量數(shù)據(jù)形心，首先通過凸包計算將葉片測量數(shù)據(jù)點分為凸包點和凹弧點，然后采用分割三角形法分別計算凸包點和凹弧點的形心，最后通過面積加權(quán)計算出葉片截面測量數(shù)據(jù)點的形心。由于陶芯縱向肋結(jié)構(gòu)會產(chǎn)生肋面數(shù)據(jù)點，可采用1.2節(jié)中提取截面輪廓點的方法去除內(nèi)部數(shù)據(jù)點，然后計算凸包。凸包算法計算陶芯測量數(shù)據(jù)形心的步驟為：（1）提取截面輪廓數(shù)據(jù)（圖 2（c））；（2）計算凸包（圖 5（a））；（3）采用三角形分割分別計算凸包點（圖5（b））和凹弧點形心（圖 5（c））；（4）計算陶芯截面形心。

設(shè)（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分別為凸包點和凹弧點個數(shù)，點（x1，y1）和點（x'1，y'1）相同，為凸包（凹?。┢瘘c；點（xn，yn）和點（x'm，y'm）相同，為凸包（凹?。┙K點。設(shè)陶芯截面形心坐標（xc，yc）為：

圖5 陶芯測量數(shù)據(jù)凸包計算Fig.5 Convex hull of measured contour points

凸包算法的扭轉(zhuǎn)變形，以CAD截面弦線與測量截面弦線（凸包多邊形的最大邊長）夾角來表示。

2.2 二維配準算法

二維配準法直接采用SVD-ICP算法對實測截面點云與CAD截面配準，在配準的過程中記錄平移及旋轉(zhuǎn)量，累加平移量為彎曲變形量，累加旋轉(zhuǎn)量為扭轉(zhuǎn)變形，其流程如圖6所示。

圖6 二維配準算法流程Fig.6 Flowchart of 2-D registration method

2.3 距離權(quán)值算法

距離權(quán)值法與凸包算法類似，不通過擬合截面線直接計算陶芯截面測量數(shù)據(jù)形心；不同點是凸包算法將陶芯截面看作片體，用三角分割法計算形心，距離權(quán)值法將陶芯截面輪廓看作線體，以相鄰點間的距離作為權(quán)值計算形心。算法步驟如下：

（1）提取陶芯截面輪廓數(shù)據(jù)（圖2（c））；

（2）對截面輪廓數(shù)據(jù)排序；

（3）計算距離加權(quán)形心。

設(shè)排序后截面輪廓數(shù)據(jù)點為p1，p2，p3，…pn，令pn+1=p1，對點序列p1，p2，p3，…，pn，pn+1計算距離加權(quán)形心pc：

式中，di為點pi和點pi+1的距離。

以測量數(shù)據(jù)點中距離形心最遠的點近似后緣點計算扭轉(zhuǎn)變形。

仿真計算

為驗證上述算法的準確性，截取某型號陶芯CAD模型Z=50截面，均勻離散數(shù)據(jù)點，將點集沿X軸和Y軸各偏移0.5mm，繞Z軸旋轉(zhuǎn)2°，如圖7，紅色實線為CAD模型，黑色點為變換后的數(shù)據(jù)點。

變換后的數(shù)據(jù)點視為仿真測量點，采用Snake模型算法、樣條擬合輪廓線法、凸包算法、二維配準算法、距離權(quán)值算法，進行彎扭變形計算分析。

算例1：對仿真數(shù)據(jù)不添加隨機誤差，計算結(jié)果如表1所示。

算例2：對數(shù)據(jù)點集添加服從正態(tài)分布N（0，0.01）的隨機法向誤差，計算結(jié)果如表2所示。

算例3：對數(shù)據(jù)點集添加服從正態(tài)分布N（0，0.1）的隨機法向誤差，計算結(jié)果如表3所示。

圖7 彎扭變形仿真數(shù)據(jù)Fig.7 Simulated points with known bending and twist

通過表1、2、3看出，當測量數(shù)據(jù)不存在誤差時，基于數(shù)據(jù)點的彎扭變形算法的計算精度明顯高于基于輪廓線的算法精度；當測量數(shù)據(jù)誤差較小時，5種算法的誤差均不大，樣條擬合算法和二維配準算法計算精度略高于其他3種算法；當測量數(shù)據(jù)誤差較大時，Snake模型、二維配準、凸包算法、距離權(quán)值算法的最大相對誤差均超過或接近10%，而樣條擬合算法的相對誤差低于這4種算法，對測量數(shù)據(jù)點誤差具有更好的計算穩(wěn)定性。

比較Snake模型算法與樣條擬合算法，樣條擬合算法對彎曲和扭曲變形的計算精度均高于Snake模型算法，因為Snake模型易在陶芯縱向肋間收斂成直線，導致提取輪廓線的光順度不是很好；凸包算法和距離權(quán)值算法的彎曲變形計算精度要高于扭曲變形計算精度，由于這兩種算法是通過所有（或大部分）點計算形心的，而扭曲變形是通過前后緣與弦線切點或后緣點計算的，所以計算扭曲變形對點的誤差更敏感。

表1 算例1計算結(jié)果

表2 算例2計算結(jié)果

表3 算例3計算結(jié)果

復雜陶蕊實測彎扭分析

對某型號渦輪葉片進行CT掃描點云提取7組水平截面，采用樣條擬合輪廓線的方法計算彎扭變形。水平截面位置如圖8（a）所示，截面線如圖8（b）所示，圖8（b）中實線為CAD輪廓線，虛線為擬合的測量截面線。彎扭變形計算結(jié)果見表4。

將表4中的彎扭變形計算結(jié)果繪制成彎曲變形、扭轉(zhuǎn)變形曲線，如圖9所示。

圖8 截面位置及擬合輪廓線Fig.8 Location and contour lines of cross-sections for ceramic core

表4 截面平移與扭轉(zhuǎn)結(jié)果

綜合表4及圖9可以得到陶芯的彎曲和扭曲變形的整體趨勢，發(fā)現(xiàn)存在彎扭不連續(xù)截面，說明陶芯在制坯、燒結(jié)和強化處理工藝過程中有局部變形，陶芯的模具設(shè)計和生產(chǎn)工藝還有待改進。下一步工作是研究陶芯的變形機理和基于變形調(diào)控的陶芯模具優(yōu)化設(shè)計及工藝改進技術(shù)。

結(jié)論

本文基于剛體位移的假設(shè)，給出了陶芯彎扭變形的定義；根據(jù)陶芯的結(jié)構(gòu)特點，研究了基于輪廓線和數(shù)據(jù)點的彎扭變形計算方法，通過仿真分析和試例驗證可得出：（1）當測量數(shù)據(jù)點與CAD模型完全對應(yīng)不包含誤差時，基于數(shù)據(jù)點的凸包算法、二維配準與距離權(quán)值算法對彎扭變形的計算精度要高于基于輪廓線的算法；（2）樣條擬合輪廓線法對數(shù)據(jù)點誤差的計算穩(wěn)定性要高于Snake模型、凸包算法、二維配準和距離權(quán)值算法；（3）樣條擬合輪廓線法對彎扭變形的計算精度要高于Snake模型算法，這是由于樣條擬合輪廓線的光順度要高于Snake模型；（4）凸包算法和距離權(quán)值算法對扭曲變形的計算精度要低于彎曲變形，這是由于兩種算法計算扭曲變形對測量數(shù)據(jù)點的誤差更敏感。

圖9 實測陶芯彎曲變形曲線Fig.9 Bending and torsion curves of measured ceramic core

樣條擬合輪廓線法能準確地計算陶芯的彎曲、扭曲變形，而且對后續(xù)陶芯變形位移場建模、輪廓度誤差分析及燒成收縮率計算都很方便，具有較大的工程實用價值。

[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]王仲奇,蘇杰先,鐘兢軍.彎扭葉片柵內(nèi)減少能量損失機理研究的新進展[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的葉片型面質(zhì)量評估[J].工具技術(shù),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]王紅霞.基于位移場的渦輪葉片變形特征分離技術(shù)研究[D].西安: 西北工業(yè)大學,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].計算機集成制造系統(tǒng),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.