鄧金秋 馮仁忠
(北京航空航天大學(xué) 數(shù)學(xué)與系統(tǒng)科學(xué)學(xué)院,北京 100191)
利于翼型優(yōu)化設(shè)計(jì)的超臨界翼型參數(shù)化方法
鄧金秋 馮仁忠
(北京航空航天大學(xué) 數(shù)學(xué)與系統(tǒng)科學(xué)學(xué)院,北京 100191)
為減少超臨界翼型優(yōu)化中的設(shè)計(jì)變量,消除優(yōu)化結(jié)果的不光順現(xiàn)象、保證C2連續(xù),在優(yōu)化過(guò)程中控制翼型幾何特性的變化范圍,設(shè)計(jì)出了由4條首尾相接的有理Bézier曲線表示的超臨界翼型的翼型參數(shù)化方法,該方法對(duì)翼型數(shù)據(jù)的參數(shù)化過(guò)程中主要運(yùn)用了Bézier曲線擬合算法與SPSA(Simultaneous Perturbation Stochastic Approximation)優(yōu)化算法,并在Bézier曲線擬合算法中使用了有別于常用方法的數(shù)據(jù)點(diǎn)參數(shù)選擇方法.將這種超臨界翼型參數(shù)化方法與優(yōu)化算法結(jié)合便可實(shí)現(xiàn)翼型優(yōu)化設(shè)計(jì),其中的設(shè)計(jì)變量為21個(gè),優(yōu)化結(jié)果不僅光順且滿足C2條件,通過(guò)設(shè)定設(shè)計(jì)變量變化范圍便可控制相應(yīng)的翼型前緣半徑、上下弦最高最低點(diǎn)的位置與曲率、尾部契角等幾何特征.
翼型;曲線;參數(shù)化;優(yōu)化
超臨界翼型具有一些優(yōu)良的氣動(dòng)特性,在民航客機(jī)和大型運(yùn)輸機(jī)的設(shè)計(jì)中有著廣泛的應(yīng)用.基于氣動(dòng)性能數(shù)值計(jì)算的超臨界翼型優(yōu)化設(shè)計(jì)具有周期短、費(fèi)用低的優(yōu)點(diǎn).因此將數(shù)值模擬計(jì)算與風(fēng)洞試驗(yàn)相結(jié)合可以縮短研發(fā)周期、降低研發(fā)費(fèi)用,這已成為目前超臨界翼型研究的發(fā)展方向[1-4].
為實(shí)現(xiàn)基于數(shù)值計(jì)算的超臨界翼型優(yōu)化,就需要設(shè)計(jì)相應(yīng)的翼型參數(shù)化方法,即將由離散數(shù)據(jù)點(diǎn)表示的待優(yōu)化翼型轉(zhuǎn)化為由含參數(shù)的翼型函數(shù)表示,再選擇翼型函數(shù)的適當(dāng)參數(shù)作為設(shè)計(jì)變量,結(jié)合優(yōu)化算法與流場(chǎng)計(jì)算實(shí)現(xiàn)優(yōu)化設(shè)計(jì).可見(jiàn)翼型參數(shù)化對(duì)于翼型優(yōu)化設(shè)計(jì)來(lái)說(shuō)是十分重要的[5-7].
現(xiàn)有的翼型參數(shù)化方法主要存在如下問(wèn)題:①設(shè)計(jì)變量過(guò)多,如通過(guò)在基準(zhǔn)翼型上添加局部擾動(dòng)函數(shù)的方法會(huì)有多達(dá)數(shù)十個(gè)設(shè)計(jì)變量;②優(yōu)化結(jié)果存在不光順現(xiàn)象,如B樣條表示方法[8]與Hicks-Henne外形函數(shù)法[9];③翼型函數(shù)不滿足C2條件,如 Figures樣條法[10];④很難在優(yōu)化過(guò)程中約束翼型的特定幾何特性的變化范圍[11].為解決上述4點(diǎn)問(wèn)題,本文結(jié)合超臨界翼型的幾何特性提出了一種基于分段有理Bézier曲線的超臨界翼型參數(shù)化方法.通過(guò)這種方法得到的翼型函數(shù)共有21個(gè)參數(shù),滿足C2條件,不存在不光順現(xiàn)象,并且翼型參數(shù)直接和翼型幾何特征相對(duì)應(yīng).本文首先介紹翼型參數(shù)化的具體過(guò)程,然后結(jié)合SPSA(Simultaneous Perturbation Stochastic Approximation)優(yōu)化算法給出具體的優(yōu)化實(shí)例.
由于超臨界翼型具有較復(fù)雜的幾何形狀,很難通過(guò)一條低次Bézier曲線表示,因此需要使用分段Bézier曲線表示翼型.另一方面為了降低翼型參數(shù)的個(gè)數(shù),應(yīng)盡量減少分段數(shù)目與每段Bézier曲線的次數(shù).綜合這兩方面因素,本文提出了如下翼型參數(shù)化方法.
設(shè)翼型函數(shù)為C(t),C(t)由4條3次有理Bézier曲線C0(t)~C3(t)首尾相連拼接而成,曲線Ci(t)的4個(gè)控制頂點(diǎn)為,見(jiàn)圖1.
圖1 翼型函數(shù)C(t)的組成形式
下面以 NASA SC(2)0712[12]為例介紹翼型參數(shù)化的具體步驟.
首先將控制頂點(diǎn)權(quán)重均置為1,并對(duì)曲線控制頂點(diǎn)做如下初步設(shè)定:
C0(t)段曲線:將坐標(biāo)值設(shè)為翼型后緣上表面最后一個(gè)數(shù)據(jù)點(diǎn)的坐標(biāo)值;將坐標(biāo)值設(shè)為翼型上弧線最高點(diǎn)坐標(biāo)值,由于數(shù)據(jù)點(diǎn)足夠密,可以直接設(shè)為數(shù)據(jù)點(diǎn)中縱坐標(biāo)最大點(diǎn)坐標(biāo);將縱坐標(biāo)值設(shè)為縱坐標(biāo)值,使得與共線且水平;不對(duì)橫坐標(biāo)值與橫縱坐標(biāo)值進(jìn)行任何設(shè)定.
C1(t)段曲線:將坐標(biāo)值設(shè)為坐標(biāo)值,使得與重合;將縱坐標(biāo)值設(shè)為縱坐標(biāo)值,使得和共線且水平;將坐標(biāo)值設(shè)為翼型前緣頂端坐標(biāo)值,由于數(shù)據(jù)點(diǎn)足夠密,可以直接設(shè)為數(shù)據(jù)點(diǎn)中橫坐標(biāo)最小點(diǎn)坐標(biāo)值;將橫坐標(biāo)值設(shè)為橫坐標(biāo)值,使得與共線且豎直;不對(duì)縱坐標(biāo)值與橫坐標(biāo)值進(jìn)行任何設(shè)定.
C2(t)段曲線:將坐標(biāo)值設(shè)為坐標(biāo)值,使得與重合;將橫坐標(biāo)值設(shè)為橫坐標(biāo)值,使得和共線且豎直;將坐標(biāo)值設(shè)為翼型下弧線最低點(diǎn)坐標(biāo)值,由于數(shù)據(jù)點(diǎn)足夠密,可以直接設(shè)為數(shù)據(jù)點(diǎn)中縱坐標(biāo)最小點(diǎn)坐標(biāo);將縱坐標(biāo)值設(shè)為縱坐標(biāo)值,使得與共線且水平,不對(duì)縱坐標(biāo)值與橫坐標(biāo)值進(jìn)行任何設(shè)定.
C3(t)段曲線:將坐標(biāo)值設(shè)為坐標(biāo)值,使得與重合;將縱坐標(biāo)值設(shè)為縱坐標(biāo)值,使得和共線且水平;將坐標(biāo)值設(shè)為翼型后緣下表面最后一個(gè)數(shù)據(jù)點(diǎn)的坐標(biāo)值;不對(duì)橫坐標(biāo)值與橫縱坐標(biāo)值進(jìn)行任何設(shè)定.
通過(guò)初步設(shè)定,可以確定一部分控制頂點(diǎn)的坐標(biāo),這些坐標(biāo)在之后參數(shù)化過(guò)程中均保持不變.而對(duì)于另一部分仍處于自由狀態(tài)的控制頂點(diǎn)坐標(biāo)則需要在下一個(gè)步驟中確定.
通過(guò)對(duì)控制頂點(diǎn)坐標(biāo)的初步約束可以發(fā)現(xiàn),4條Bézier曲線的端點(diǎn)坐標(biāo)均已確定了下來(lái),因此每條曲線也自然對(duì)應(yīng)了一段翼型數(shù)據(jù)點(diǎn).下面就要確定仍處于自由狀態(tài)的控制頂點(diǎn)坐標(biāo),使每段曲線最大限度貼合到其對(duì)應(yīng)的一組數(shù)據(jù)點(diǎn)組.由于改變?nèi)我庖粭l曲線不會(huì)對(duì)其它3條曲線產(chǎn)生影響,因此可以將4段Bézier曲線分別擬合到對(duì)應(yīng)的數(shù)據(jù)點(diǎn)組上,進(jìn)而求出所有的控制頂點(diǎn)坐標(biāo),下面以C1(t)為例說(shuō)明擬合方法.對(duì)應(yīng)的數(shù)據(jù)點(diǎn)組如下:
設(shè)C1(t)對(duì)應(yīng)的數(shù)據(jù)點(diǎn)組為,用(xpj,ypj)表示 pj的坐標(biāo)值,即設(shè)1],為Bernstein基函數(shù),為C1(t)控制頂點(diǎn),用表示的坐標(biāo)值,即
由累積弦長(zhǎng)法[13]得到對(duì)應(yīng)參數(shù),并將等比例變換到區(qū)間[0,1],使.下面使用最小二乘法確定使得最小.
圖2 初步計(jì)算結(jié)果
圖3 由累積弦長(zhǎng)法得到的對(duì)應(yīng)關(guān)系
圖4 數(shù)據(jù)點(diǎn)與曲線上與其橫坐標(biāo)相等點(diǎn)的對(duì)應(yīng)關(guān)系
圖5 以為參數(shù)進(jìn)行最小二乘得到的結(jié)果
表1 誤差降低速度
至此已完全確定了控制頂點(diǎn)坐標(biāo),在之后的參數(shù)化過(guò)程中所有控制頂點(diǎn)坐標(biāo)均保持不變.
通過(guò)之前計(jì)算得到的翼型函數(shù)在Bézier曲線拼接點(diǎn)處僅滿足C1條件,其他部分滿足C2條件.拼接點(diǎn)曲率見(jiàn)表2.
表2 拼接點(diǎn)曲率
3 次有理 Bézier曲線有如下性質(zhì)[13]:
1) 設(shè)4 個(gè)控制頂點(diǎn)權(quán)重為 w0,w1,w2,w3,那么可以在不改變曲線外形的前提下將4個(gè)控制頂點(diǎn)權(quán)重調(diào)整為 1,v1,v2,1,其中 v1,v2相互獨(dú)立且由 w0,w1,w2,w3唯一確定.
2)保持控制頂點(diǎn)位置不變,設(shè)4個(gè)控制頂點(diǎn)權(quán)重為1,v1,v2,1 則 v1,v2與兩個(gè)端點(diǎn)的曲率一一對(duì)應(yīng).
由上述兩條性質(zhì)可知3次有理Bézier曲線可以由4個(gè)控制頂點(diǎn)坐標(biāo)和兩個(gè)端點(diǎn)曲率唯一確定.為使拼接點(diǎn)處的曲率連續(xù),只需固定控制頂點(diǎn)坐標(biāo),將相鄰Bézier曲線在拼接點(diǎn)處曲率設(shè)為相同值.
顯然人為給定拼接點(diǎn)處曲率,如設(shè)為相鄰Bézier曲線在拼接點(diǎn)曲率的平均值,會(huì)導(dǎo)致之前得到的翼型函數(shù)出現(xiàn)較大變化,增大翼型函數(shù)與數(shù)據(jù)點(diǎn)之間的誤差.所以不妨使用優(yōu)化算法,以翼型函數(shù)與數(shù)據(jù)點(diǎn)的誤差為優(yōu)化目標(biāo),以3個(gè)拼接點(diǎn)的曲率為設(shè)計(jì)變量,使用 SPSA[14-15]方法進(jìn)行優(yōu)化計(jì)算.SPSA方法的主要思想是通過(guò)隨機(jī)擾動(dòng)各個(gè)設(shè)計(jì)變量使目標(biāo)函數(shù)逐漸收斂到極小值,其算法如下:
Input:
SPSA 參數(shù) a,c,A,α,γ;
/* 分別設(shè)為1,0.1,80,0.6,0.7*/
目標(biāo)函數(shù) y(θ0,θ1,θ2);
/*y(θ0,θ1,θ2)為翼型函數(shù)與數(shù)據(jù)點(diǎn)的誤差,在保持控制頂點(diǎn)坐標(biāo)不變的前提下僅與3個(gè)拼接點(diǎn)的曲率有關(guān)*/
Output:
Rθ0,Rθ1,Rθ2.
/*使得誤差最小的3個(gè)拼接點(diǎn)曲率*/
Function:
for(k=0;;k++)
{
Δk=[Δk0,Δk1,Δk2]T;
//Bernoulli01(1,-1)
計(jì)算后誤差為0.002 06,最終得到的翼型函數(shù)見(jiàn)圖6.
圖6 最終得到的翼型函數(shù)
將這種參數(shù)化方法應(yīng)用于NASA SC(2)翼型庫(kù),翼型函數(shù)與數(shù)據(jù)點(diǎn)誤差見(jiàn)表3.
表3 翼型函數(shù)與數(shù)據(jù)點(diǎn)誤差
表3表明這種翼型參數(shù)化方法可以較精確地表示常見(jiàn)的超臨界翼型,因此可以為優(yōu)化提供足夠豐富的翼型搜索空間.
翼型參數(shù)應(yīng)滿足如下兩條性質(zhì):①通過(guò)翼型參數(shù)可以還原翼型函數(shù);②各參數(shù)之間應(yīng)相互獨(dú)立.由于3次有理Bézier曲線可以由4個(gè)控制頂點(diǎn)坐標(biāo)和兩個(gè)端點(diǎn)曲率唯一確定,所以應(yīng)選擇4條有理Bézier曲線的控制頂點(diǎn)和端點(diǎn)曲率組成翼型參數(shù){Pn}.但由于拼接點(diǎn)連續(xù)且曲率相等,上下兩條控制邊水平,前緣控制邊豎直導(dǎo)致了{(lán)Pn}中含有相關(guān)聯(lián)參數(shù).
因此去除相關(guān)聯(lián)參數(shù),最終得到21個(gè)相互獨(dú)立的參數(shù):參數(shù)1(縱坐標(biāo))、參數(shù)2(橫坐標(biāo))、參數(shù)3(縱坐標(biāo))、參數(shù)4(橫坐標(biāo))、參數(shù)5(縱坐標(biāo))、參數(shù) 6(橫坐標(biāo))、參數(shù)7(橫坐標(biāo))、參數(shù)8(縱坐標(biāo))、參數(shù)9(縱坐標(biāo))、參數(shù)10(橫坐標(biāo))、參數(shù)11(縱坐標(biāo))、參數(shù) 12橫坐標(biāo))、參數(shù)13(橫坐標(biāo))、參數(shù)14(橫坐標(biāo))、參數(shù)15(縱坐標(biāo))、參數(shù)16(縱坐標(biāo))、參數(shù)17(處曲率)、參數(shù)18(處曲率)、參數(shù) 19(處曲率)、參數(shù)20(處曲率)、參數(shù)21(處曲率).
目前常用的翼型參數(shù)化方法為Hick-Henne外型函數(shù)法,該方法通過(guò)大量局部擾動(dòng)函數(shù)的疊加來(lái)組合出翼型曲線.由于任意一段翼型曲線均為多個(gè)局部擾動(dòng)函數(shù)疊加的結(jié)果,因此極易出現(xiàn)凹凸不平等不光順現(xiàn)象.而通過(guò)本方法得到的翼型曲線是由3次有理Bézier曲線拼接而成,由于每段Bézier曲線的內(nèi)部均是光順的且拼接處滿足C2條件,因此得到的翼型曲線具有很好的光滑度.與Hick-Henne外型函數(shù)法相比,本方法得到的翼型曲線不會(huì)出現(xiàn)不光順現(xiàn)象,因此更有利于流場(chǎng)計(jì)算.
另一種經(jīng)典的翼型參數(shù)化方法為Figures樣條法.這種方法同樣應(yīng)用分段Bézier曲線表示翼型,但在兩段曲線的拼接點(diǎn)只滿足C1條件而不滿足C2條件,因此所得翼型函數(shù)的曲率不連續(xù).這就導(dǎo)致在流場(chǎng)計(jì)算的過(guò)程中會(huì)出現(xiàn)湍流現(xiàn)象,影響結(jié)果的準(zhǔn)確性.而本方法保證了拼接點(diǎn)的C2連續(xù)性,進(jìn)而確保了流場(chǎng)計(jì)算的準(zhǔn)確性.
此外,現(xiàn)有的絕大部分翼型參數(shù)化方法都無(wú)法將翼型參數(shù)與翼型幾何性質(zhì)緊密的聯(lián)系在一起,這就導(dǎo)致了翼型優(yōu)化的盲目性.本方法所選取的翼型參數(shù)與翼型幾何特性緊密相關(guān),因此大大增強(qiáng)了翼型優(yōu)化的目的性與針對(duì)性.
為說(shuō)明將該翼型參數(shù)化方法應(yīng)用與翼型優(yōu)化設(shè)計(jì)的可行性,本文給出了具體的優(yōu)化實(shí)例.
由于本文討論的重點(diǎn)在于構(gòu)造翼型參數(shù)化方法,所以沒(méi)有涉及流場(chǎng)計(jì)算與高效優(yōu)化算法的設(shè)計(jì).這里通過(guò)使用SPSA優(yōu)化算法將基準(zhǔn)翼型優(yōu)化為目標(biāo)翼型的方法來(lái)證明本文提出的翼型參數(shù)化方法可以應(yīng)用于翼型優(yōu)化設(shè)計(jì).優(yōu)化結(jié)果見(jiàn)表4、圖7和圖8.
表4 優(yōu)化效果
圖7所示優(yōu)化算例中的基準(zhǔn)翼型為SC(2)0503,目標(biāo)翼型為SC(2)0706.圖7a中離散點(diǎn)為目標(biāo)翼型數(shù)據(jù)點(diǎn),曲線為基準(zhǔn)翼型,圖7b中離散點(diǎn)為目標(biāo)翼型數(shù)據(jù)點(diǎn),曲線為優(yōu)化結(jié)果.優(yōu)化結(jié)果與目標(biāo)翼型誤差為0.001450.
圖8所示優(yōu)化算例中的基準(zhǔn)翼型為SC(2)0614,目標(biāo)翼型為SC(2)1010.圖8a中離散點(diǎn)為目標(biāo)翼型數(shù)據(jù)點(diǎn),曲線為基準(zhǔn)翼型,圖8b中離散點(diǎn)為目標(biāo)翼型數(shù)據(jù)點(diǎn),曲線為優(yōu)化結(jié)果.優(yōu)化結(jié)果與目標(biāo)翼型誤差為0.001218.
圖8 將基準(zhǔn)翼型SC(2)0614優(yōu)化為目標(biāo)翼型SC(2)1010
本文給出的翼型參數(shù)化方法所含設(shè)計(jì)變量較少,共21個(gè);設(shè)計(jì)的基于有理Bézier曲線的翼型函數(shù)不僅光順且滿足C2條件;參數(shù)與翼型幾何性質(zhì)對(duì)應(yīng),通過(guò)限制參數(shù)變化范圍便可控制對(duì)應(yīng)的幾何參數(shù);該方法可以為翼型優(yōu)化提供足夠豐富的翼型搜索空間.這些特性表明,該翼型參數(shù)化方法會(huì)為下一步的翼型優(yōu)化提供很大的便利.此外分段有理 Bézier曲線可以很方便地轉(zhuǎn)化為NURBS曲線,便于工程應(yīng)用.
References)
[1]劉周,朱自強(qiáng),付鴻雁,等.高升阻比翼型的設(shè)計(jì)[J].空氣動(dòng)力學(xué)學(xué)報(bào),2004,22(4):410 -414 Liu Zhou,Zhu Ziqiang,F(xiàn)u Hongyan,et al.Design of airfoil with high ratio of lift over drag[J].Acta Aerodynamica Sinica,2004,22(4):410-414(in Chinese)
[2]王曉璐,朱自強(qiáng),劉周.基于N-S方程的翼型雙設(shè)計(jì)點(diǎn)雙目標(biāo)優(yōu)化設(shè)計(jì)[J].北京航空航天大學(xué)學(xué)報(bào),2006,32(5):503-507 Wang Xiaolu,Zhu Ziqiang,Liu Zhou.Bi-point/bi-objective optimization design of ailfoil using N-Sequations[J].Journal of Beijing University of Aeronautics and Astronautics,2006,32(5):503-507(in Chinese)
[3]周濤,張淼,李亞林.基于全速勢(shì)方程的超臨界翼型設(shè)計(jì)[J].航空計(jì)算技術(shù),2009,39(4):58 -64 Zhou Tao,Zhang Miao,Li Yalin.Supercritical airfoil design based on full potential equations[J].Aeronautical Computer Technique,2009,39(4):58 -64(in Chinese)
[4] Khurana M,Winar to H.Application of swarm approach and artificial neural networks for airfoil shape optimization[R].AIAA-2008-5954,2008
[5] Heine B,Mack S.Aerodynamic scaling of general aviation airfoil for low Reynolds number application[R].AIAA-2008-4410,2008
[6] Padulo M,Maginot J.Airfoil design under uncertainty with robust geometric parameterization[R].AIAA-2009-2270,2009
[7] Haderlie J,Crossley W.A parametric approach to supercritical airfoil design optimization[R].AIAA-2009-6950,2009
[8] Lepine J,Guibault F,Trepanier JY,et al.Optimized nonuniform rational B-spline geometrical representation for aerodynamic design of wings[J].AIAA,2001,39(11):570 - 577
[9] Hicks R M,Hennet P A.Wing design by numerical optimization[J].Journal of Aircraft,1978,15(7):407 -412
[10] Sobester A,Keane A J.Airfoil design via Cubic splines fer guson's curves revisited[R].AIAA 2007-2881,2007
[11] Song Wenbin,Keane A J.A study of shape parameterization methods for airfoil optimization[R].AIAA 2003-4482,2004
[12] Harris C D.NASA supercritical airfoils a matrix of family-related airfoils[R].NASA Technical Paper 2969,1990
[13] Farin G E.Curves and surfaces for computer-aided geometric design[M].4.San Diego:Academic Press,1997:53 - 56,172-224
[14] Spall JC.Implementation of the simultaneous perturbation algorithm for stochastic optimization[J].IEEE Transactions on Aerospace and Electronic Systems,1998,34(3):817 -823
[15] Song Qing,Spall JC,Ni Jie.Robust neural network tracking controller using simultaneous perturbation stochastic approximation[J].Neural Networks IEEE Transactions,2008,19(5):817-835
(編 輯:趙海容)
Supercritical airfoil parameterization method feasible to optimum design
Deng Jinqiu Feng Renzhong
(School of Mathematics and Systems Science,Beijing University of Aeronautics and Astronautics,Beijing 100191,China)
To reduce the number of design variables in the supercritical airfoil optimization,eliminate the unfairness phenomenon,ensure C2condition,control the geometric characteristics of the airfoil in the optimization process,a parameterization method for supercritical airfoil based on four rational Bézier curve was presented.In the parametric process to the airfoil data,the Bézier curve approximation algorithm and SPSA(simultaneous perturbation stochastic approximation)algorithm were used and in the Bézier curve approximation algorithm the way to choose the parameter to the data points is different from the common method.Supercritical airfoil shape optimization can be achieved by combining this method with optimization algorithms.It contents 21 design variables,the optimization result is fair and satisfies the C2condition,the geometric characteristics of the airfoil such as leading edge radius,upper and lower crest location including curvature there,the boat tail angle can be controlled in the optimization process by setting the range of the design variables.
airfoils;curves;parameterization;optimization
V 211.3
A
1001-5965(2011)03-0368-06
2010-01-25
鄧金秋(1988 -),男,北京人,碩士生,deng.jinqiu@gmail.com.