?

基于本征正交分解的氣動外形設(shè)計空間重構(gòu)方法研究

劉南1，白俊強1，邱亞松1，華俊2

(1．西北工業(yè)大學航空學院，陜西西安710072; 2．中國航空研究院，北京100012)

摘要:在飛行器設(shè)計過程中為了提高優(yōu)化設(shè)計的尋優(yōu)精度，設(shè)計變量不斷增加，從而使整個過程更加復雜且大幅延長設(shè)計周期。針對這一問題，基于本征正交分解降階方法開展氣動外形設(shè)計空間重構(gòu)方面的研究。工作針對二維翼型開展，主要目標分為2個方面:①減少優(yōu)化過程中的氣動外形設(shè)計參數(shù);②提高設(shè)計空間中滿足設(shè)計約束的樣本比例。在Hicks-Henne參數(shù)化和POD重構(gòu)得到設(shè)計空間內(nèi)隨機選擇20 000個樣本發(fā)現(xiàn)，Hicks-Henne參數(shù)化空間中滿足設(shè)計約束的樣本比例不足25%，而重構(gòu)之后的空間則超過70%。因此，采用POD方法對設(shè)計空間進行重構(gòu)大大提高了樣本質(zhì)量，同時減少了優(yōu)化設(shè)計參數(shù)。以RAE2822進行厚度約束下的單目標升阻比增大優(yōu)化設(shè)計為例分別研究傳統(tǒng)的約束處理方法和設(shè)計空間重構(gòu)對優(yōu)化結(jié)果的影響。傳統(tǒng)約束處理方法中包括罰函數(shù)法和拒絕策略，優(yōu)化結(jié)果表明拒絕策略略優(yōu)于罰函數(shù)法，且無須設(shè)置懲罰權(quán)重，使用方便。對比重構(gòu)前后設(shè)計空間的優(yōu)化結(jié)果可見，2種傳統(tǒng)約束處理方法在32個Hicks-Henne參數(shù)化空間中最優(yōu)設(shè)計結(jié)果升阻比增加分別為27. 61%和28. 20%，采用POD方法重構(gòu)后的設(shè)計空間得到的升阻比提升分別為28. 20% 和30. 63%。因此，設(shè)計空間重構(gòu)前后的優(yōu)化精度基本類似，而且設(shè)計空間重構(gòu)之后優(yōu)化設(shè)計參數(shù)大大減少，設(shè)計效率得到明顯提升。

關(guān)鍵詞:優(yōu)化;參數(shù)化;效率;升阻比;本征正交分解方法;樣本質(zhì)量

隨著計算機硬件水平的提高和智能優(yōu)化算法的不斷完善，基于計算流體力學(computational fluid dynamics，CFD)的氣動外形智能優(yōu)化已成為現(xiàn)代飛機設(shè)計的重要手段。為了擴大設(shè)計空間和提高設(shè)計精度，往往需要大量的優(yōu)化設(shè)計變量。但是設(shè)計變量的增加不僅使整個優(yōu)化設(shè)計過程更加復雜，而且會導致優(yōu)化所需時間大大增加。

近年來，針對減少設(shè)計變量的參數(shù)化方法研究也逐漸引起關(guān)注。Chang等［1］基于正交函數(shù)對NACA系列的傳統(tǒng)翼型和超臨界翼型進行分析，結(jié)果表明采用10個正交函數(shù)能夠?qū)ΤR界翼型進行很好的描述; Robinson等［2］針對超臨界翼型采用Gram-Schmidt正交化方法得到一組正交基函數(shù)，僅使用前兩階正交基函數(shù)描述翼型，從而大大減少了設(shè)計變量; Toal等［3］使用本征正交分解方法(proper orthogonal decomposition，POD)在優(yōu)化過程中對較優(yōu)的設(shè)計樣本進行幾何層次上的過濾，重新構(gòu)建了設(shè)計空間，明顯縮短了設(shè)計周期，但是該方法以損失設(shè)計精度為代價，可能會導致優(yōu)化收斂到局部最優(yōu); Ghoman等［4］基于POD方法開展了外形參數(shù)減少的研究，結(jié)果表明本方法可以有效地減少設(shè)計變量數(shù)目并避開典型的缺陷(比如型函數(shù)選取的困難、計算效率低等原始參數(shù)化方法常見的問題)，但是沒有對設(shè)計空間的變換進行詳細討論。

與此同時，約束的處理對優(yōu)化設(shè)計的效率和精度都有很大影響，目前主要的約束處理方法有——懲罰函數(shù)方法、尋找可行解方法、保留解可行性的方法以及混合方法等［5］，但是這些方法在實際應(yīng)用中存在很多難題，比如傳統(tǒng)的罰函數(shù)方法中的懲罰系數(shù)難以確定、尋找可行解的方法很多僅適用于某一類問題等。目前在氣動外形設(shè)計中多采用罰函數(shù)一類的方法來處理約束，但是如果設(shè)計空間中不滿足約束的樣本比例很大，就會造成計算資源的浪費和設(shè)計精度的下降。

綜上所述，本文基于POD方法開展了氣動外形優(yōu)化設(shè)計過程中的設(shè)計空間降維和非法樣本點過濾方面的研究。設(shè)計算例表明，相對于不作降維和過濾處理的優(yōu)化設(shè)計，本文所提出的方法在保證尋優(yōu)精度的前提下，有效提高了優(yōu)化設(shè)計效率。

1　Hicks-Henne參數(shù)化方法

本文中對翼型的描述采用Hicks-Henne參數(shù)化方法，上下表面各采用相同數(shù)量的Hicks-Henne鼓包［6］，鼓包函數(shù)分別如下:

式中: b為鼓包高度，p為Hicks-Henne鼓包位置，x為翼型的x軸坐標點位置。本文中參數(shù)化翼型的各個鼓包高度變化范圍通過翼型所需的變形量來確定。很顯然，通過增加設(shè)計變量可以擴大優(yōu)化設(shè)計空間，從而得到較好的優(yōu)化結(jié)果，但是會延長設(shè)計周期，降低優(yōu)化效率。

2　Kriging代理模型

為了提高氣動外形優(yōu)化設(shè)計效率，本文采用代理模型代替耗時的N-S方程計算。常用的代理模型主要包括多項式響應(yīng)面模型、人工神經(jīng)元網(wǎng)格模型、徑向基函數(shù)模型及Kriging模型等，其中Kriging模型具有訓練樣本點處無偏估計、良好的高度非線性近似能力，非常適合作為代理模型使用，目前Kriging模型在工程優(yōu)化設(shè)計領(lǐng)域得到了廣泛應(yīng)用［7］。

Kriging模型將未知的函數(shù)用一個回歸函數(shù)B(x)和一個均值為零和方差為σ2的高斯隨機過程Z(x)組成，因此未知點的函數(shù)值^y為:

通過無偏估計和最大似然估計的方法，可以得到:

式中: B(x) = fT^B，f是已知點優(yōu)化變量的函數(shù)，對于常用的零階回歸函數(shù)，f是一維數(shù)組，其值為1，^B為回歸參數(shù)，通過回歸分析可以得到。Z(x) = rT(x) R－1(y－f^B)，y是已知樣本點的函數(shù)值，R是已知樣本點處的相關(guān)矩陣，r是未知點和已知樣本點之間的相關(guān)矢量，最常采用的相關(guān)函數(shù)為Gauss函數(shù)，如下所示:

式中:θ= (θ1，θ2，…，θn)T是空間相關(guān)參數(shù)矢量，可以通過最大似然估計得到θ。

3　POD降階方法

文獻［3-4］所采用的POD是一種模型降階方法，基于“快照”的思想，從大量的樣本中提取出主要特征(基模態(tài))［8］。

POD降階方法的主要思想是尋找一個子空間使所有“快照”{ Yi，i = 1，2，…，N} (Yi∈Rn)在該空間中的投影誤差最小。若向模態(tài)空間的正交投影關(guān)系定義為Φ: Rn→Rn，即

(5)式等價于尋找各階基模態(tài){φj，j = 1，2，…，r}使

對于快照集中的所有快照均成立，其中〈·，·〉是兩向量的內(nèi)積。本文選擇常用的實數(shù)域內(nèi)歐式空間中的內(nèi)積和范數(shù)，具體表達式如下:

計算POD基模態(tài)共有2種方法:①傳統(tǒng)的特征值分解降階方法，首先構(gòu)造自相關(guān)矩陣并對其僅需特征分解求得各個特征值及其對應(yīng)的特征向量，通過特征向量和各個“快照”即可求得各階基模態(tài)，而其所對應(yīng)的特征值大小表征各階基模態(tài)所含“能量”;②奇異值分解(SVD)方法計算各階基模態(tài)。SVD方法不僅計算效率高于特征值分解POD降階方法，而且在計算高階模態(tài)時比特征值分解方法更為精確［9］，因此本文采用SVD方法，具體過程如下:

假設(shè)“快照”集為{ Yi，i = 1，2，…，N}，首先得到所有快照的均值，如下:

通過快照的脈動Y'i= Yi－珘Y構(gòu)造矩陣A如下:

將A矩陣進行SVD分解得到:

式中，U∈RN×N，∑∈RN×n，V∈Rn×n，∑僅有對角線元素∑ii=σi，且滿足σ1≥σ2≥…≥σr≥0，其余元素全部為0。這樣就可以得到各階基模態(tài):

第i個“快照”可由所有基模態(tài)的線性組合得到:

通過“能量準則”可以選擇基模態(tài)數(shù)量以近似表達各個“快照”，比如選取p個基模態(tài)使前p個基模態(tài)的“能量”之和占所有基模態(tài)“能量”之和的99. 9%以上，其中基模態(tài)的“能量”的大小以其所對應(yīng)的奇異值平方表征，即

最終得到各個“快照”的近似解如下:

式中，p＜r。

4　二維翼型氣動外形優(yōu)化問題

采用粒子群優(yōu)化算法對RAE2822翼型進行氣動外形優(yōu)化設(shè)計，設(shè)計狀態(tài)為: Ma = 0. 729，α= 2. 31°，Re = 6. 5×106。為了保證CFD計算的可靠性，計算結(jié)果和風洞試驗結(jié)果的壓力分布對比如圖1所示，計算結(jié)果對前緣吸力峰值處的小凸起、上表面壓力平臺區(qū)、壓力恢復區(qū)等吻合較好，對激波的捕捉有所欠缺，基本達到工程所需的精度要求。

設(shè)計目標和約束如下:

圖1　CFD計算結(jié)果與風洞試驗結(jié)果對比

式中，Cd為阻力系數(shù)，Cl為升力系數(shù)，目標函數(shù)為升阻比的倒數(shù)Cd/Cl。thickness constraints為:前梁位于0. 16c處，后梁位于0. 6c處，前梁厚度、后梁厚度和最大厚度均不小于RAE2822翼型的相應(yīng)厚度。對于厚度約束，分別采用懲罰策略、拒絕策略和基于POD方法的設(shè)計空間過濾3種措施進行處理。

4. 1懲罰措施

懲罰措施是對約束進行處理的最一般的方式，是通過對不可行解的懲罰來將約束問題轉(zhuǎn)化為無約束問題。任何對于約束的違反都要在目標函數(shù)中添加懲罰項。但是在實際應(yīng)用中很難確定懲罰函數(shù)的形式和懲罰權(quán)重，本文中的適應(yīng)值函數(shù)如下:

式中:θ1，θ2，θ3為翼型的前梁厚度、后梁厚度和最大厚度，θo1，θo2，θo3為RAE2822翼型相應(yīng)位置處的厚度。a1，a2，a3為相應(yīng)的懲罰權(quán)重，為了保證懲罰強度，本文采用懲罰權(quán)重隨優(yōu)化代數(shù)線性遞增的方法，a1，a2，a3的初始值均設(shè)置為0. 2，最后一代的懲罰權(quán)重為100。

4. 2拒絕策略

拒絕策略也稱為死亡懲罰，其做法是直接拒絕優(yōu)化過程中所有不可行解，減小了搜索范圍。由于無需設(shè)置格外的參數(shù)，所以拒絕策略是處理約束最簡單的方法。在種群初始化中，為了保證初始種群全部滿足約束，首先采用拉丁超立方取樣方法(Latin hypercube sampling，LHS)選取大量樣本，將其中滿足約束的樣本作為初始種群;優(yōu)化過程中生成新的樣本之后首先判斷其是否滿足約束，若滿足則加入種群中，否則重新生成新的樣本，直到滿足約束為止。

4. 3基于POD方法的設(shè)計空間過濾策略

為提高設(shè)計空間中滿足約束的樣本比例，對二維翼型從Hicks-Henne參數(shù)化空間向POD模態(tài)空間的轉(zhuǎn)換進行研究。具體過程如下:

1)采用LHS方法在Hicks-Henne參數(shù)空間取出大量的樣本，并判斷樣本是否滿足厚度約束，滿足則放入“快照”集中，不滿足則舍棄該樣本;

2)采用第2節(jié)中介紹的降階方法從“快照”集中提取出各階基模態(tài)，并求得所有“快照”在各階基模態(tài)上的投影系數(shù);

3)從投影系數(shù)中提取出所有“快照”在各階基模態(tài)上投影系數(shù)的上下限，得到快照在模態(tài)空間上的投影系數(shù)的變化范圍。

這樣就把優(yōu)化問題從Hicks-Henne參數(shù)化空間轉(zhuǎn)換到模態(tài)空間。具體步驟如下:首先采用Hicks-Henne參數(shù)化方法和拉丁超立方取樣方法［10］，上下表面各選16個鼓包，共有32個參數(shù)，得到20 000個翼型，其中滿足幾何約束的約占24%左右，對這些翼型進行SVD分解得到各階基模態(tài)及其對應(yīng)的奇異值。圖2是得到的各階基模態(tài)的“能量”的變化情況，可見“能量”主要集中在前若干個基模態(tài)中。

圖2　POD各階模態(tài)的“能量”大小

所選取基向量的“能量”之和超過所有基向量“能量”之和的99. 9%，最終確定的基向量數(shù)目為16。

多次采用LHS方法在32維的Hicks-Henne參數(shù)化空間和16維的POD基模態(tài)空間分別獲得隨機不同的20 000個翼型樣本發(fā)現(xiàn)，結(jié)果如圖3所示。由圖可見，Hicks-Henne參數(shù)空間中滿足厚度約束的樣本比例約占24%，而在POD基模態(tài)空間中該比例則超過70%。由此可見，從參數(shù)化空間轉(zhuǎn)換到基模態(tài)空間之后，設(shè)計空間有所減小，許多非法的設(shè)計空間被剔除，有利于縮短優(yōu)化耗時。

圖3　不同設(shè)計空間中取樣結(jié)果對比

5　不同設(shè)計空間優(yōu)化結(jié)果對比

為了體現(xiàn)POD空間重構(gòu)方法在優(yōu)化設(shè)計中的可行性，本文采用32個Hicks-Henne鼓包構(gòu)成的參數(shù)化空間和前16個POD基向量對應(yīng)的設(shè)計空間分別對RAE2822翼型在厚度約束下的升阻比進行優(yōu)化。

5. 1訓練Kriging代理模型

在氣動外形優(yōu)化設(shè)計之前，必須首先訓練具有高可信度的Kriging代理模型，從而減少優(yōu)化設(shè)計耗時，提高效率。在Hicks-Henne參數(shù)化空間和POD基模態(tài)空間中分別采用LHS方法得到500個和300個訓練樣本，以此訓練Kriging代理模型。

然后從兩個空間中隨機抽取50個測試樣本，使用已建立的代理模型對測試樣本的升力和阻力系數(shù)進行預測。在Hicks-Henne參數(shù)化空間和POD基模態(tài)空間中建立的代理模型對測試樣本升力系數(shù)預測均方根誤差分別為1. 25%、1. 46%，阻力系數(shù)預測誤差分別為2. 77%和3. 77%。其中升力系數(shù)的計算結(jié)果和預測結(jié)果對比如圖4所示。

在優(yōu)化過程中適當?shù)貙Υ砟Ｐ瓦M行更新以保證其精度和可靠性，本文采用的方法是每代選取3個適應(yīng)值最高的樣本，采用CFD計算得到其氣動力系數(shù)，并將其加入代理模型樣本庫，對其重新進行訓練，得到新的Kriging代理模型。

圖4　計算和預測的結(jié)果對比

5. 2優(yōu)化結(jié)果

為了對比翼型優(yōu)化設(shè)計在Hicks-Henne參數(shù)化空間和POD基模態(tài)空間中的不同表現(xiàn)，采用粒子群智能優(yōu)化算法對翼型進行優(yōu)化，優(yōu)化代數(shù)為50代，每代種群由100個樣本構(gòu)成。每種策略(Hicks-Henne參數(shù)化空間+罰函數(shù)、Hicks-Henne參數(shù)化空間+拒絕策略、POD基模態(tài)空間+罰函數(shù)、POD基模態(tài)空間+拒絕策略)優(yōu)化5次，表1為每次優(yōu)化得到最優(yōu)構(gòu)型的適應(yīng)值對比。

表1　優(yōu)化結(jié)果對比

由表1可見，16個POD基向量系數(shù)的優(yōu)化結(jié)果略優(yōu)于32個Hicks-Henne參數(shù)化方法的優(yōu)化結(jié)果，而拒絕策略的結(jié)果也略好于罰函數(shù)方法。其中每種策略的最優(yōu)收斂過程對比如圖5所示。

采用12核Intel(R) Xeon(R) CPU X5650 @ 2. 67GHz圖站開展優(yōu)化設(shè)計，對比采用拒絕策略的Hicks-Henne方法和POD方法的優(yōu)化耗時如表2所示，表中時間單位為h。優(yōu)化過程中的優(yōu)化耗時主要用于Kriging代理模型的更新。更新代理模型過程中需要大量調(diào)用BLAS庫函數(shù)中的向量運算程序，較為耗時。相比于Hicks-Henne方法，采用POD方法節(jié)省了約1/3的優(yōu)化耗時。

圖5　收斂結(jié)果對比

表2　取樣與優(yōu)化耗時對比

最終優(yōu)化結(jié)果的幾何和壓力分布對比如圖6和圖7，升阻力系數(shù)見表3。4個構(gòu)型的升阻比相差在2. 0以內(nèi)，前緣吸力峰值壓力系數(shù)差量在0. 1以內(nèi)，激波均有了明顯減弱。其中采用Hicks-Henne方法優(yōu)化結(jié)果在翼型上表面70%至90%弦長位置有明顯的凸起，該凸起能夠降低逆壓梯度，采用罰函數(shù)策略的最為明顯，該優(yōu)化翼型上表面基本無激波;采用POD方法得到的2個最優(yōu)構(gòu)型均有明顯的前加載特征;而4個優(yōu)化結(jié)果在翼型下表面后緣處的后加載強度相差不大。

表3　最優(yōu)結(jié)果氣動力系數(shù)對比

圖6　優(yōu)化前后幾何外形對比

圖7　優(yōu)化前后壓力分布對比

6　結(jié)論

1)對于厚度約束情況，Hicks-Henne參數(shù)化方法得到的設(shè)計空間中滿足約束的樣本約占24. 36%，而基于POD方法重構(gòu)的設(shè)計空間中滿足約束的樣本約占70. 26%，大大提高了合法樣本的比例;

2)在基本不降低設(shè)計精度的前提下，基于本征正交分解方法的氣動外形優(yōu)化設(shè)計空間重構(gòu)能夠有效地減少了氣動外形優(yōu)化設(shè)計參數(shù)，縮短設(shè)計周期，提高尋優(yōu)效率;

3)目前最常用的2種懲罰策略——罰函數(shù)法和拒絕策略，罰函數(shù)方法要反復試驗方能得到較好的懲罰因子，從而提高設(shè)計精度，但是拒絕策略較為簡單。在厚度約束下RAE2822翼型單目標增加升阻比優(yōu)化中拒絕策略得到的優(yōu)化結(jié)果略優(yōu)于罰函數(shù)法。

參考文獻:

［1］Chang I C，Torres F J，Tung C．Geometric Analysis of Wing Sections［R］．NASA Technical Memorandum 110346，1995

［2］Robinson G M，Keane A J．Concise Orthogonal Representation of Supercritical Airfoils［J］．Journal of Aircraft，2001，38(3) : 580-583

［3］Toal D J J，Bressloff N W，Keane A J，et al．Geometric Filtration Using Proper Orthogonal Decomposition for Aerodynamic Design Optimization［J］．AIAA Journal，2010，48(5) : 918-928

［4］Ghoman S S，Wang Z，Chen P C，et al．A POD-Based Reduced Order Design Scheme for Shape Optimization of Air Vehicles ［R］．AIAA-2012-1808

［5］Coello C A C．Theoretical and Numerical Constraint-Handling Techniques Used with Evolutionary Algorithms: A Survey of the State of the Art［J］．Computer Methods in Applied Mechanics and Engineering，2002，191(11/12) : 1245-1287

［6］Hicks R，Henne P．Wing Design by Numerical Optimization［J］．Journal of Aircraft，1978，15(7) : 407-413

［7］Shao T F，Krishnamurty S，Wilmes G C．Preference-Based Surrogated Modeling in Engineering Design［J］．AIAA Journal，2007，45(11) : 2688-2701

［8］Paksoy A，Apacoglu B，Aradag S．Reduced Order Modeling of Turbulent Two Dimensional Cylinder Wake with Filtered POD and Artificial Neural Networks［R］．AIAA-2011-58

［9］Bryan R，Mohan P S，Hopkins A，et al．Statistical Modelling of the Whole Human Femur Incorporating Geometric and Material Properties［J］．Medical Engineering＆Physics，2010，32: 57-65

［10］Husslage B G M，Rennen G，van Dam E R，et al．Space-Filling Latin Hypercube Designs for Computer Experiments［J］．Optim Eng，2011，12: 611-630

Investigating Aerodynamic Shape Design Space Reconstruction Using Proper Orthogonal Decomposition (POD)

Liu Nan1，Bai Junqiang1，Qiu Yasong1，Hua Jun2

(1．College of Aeronautics，Northwestern Polytechnical University，Xi'an 710072，China 2．Chinese Aeronautical Establishment，Beijing 100012，China)

Abstract:In order to increase optimization precision in the design process of aircraft，the design parameters have augmented; this makes the design process more complicated and extends the design cycle．Aiming at this issue，reconstruction method of aerodynamic design space using Proper Orthogonal Decomposition is investigated．Main objective is divided into two aspects: one is reduction of the number of aerodynamic design parameters in design process，the other is increase of the ratio of samples which satisfy design constraints．In the Hicks-Henne parameterization space and reconstructed space by POD，20000 samples are chosen by LHS method．It is illustrated that the ratio of Hicks-Henne parameterization space satisfying thickness constraints is less than 25%，while the ratio of reconstructed space by POD is about 70%．Therefore，after space reconstruction by POD，sample quality is enhanced significantly．Moreover，the number of optimization parameters is decreased．The influence of traditional constraint handling methods and reconstructed design space is investigated with the optimization of Cl/Cd of RAE2822 airfoil．Traditional constraint handling methods include penalty function method and death penalty．It is shown that death penalty method is slightly better than penalty function method from optimization results．Additionally，there are no penalty factors in death penalty process，which is convenience to use．Comparing the optimization results before and after reconstructed design space，the improvements of Cl/Cd are 27. 61% and 28. 20% respectively for two traditional penalty handing methods in Hicks-Henne parameterization space，while the improvements are 28. 51% and 30. 63% respectively in reconstructed design space by POD．Therefore，the optimization precision is almost the same before and after space reconstruction，and the number of design parameters is decreased by half; this makes design efficiency improve significantly．

Key words:aerodynamic configurations; aircraft; angle of attack; computational fluid dynamics; design; drag coefficient; efficiency; experiments; lift drag ratio; Mach number; maximum likelihood estimation; optimization; parameterization; Reynolds number; singular value decomposition; wind tunnels; Proper Orthogonal Decomposition; sample quality

作者簡介:劉南(1989—)，西北工業(yè)大學博士研究生，主要從事降階模型、氣動優(yōu)化設(shè)計及氣動彈性力學研究。

收稿日期:2014-09-18

文章編號:1000-2758(2015) 02-0171-07

文獻標志碼:A

中圖分類號:V211．3