王鋼林 劉沛清 / WANG Ganglin LIU Peiqing

(1. 中航工業(yè)航空研究院，北京 119001； 2. 北京航空航天大學(xué)，北京 100191)

?

基于B樣條的機翼外形參數(shù)化方法研究

王鋼林1劉沛清2/ WANG Ganglin1LIU Peiqing2

(1. 中航工業(yè)航空研究院，北京 119001； 2. 北京航空航天大學(xué)，北京 100191)

描述了基于B樣條的機翼外形參數(shù)化方法在超臨界翼型和機翼氣動設(shè)計中的應(yīng)用研究。在機翼外形參數(shù)化過程中，通過B樣條曲線插值和擬合型值點生成機翼截面曲線和展向引導(dǎo)線，采用Coons-Gordon方法通過自由放樣和點陣插值獲得B樣條張量積曲面并進行疊加，進而獲得插值于機翼表面曲線網(wǎng)的B樣條曲面。還運用CFD工具對不同機翼放樣面的氣動性能進行了分析對比。通過研究，B樣條控制點可作為設(shè)計變量傳遞給氣動優(yōu)化程序，采用B樣條曲線或曲面表達式為各設(shè)計變量的線性組合，可以方便快捷獲取幾何梯度信息。

B樣條；參數(shù)化造型；Coons-Gordon；氣動優(yōu)化

0 前言

超臨界機翼設(shè)計是大型民用飛機設(shè)計的重點，精細而繁多的變量和約束被引入設(shè)計工作中，給單純依賴經(jīng)驗的設(shè)計方法帶來挑戰(zhàn)，自動化的優(yōu)化設(shè)計方法成為一種自然選擇。隨著多目標多學(xué)科優(yōu)化逐漸興起，快速精確地生成翼型與機翼形狀，并實現(xiàn)參數(shù)化造型與優(yōu)化過程的整合，已成為當前飛行器多目標優(yōu)化設(shè)計中的一個重要問題。使用遺傳算法等隨機優(yōu)化方法需要計算大量樣本，若采用精度較高的CFD模擬方法則總時間成本難以接受。使用梯度優(yōu)化算法則可大幅提高搜索效率，如離散共軛方法等[1]先進方法已被應(yīng)用于各種氣動優(yōu)化設(shè)計平臺。在基于梯度的氣動優(yōu)化設(shè)計中，為了求解目標函數(shù)對設(shè)計變量的偏導(dǎo)數(shù)[2]，首先需完成氣動外形參數(shù)化并提取設(shè)計變量(如B樣條函數(shù)控制點坐標)，而后根據(jù)鏈式求導(dǎo)方法借助CFD計算、網(wǎng)格變形等模塊進行計算。

B樣條曲線和曲面具有良好的局部控制性質(zhì)和可微性。在大型飛機增升裝置設(shè)計中，B樣條曲線還可用于生成多段翼型、運動機構(gòu)整流罩外形等，相對于傳統(tǒng)分段橢圓曲線、拋物線或三次曲線以及型函數(shù)法等造型手段，B樣條方法具有高度靈活、適用面廣、易于向三維推廣等優(yōu)勢，更加適用于計算機優(yōu)化設(shè)計平臺中的自動化成型。B樣條曲線可視為NURBS曲線的特例。國外已發(fā)表研究結(jié)果顯示，使用NURBS曲線逼近任意翼型，至多需要13個控制點便可以達到流場求解器所要求的幾何逼近精度。[3]但因加入控制點的權(quán)系數(shù)這一組變量，所以NURBS曲線逼近的求解過程變得更加復(fù)雜，而基于梯度的優(yōu)化算法的優(yōu)化效率對設(shè)計變量的個數(shù)不敏感，亦不追求對初始特定外形的精確逼近。

一些研究者已將B樣條方法應(yīng)用于氣動優(yōu)化中[3-4]，但多各表一處，缺乏對氣動外形參數(shù)化過程中各種B樣條算法的梳理總結(jié)。作者基于B樣條理論，編寫程序?qū)崿F(xiàn)了氣動外形曲線和曲面的擬合、插值等算法，將曲線網(wǎng)蒙面法用于超臨界機翼初始設(shè)計外形的參數(shù)化，以獲得優(yōu)化平臺所需的設(shè)計變量和梯度信息。

1 B樣條基礎(chǔ)

簡要介紹本文用到的B樣條基礎(chǔ)知識[5，8]。設(shè)U={u0,u1,…,un+p+1}是遞增的實數(shù)序列，則以U為樣條節(jié)點向量的第i個p次B樣條基函數(shù)Ni,p(u)由下面的遞歸形式給出：

(1)

當計算中出現(xiàn)0/0型時約定其比值為0。

設(shè){Pi|0≤i≤n}?Rd為平面d=2或空間d=3內(nèi)的控制頂點序列，則p次B樣條曲線定義為：

為了統(tǒng)一記號的方便，本文中除特別注明的情況外，所有節(jié)點向量都定義在區(qū)間[0,1]上，并約定p次B樣條采用式(4)定義的節(jié)點向量：

其中兩端的0和1各出現(xiàn)p+1次，目的是使曲線插值首尾控制點P0和Pn。

將B樣條曲線以張量積形式推廣，可得B樣條曲面的表達式：

2 翼型曲線的B樣條參數(shù)化

2.1 曲線插值和最佳平方逼近

m+1個控制點Pi是未知量，通過全局曲線插值，離散的型值點信息被翻譯成了具有良好連續(xù)性的參數(shù)曲線。為保證精確描述翼型幾何外形，型值點個數(shù)一般不少于100個，因此方程組規(guī)模較大，控制點數(shù)較多，若希望用更少的控制點完成翼型的參數(shù)化，可采用最佳平方逼近的方法。

對于三次B樣條基函數(shù)而言，記為：

A=[N0,N1,…,Nn]

c=(P0,P1,…,Pn)T，

確定的

c=c*=(P0*,P1*,…,Pn*)T

圖1 目標曲線與B樣條擬合曲線(分別有13和31個控制點)

圖2 弦長參數(shù)化法和向心參數(shù)化法擬合曲線在不同控制點數(shù)下的誤差平方和對比

在一定精度范圍內(nèi)，幾何精度更高不代表氣動性能更加接近。諸如局部最大誤差、頭部曲率半徑、切線或曲率連續(xù)性等幾何參數(shù)對翼型在不同飛行條件下的氣動性能都有重要影響。使用RANS-CFD數(shù)值模擬方法計算不同翼型在巡航狀態(tài)下的升阻力系數(shù)，計算用四邊形結(jié)構(gòu)網(wǎng)格數(shù)量約4萬，圓形計算域半徑約100倍翼型弦長。計算網(wǎng)格如圖3所示，計算結(jié)果見表1。

表1 不同擬合方案與目標翼型的升阻力系數(shù)和升阻比(α=2°,Ma=0.71，Re=6.5×106)

圖4對比了部分翼型的表面壓力系數(shù)分布。目標翼型上翼面超音速區(qū)壓力系數(shù)平順，13個控制點的擬合曲線由于上翼面外形不夠飽滿，翼型上最大負壓位置和激波位置略向后移，貢獻了更多的升力，但阻力系數(shù)增大，使得升阻比明顯低于其他翼型。使用31個控制點擬合得到的翼型升阻力系數(shù)更加靠近目標翼型，但是壓力系數(shù)曲線具有局部波動。因為三次B樣條曲線是能確保整體二階導(dǎo)數(shù)連續(xù)的，所以壓力系數(shù)的波動是由于曲率的波動導(dǎo)致的。隨著控制點數(shù)增加，節(jié)點間距離靠近，抑制了局部曲率變化對擬合精度的影響，使用55個控制點后壓力系數(shù)曲線波動幅度減弱。

2.1 帶約束的曲線擬合

則問題轉(zhuǎn)化為求解帶S個插值約束條件的優(yōu)化問題。

記點集

λ=(λ0,λ1,…,λs)T

并令

分別令式(8)對Pi(i=0,1,…,n)和λj(j=0,1,…,s)的各偏導(dǎo)數(shù)等于0，將得到的線性方程組整理后得:

在原擬合問題上添加插值約束為局部形狀控制帶來方便，但需要添加與約束條件同等數(shù)量的控制點使方程組封閉。在控制點個數(shù)足夠多時，添加插值約束后整體擬合精度基本不變，但是會導(dǎo)致局部曲率波動增大。在使用31個控制點擬合的情況下，相對于無約束的擬合效果，添加約束后擬合翼型的升阻力系數(shù)與目標翼型的差別略有增加(見表1)。通過適當增加控制點，聯(lián)合應(yīng)用曲線擬合、插值以及切線約束等方法，可以更好地照顧對于局部精度的要求。在插值約束的同時可進一步增加曲線端點處的切向約束，以精確刻畫翼型尾緣厚度與切角。但如非必要，為最佳平方逼近施加過多的約束條件是得不償失的。一方面，添加約束條件會帶來計算成本的增加；另一方面，插值約束對整體逼近造成干擾，使得難以對“不純粹”擬合方法的擬合精度做出合理預(yù)期

在使用基于梯度優(yōu)化算法的進行翼型優(yōu)化設(shè)計時，在完成對初始外形的參數(shù)化后，選用B樣條曲線控制點坐標作為設(shè)計變量，此時B樣條曲線可視為關(guān)于設(shè)計變量的線性表達式，則求解幾何外形對設(shè)計變量的梯度變得非常方便。對于B樣條曲面來說情況是類似的，所以易于向三維情形推廣，這也是B樣條幾何表示方法在先進的氣動優(yōu)化研究中獲得關(guān)注和應(yīng)用的原因。

3 機翼的B樣條參數(shù)化造型

3.1 B樣條曲面自由放樣

翼型參數(shù)化是機翼參數(shù)化的基礎(chǔ)，在完成若干展向截面線的二維參數(shù)化后，下一步即需要生成機翼曲面，以在展向的相應(yīng)位置插值各個截面線，即通常說的放樣。對于CAD軟件而言，放樣方法被普遍采用，且被認為最適用于交互設(shè)計[7]。當能夠用B樣條曲面來完成放樣時，便完成了機翼的B樣條參數(shù)化。

截面線控制點是已知的，又可表示為

解出曲面控制點Pk,i，即得到機翼放樣曲面：

傳統(tǒng)飛機設(shè)計受制于圖形理論、計算機技術(shù)和制造工藝的限制，主要采用直紋曲面生成機翼外形。圖5是兩截面間的放樣。更多的截面線意味著可以選擇更高階的B樣條曲面進行放樣，但v方向B樣條基函數(shù)次數(shù)應(yīng)滿足q≤r。現(xiàn)代大型運輸機采用的超臨界機翼為滿足巡航性能需求，采用復(fù)雜的非線性曲面，且機翼尾緣存在彎折，簡單的線性放樣方法已不適用。圖6所示為一個截掉了翼梢小翼和翼根鼓包的超臨界機翼，由曲線網(wǎng)的形式呈現(xiàn)，如果選擇對其展向共七條截面線進行放樣(三次B樣條)，得到的結(jié)果如圖7所示。

圖5 兩截面間的自由放樣曲面

圖6 超臨界機翼表面曲線網(wǎng)

圖7 超臨界機翼截面間的自由放樣

圖8中顯示放樣得到的機翼尾緣外側(cè)曲線向后突出，與實際平直的尾緣形狀不符。這是由于采用三次B樣條進行放樣時機翼轉(zhuǎn)折處附近截面位置的變化影響了外側(cè)的放樣形狀。解決方法有兩個：(1)在機翼轉(zhuǎn)折附近增加兩到三個截面以屏蔽掉轉(zhuǎn)折信息的展向傳遞；(2)使用曲線網(wǎng)蒙面法，添加引導(dǎo)線約束曲面的輪廓。

圖8 自由放樣機翼的尾緣彎折現(xiàn)象

3.2 曲線網(wǎng)蒙面法的實現(xiàn)

本節(jié)將曲線網(wǎng)蒙面算法[8]應(yīng)用于帶引導(dǎo)線的曲面放樣，完成機翼外形的參數(shù)化表達，以便于向優(yōu)化平臺傳遞設(shè)計變量和變形的梯度信息。圖9是基于某型飛機翼型和機翼平面形狀初始設(shè)計由CAD軟件導(dǎo)出的雙向曲線網(wǎng)，由r+1=7條截面線(u曲線于弦向)和s+1=11條引導(dǎo)線(v曲線于展向)組成。分別記：

為表示截面線和引導(dǎo)線的B樣條曲線向量，其中每個參數(shù)式通過最佳平方逼近得到，且滿足插值曲線交點和端點切向約束條件，算法已在上文給出。

需指出的是，對所有截面線(或引導(dǎo)線)進行擬合時，要統(tǒng)一以曲線上交點為界分段進行型值點位置的參數(shù)化，以使任意一條截面線(或引導(dǎo)線)在與同一條引導(dǎo)線(或截面線)的交點位置的參數(shù)值相等，本文稱之為“參數(shù)一致化規(guī)定”。

使用三次B樣條得到的結(jié)果已經(jīng)在圖7給出。為便于觀察，圖中選擇各個曲面上三條截面線和五條引導(dǎo)線位置的曲線進行顯示，以對比各個曲面對雙向曲線網(wǎng)和控制點陣的插值情況(原曲線網(wǎng)用虛線表示，實線為曲面在相應(yīng)位置的截線，實線覆蓋虛線即表示曲面插值曲線)。圖9(a)中，曲面L1插值所有截面線，而由于放樣時未考慮引導(dǎo)線的存在，故曲面L1在引導(dǎo)線位置存在偏差。類似地有所有引導(dǎo)線生成的自由放樣曲面，如圖9(b)所示，曲面插值所有引導(dǎo)線，由于未考慮截面線，所以在截面線位置存在較大偏差。

圖9 每個方向的放樣面、交點陣插值面和疊加后的曲面對部分引導(dǎo)線和截面線的插值情況

下一步，構(gòu)造曲面T(u,v)以插值所有的曲線網(wǎng)交點Ql,k。不妨先構(gòu)造r+1條曲線分別插值每截面線上的各s+1個交點，此曲線插值問題可作為式(6)的一個應(yīng)用加以解決。而后依式(10)～式(12)所述方法得到一個放樣面，此曲面插值于曲線網(wǎng)交點陣，但不插值各條曲線(圖9(c))

L1(u,v)、L2(u,v)和T(u,v)是定義在不同的節(jié)點向量上，通過節(jié)點加密使得L1(u,v)、L2(u,v)和T(u,v)在B樣條意義下是兼容的，即定義在相同的節(jié)點向量上，且有相同規(guī)模的控制點陣(有關(guān)節(jié)點插入和加密的實現(xiàn)方法，參見文獻[8])，將以上三個曲面做布爾運算，用兩個放樣曲面之和減去一個插值點陣曲面，所得曲面即為要求的曲線網(wǎng)插值面S(u,v)，即，

S(u,v)的控制點Rj,i可以由下式來計算：

S(u,v)被稱作Gordon曲面，本節(jié)實現(xiàn)的曲線網(wǎng)蒙面法稱為Coons-Gordon造型方法[7-8]。最終，曲線網(wǎng)插值面以B樣條參數(shù)曲面的形式存在，其控制點可作為設(shè)計變量傳遞給優(yōu)化平臺，其空間坐標可供后續(xù)模塊完成CFD前處理，B樣條曲面的導(dǎo)數(shù)可作為外形梯度信息供優(yōu)化算法調(diào)用。

3.3 不同放樣面對比分析

用上文描述的方法分別在圖6所示截面間完成自由放樣、帶前后緣引導(dǎo)線(共3條)放樣和帶全部11條引導(dǎo)線放樣，以生成B樣條曲面表示的機翼外形曲面。本節(jié)采用CFD數(shù)值分析方法，模擬上述三個機翼在巡航狀態(tài)下的流場，對比不同引導(dǎo)線放樣面的幾何外形和氣動力系數(shù)。計算用六邊形結(jié)構(gòu)網(wǎng)格數(shù)量約250萬，如圖10所示。計算邊界與機翼距離：弦向大于40倍弦長，展向大于10倍展長，上下大于40倍弦長。來流條件：Ma=0.785，Re= 2.5×107。

圖10 機翼附近的網(wǎng)格

表2是計算得到的三個機翼放樣面在1°、2°和3°時的升阻力系數(shù)和升阻比，氣動力系數(shù)參考平均氣動弦長為4.254m，參考面積63m2。自由放樣得到的機翼在各迎角下的升力系數(shù)和阻力系數(shù)均比其他兩個機翼高約1%～2%。3條引導(dǎo)線放樣得到的機翼升力系數(shù)最小，升阻比最小。11條引導(dǎo)線放樣得到的機翼在各個迎角下具有最小的阻力系數(shù)和最大的升阻比。

自由放樣得到的機翼平面形狀與帶引導(dǎo)線放樣得到的外形略有差異，主要表現(xiàn)為自由放樣曲面尾緣外側(cè)位置存在局部外凸，如圖11所示(亦如圖8所示)。因此，考慮到自由放樣機翼實際上具有更大的機翼面積，根據(jù)表2可以認為11條引導(dǎo)線放樣

表2 三個機翼放樣面的升阻力系數(shù)和升阻比

得到的機翼巡航狀態(tài)下的氣動性能更好。為便于進一步分析，沿機翼展向(圖示外露展長為b)，各放樣截面之間，自內(nèi)而外分別在7%b、17%b、23%b、32%b、50%b、80%b共六個位置截取剖面，以對比各機翼曲面在截面位置的幾何外形和表面壓力系數(shù)。截面位置標示在圖11中。

圖11 自由放樣和帶引導(dǎo)線放樣機翼平面形狀對比

圖12～圖14為圖11所示各截面外形和表面壓力系數(shù)對比。

如圖12所示，7%b位置處，11引導(dǎo)線曲面的截線上翼面最高位置相對較低，因而上翼面壓力形態(tài)比較平緩，而下翼面最低位置附近相對突出，正壓力系數(shù)略??；自由放樣面截線最大厚度位置和激波位置略靠后，是其升力系數(shù)的貢獻因素。7%b截線位置靠近翼根，機翼展向變化梯度較大，自由放樣和不同引導(dǎo)線放樣對外形的引導(dǎo)和約束比較明顯，使得附近截面外形存在差異。

如圖13所示，在17%b、23%b、32%b位置，三個曲面的截線外形沒有明顯差別，原因是在這一段區(qū)域展向變形已經(jīng)趨于平緩，且附近布置了四個放樣截面，密集的放樣截面提供了更多的幾何約束，引導(dǎo)線對外形的改變已不明顯。另外，雖然中間的三個截面外形差別不大，但自由放樣面激波位置偏后這一現(xiàn)象在所有截面位置都存在，這可能是在展向受到外側(cè)截面激波位置偏后的影響，如圖14所示。

圖12 7%b處截面外形和表面壓力系數(shù)對比

圖13 三個截面外形和表面壓力系數(shù)對比

如圖14所示，在50%b位置，自由放樣面截面翼型厚度在下翼面略減小，弦長也小于其他兩個翼型。結(jié)合圖11中的機翼平面形狀來看，自轉(zhuǎn)折位置向后，由于提供的放樣截面少、間隔大，前文已經(jīng)指出，自由放樣曲面的尾緣發(fā)生彎折。發(fā)展至80%b位置，自由放樣面截面弦長明顯增加，翼型厚度明顯增大，翼型后部變化明顯，上翼面負壓力系數(shù)增加，下翼面正壓力系數(shù)減小。不同于翼根附近，在機翼展向外側(cè)，幾何變化呈現(xiàn)出明顯的線性特征，所給翼面引導(dǎo)線也近似直線，因此，在保持前后緣外形相同的情況下，使用11條直引導(dǎo)線與只用3條引導(dǎo)線的放樣效果幾乎相同。

引導(dǎo)線的布置直接影響到各截面型值點對應(yīng)的參數(shù)值和任意展向截面上節(jié)點參數(shù)的分布形態(tài)。反映到機翼的三維設(shè)計中，在不同展向位置上，翼型變形和扭轉(zhuǎn)可產(chǎn)生豐富的外形變化，但需要對放樣效果進行更精細的控制以得到最優(yōu)的氣動外形。

圖14 兩個截面外形和表面壓力系數(shù)對比

4 結(jié)論

本文根據(jù)基于梯度優(yōu)化算法氣動優(yōu)化程序?qū)缀螀?shù)化造型和求解外形梯度信息的需求，選擇具有良好局部控制性和幾何外形表達能力的B樣條方法，對B樣條曲線帶約束擬合、B樣條自由放樣法以及基于B樣條曲面的Coons-Gordon曲線網(wǎng)蒙面法等造型方法在翼型和機翼參數(shù)化造型中的應(yīng)用進行了研究。

對于一般的超臨界翼型而言，使用30～50個控制點可使得對型值點的最大擬合誤差控制在翼型弦長的10-4量級。將要求局部插值的曲線擬合問題視為帶約束條件的優(yōu)化問題加以解決，可以在不大量增加控制點情況下做到在造型時對局部形狀的嚴格控制。B樣條曲面自由放樣算法是二維參數(shù)化向三維參數(shù)化過渡的橋梁。超臨界機翼的非線性外形使其有別于傳統(tǒng)直放樣機翼，用Coons-Gordon曲線網(wǎng)蒙面法復(fù)現(xiàn)了某型飛機超臨界機翼曲面。增加機翼前后緣和翼面引導(dǎo)線可以有效約束和引導(dǎo)機翼曲面的空間輪廓和表面脈絡(luò)。

B樣條方法具有良好的控制能力和交互性能，且便于實現(xiàn)和幾何文件的信息轉(zhuǎn)換。在完成初始外形的參數(shù)化后，B樣條曲線或曲面控制點可以直接作為氣動優(yōu)化的設(shè)計變量，B樣條參數(shù)式成為設(shè)計變量的線性表達式，簡化了氣動優(yōu)化中梯度信息的傳遞和計算。B樣條方法可以運用于更復(fù)雜的氣動外形參數(shù)化造型中，例如可以向翼梢小翼、翼根鼓包、發(fā)動機短艙甚至飛翼等進行推廣，并進行氣動精度與幾何精度的綜合研究。另外，作者認為用NURBS方法完成參數(shù)化是今后最有希望的選擇之一，但是需要克服三維造型時算法的改進和實現(xiàn)問題。

[1] Nield B N. An overview of the Boeing 777 high lift aerodynamic design[J].AeronauticalJournal，1995，99(989): 361-371.

[2] Samareh J A. Survey of shape parameterization techniques for high-fidelity multidisciplinary shape optimization[J]. AIAA journal，2001，39(5): 877-884.

[3] Hatanaka K，Obayashi S，Jeong S. Application of the variable-fidelity MDO tools to a jet aircraft design[C]. The 25th International Congress of the Aeronautical Science， 2006.

[4] Jé, r-iuml，Lé，et al. Optimized nonuniform rational B-spline geometrical representationfor aerodynamic design of wings[J]. AIAA journal，2001，39(11): 2033-2041.

[5] 羅鐘鉉，孟兆良，劉成明. 計算幾何——曲面表示論及應(yīng)用[M]. 北京:科學(xué)出版社， 2010.

[6] 顏慶津. 數(shù)值分析[M]. 北京:北京航空航天大學(xué)出版社，2006.

[7] Woodward C D. Skinning techniques for interactive B-spline surface interpolation[J]. Computer-Aided Design，1988，20(8): 441-451.

[8] Piegl L A，Tiller W. The NURBS book[M]. Springer，1997.

A Method of Shape Parameterization Based on B-spline for Wing Design

(1.Chinese Aeronautical Establishment，Beijing 119001，China；2. Beijing University of Aeronautics and Astronautics， Beijing 100191， China)

The paper describes the application of B-spline parameterization method on the supercritical airfoil and wing design. Wing section curves and span wise guide lines are generated by B-spline curve fitting and interpolation. Two lofting surfaces and a points-interpolation surface were overlaid by Coons-Gordon method, and bidirectional curve network on wing surface is interpolated.The aerodynamic characteristics of different loft surfaces were compared with CFD analysis. The research results show B-spline control points are passed to aerodynamic optimization process as design variables, and the B-spline expression can be treated as a linear combination of the design variables, which provides an easy access to the geometry gradient information.

B-spline; shape parameterization; Coons-Gordon; aerodynamic optimization

10.19416/j.cnki.1674-9804.2016.03.002

V211.41

A