幾何不確定性區間分析及魯棒氣動優化設計

宋鑫， 鄭冠男,*， 楊國偉， 姜倩

(1. 中國科學院力學研究所流固耦合系統力學重點實驗室, 北京 100080；2. 中國科學院大學工程科學學院, 北京 100049)

航空航天工業中，飛行器設計的各參數往往被看作確定的量，但實際上，工程中是無法避免不確定性存在的。例如飛行器結構設計由于制造水平、測量誤差等會造成參數的不確定性，飛行器飛行環境及載荷作用也會存在更強的不確定性等[1]。不確定性因素嚴重影響了與其相對的確定性設計的精度和可靠度，因此，想要獲得更可靠更穩定的設計，必須考慮不確定性因素的影響。考慮不確定性的設計方法主要包括魯棒(穩健)優化設計方法[2]及可靠性優化設計方法[3]。魯棒優化設計旨在降低飛行器性能對不確定因素的靈敏度，而可靠性設計旨在降低發生故障或失效的概率，提高飛行器的可靠性。本文主要研究翼型的魯棒優化設計方法。

與傳統的確定性優化相比，魯棒優化設計需要在優化過程中不斷進行不確定性分析。其主要包括概率方法和非概率方法。概率方法又包括蒙特卡羅(Monte Carlo)方法、泰勒展開、隨機展開等，如徐明等[4]采用蒙特卡羅方法對旋翼的最優轉速進行了不確定性分析；Salehi等[5]和鄔曉敬等[6]采用非浸入式混沌多項式方法對離心泵的流場及翼型跨聲速隨機氣動特性等分別進行了不確定性分析；戴玉婷和楊超[7]考慮廣義剛度的隨機不確定性，基于浸入式隨機展開方法分析了不確定性對顫振邊界的影響；Papadimitriou等[8]考慮流動相關參數及幾何不確定性，并采用離散網格技術進行不確定性分析，對阻力系數目標進行了魯棒優化設計。概率方法一般要事先給定不確定性變量的概率分布函數，但是工程中很難獲得充足的樣本來估計分布函數的參數或驗證概率分布函數的合理性，而人為假設又會導致較大的分析誤差[9]，相比之下獲得不確定性因素的取值范圍較容易，因此研究人員發展了“基于邊界表征”的非概率不確定性方法，其中區間分析方法應用較多。屈小章等[10]采用區間方法描述風機翼型的幾何不確定性變量，對其氣動性能進行了穩健優化設計；Zheng和Qiu[11]采用區間分析的方法對考慮不確定性的氣動力和氣動熱進行了分析；張軍紅和韓景龍[12]提出了一種考慮區間不確定性的機翼顫振優化方法。

飛行器生產時本身存在制造誤差，在飛行過程中也會由于載荷的作用以及結冰、燒蝕等現象使飛行器氣動外形產生局部變化，進而影響氣動性能。另外，飛行器跨聲速飛行時，流場具有強非線性，幾何外形變化對氣動性能的影響會更大、更復雜。以往研究中，對幾何不確定性變量或采用假設的概率方法描述[8,13]，容易造成較大的分析誤差；或采用基于特征幾何參數的區間描述[10-11]，無法表達翼型任意的局部變化。且大多數非概率方法研究均基于區間數運算分析方法，該方法適用于分析函數在不確定區間內為單調函數的情況，不適用于復雜的非線性氣動問題。

針對以上問題，本文研究了幾何不確定性的通用區間描述及快速非線性區間分析方法。在此基礎上建立了魯棒優化設計流程，并通過標準算例驗證了方法的有效性。

1 區間不確定性分析

1.1 直接操作自由變形參數化方法

采用區間方法進行幾何不確定性分析，首先要通過合適的參數化方法建立翼型幾何不確定性的區間描述。工程中現已發展出多種翼型參數化方法[14]，如類別形狀函數變換(Class Shape Transformation, CST)方法、Hicks-Henne型函數法、PARSEC方法、B樣條法等。上述方法或通過間接變量描述翼型，或通過如前緣半徑、最大厚度位置等特征參數描述，但均無法直接控制翼型上某點的變化，不能靈活表達翼型表面凹凸等局部變形，故不適合描述任意的幾何不確定性。本文采用直接操作自由變形(DFFD)方法[15]對翼型參數化，以翼型表面點的直接變形為參數建立幾何不確定性量區間描述。

對于一般自由變形(FFD)方法，首先在變形物體周圍布置控制體，如圖1(a)所示，圖中，x/c為弦向位置,y/c為縱向位置，x、y分別為實際弦向及縱向坐標，c為翼型弦長。然后定義變形物體上各點實際坐標關于局部坐標的函數為

圖1 FFD方法及DFFD方法控制翼型變形Fig.1 Airfoil deformations controlled by FFD and DFFD methods

(1)

在定義的局部坐標系中，若FFD控制點的位移量為ΔPi,j,k，相應可得變形物體上各點的位移量為

(2)

則變形物體上各點的最終位置為

x′(s,t,u)=x(s,t,u)+Δx(s,t,u)

(3)

FFD方法雖然變形簡單，但變形物面的點和FFD控制點并沒有直觀的聯系，無法通過FFD控制點的變形精確控制物面某點的變形。DFFD方法在FFD方法的基礎上發展而來，其不同于一般FFD方法的是除在翼型周圍布置控制體和控制點外，在翼型上也布置直接操作點(pilot points)。通過直接操作點的位移反求FFD控制體上控制點的位移，進而通過FFD方法控制翼型的變化。由于將直接操作點的位移作為控制參數，因此可以使用高階的FFD控制體而不增加變量維數。反求的過程可看作求解滿足直接操作點位移約束的FFD控制點的位移組合，假設直接操作點的位移相互獨立，Sf(f=1,2,…,h) 為變形前直接操作點的原始坐標，Tf(f=1,2,…,h)為變形后各直接操作點的坐標。由式(2)可得

(4)

式中：δi,j,k為待求的FFD控制點的位移。

1.2 非線性區間分析方法

區間分析方法[16]采用區間數學的概念，將每一維不確定量v表示為區間形式：

AI=[AL,AR]={v|AL≤v≤AR,v∈R}

(5)

式中：上標I、L、R分別表示區間、區間下界、區間上界。由區間2個端點組成的一對有序實數稱為區間數。區間數的中值記作AC=(AL+AR)/2，區間半徑或離差記作Aw=(AR-AL)/2。

普通的線性區間函數問題可以通過區間數運算以及區間擴張獲得關于自變量的區間函數值，但跨聲速條件下的氣動問題屬于復雜的非線性問題，且沒有顯式的函數關系。本文通過優化方法進行非線性區間分析，即采用單目標粒子群進化算法，分別進行最大化目標、最小化目標共兩次尋優獲得目標函數區間的上界和下界。

除采用直接優化進行不確定分析外，本文同時采用構建代理模型結合優化方法進行不確定分析，以提高效率。本文選擇具有優異的非線性函數近似能力的Kriging模型[17]作為代理模型，建立不確定變量與氣動性能參數的函數關系。即通過拉丁超立方設計在不確定性變量取值空間選取一定量樣本并進行CFD計算，通過樣本集建立Kriging模型，對模型分別進行最大最小化目標尋優。為進一步提高代理模型精度，在初始樣本的基礎上采用最大均方差準則(MSE)進行加點，即選取當前模型在建模空間內預測均方差最大的點進行精確CFD分析并加入到樣本集內，重新構建代理模型。

1.3 不確定性分析算例

本文針對NACA0012翼型，進行了考慮幾何不確定性的阻力性能分析。為簡化問題，考慮不確定性時仍保證翼型對稱。首先采用DFFD方法對上半翼型進行參數化，建立10×8階FFD控制體，并在翼型上選取7個直接操作點，其弦向位置如表1所示。為保證前后緣固定，首尾2個直接操作點固定，改變其余5個直接操作點的y向位移以控制翼型變化，因此共有5個不確定性變量，其變化范圍均為[-0.001 5,0.001 5]。采用粒子群算法優化求解時，除粒子的位置和速度約束在變化范圍內，無其他約束。

翼型的阻力預測采用課題組自主研發的基于Navier-Stokes方程的并行CFD求解器[18]，湍流模型采用k-ω兩方程模型，空間離散采用二階格式，時間推進采用LU-SGS方法，網格采用美國國家航空航天局(NASA)提供的449×129標準網格。圖2為馬赫數0.7、迎角4°、雷諾數4.06×106條件下，計算壓力系數Cp分布與試驗數據[6]的比較，計算精度可滿足要求。

不確定分析算例采用工況為迎角0°，馬赫數0.85，該條件下計算得阻力系數為0.0538。首先采用20個初始樣本建立了阻力系數與不確定變量的Kriging模型(Kriging模型1)。在此基礎上又加點7次，重新構建模型(Kriging模型2)。在不確定性變化區間內隨機選取80個樣本使用模型進行預測并采用精確CFD分析進行驗證，以檢驗代理模型精度，樣本阻力系數的相對誤差如圖3所示，采用20個初始樣本建立代理模型，相對誤差最大為5.47%，加點7次后，最大相對誤差減小為2.24%，該精度可以滿足工程應用。

表1 直接操作點x方向位置Table 1 Position of pilot points in x direction

圖2 NACA0012翼型計算與試驗壓力分布比較Fig.2 Comparison of pressure distribution between computation and experiment for NACA0012 airfoil

為驗證不確定性分析的準確性，采用蒙特卡羅方法進行不確定性分析，即在設計空間內選取200個隨機樣本進行精確CFD分析并統計，并與上述3種結果對比，直接優化、Kriging模型優化與蒙特卡羅方法結果比較如圖4所示。結果表明直接優化方法具備較高的精度，只通過20個初始樣本構建Kriging模型，得到的結果與直接優化結果存在一定誤差，通過加點重建Kriging模型，誤差進一步減小。直接優化中每一步所有CFD分析可并行進行，建立代理模型的20個初始樣本也可并行進行CFD分析，但加點過程無法并行進行，故3種方法的分析結果、相對誤差、CFD計算次數、并行計算時間比較如表2所示。其中，采用Kringing模型1最大相對誤差小于5%，而分析效率可提高95%。3種方法得到的阻力系數變化區間上下界對應的翼型及壓力分布如圖5所示，由圖中可以看出，各方法得到的翼型相似。

圖3 Kriging模型阻力系數相對誤差Fig.3 Drag relative coefficient errors of Kriging models

圖4 直接優化、Kriging模型優化與蒙特卡羅方法結果比較Fig.4 Comparison of results among direct optimization, Kriging model optimization and Monte Carle method

表2 3種方法的分析結果、誤差、CFD計算次數、計算時間比較Table 2 Comparison of analysis results, errors, CFD calculation times and computing time among three methods

圖5 阻力系數變化區間上下界對應的翼型及壓力分布Fig.5 Airfoils and pressure distribution corresponding to upper and lower bounds of drag coefficient variation interval

2 魯棒優化設計

2.1 優化算例

本文采用RAE2822翼型優化設計標準算例，流場求解采用前述的CFD求解器。優化問題描述為

(6)

式中：CD為阻力系數目標；X為描述翼型幾何的n維設計變量；V為設計空間；考慮升力系數CL、俯仰力矩系數CM以及面積Aa約束，Ainitial為RAE2822翼型的面積，大小為0.077 873。計算工況為馬赫數0.734，雷諾數6.5×106。圖6為迎角2.79°條件下該翼型的計算與試驗壓力分布比較結果。

首先對該問題進行傳統的確定性優化。通過試驗，采用DFFD方法對翼型進行參數化時，若參數的變化范圍較大，極易使翼型產生大的波動，從而不具備翼型的基本幾何特征，類似問題也出現在其他直接以物面點位移參數化的方法中[19]。因此，DFFD方法雖可以直接控制翼型表面變形，但并不適合較大空間內的翼型參數化建模，故本文采用CST方法對翼型進行參數化，作為確定性的設計變量，翼型上下部分分別采用5階CST參數化。由于進行CFD計算時直接改變迎角進行定升力計算，因此假設升力系數等式約束直接成立，實現時，給定閾值[0.823,0.825]，檢測升力系數穩定收斂到閾值內結束計算。不等式約束處理采用罰函數法，約束函數值計算式為

(7)

優化算法采用單目標粒子群算法，優化結果以及RAE2822翼型計算結果如表3所示，α為迎角。對本文的優化結果與其他文獻[19-21]進行了對比，數值結果如表4所示，圖7為優化后的翼型對比，圖8為優化后壓力分布對比，本文優化后阻力系數下降35.9%，且滿足各項約束。與其他文獻結果對比，在基本相同的設計變量個數下，本文取得了較好的結果。

圖6 RAE2822翼型計算與試驗壓力分布比較Fig.6 Comparison of pressure distribution between computation and experiment for RAE2822 airfoil

表3 RAE2822翼型及優化翼型的計算結果Table 3 Computing results of RAE2822 and optimized airfoils

表4 優化翼型與其他文獻結果對比

圖7 優化后翼型對比Fig.7 Comparsion of optimized airfoils

圖8 優化后翼型壓力分布對比Fig.8 Pressure distribution comparison of optimized airfoils

2.2 魯棒優化設計方法

在2.1節算例的基礎上考慮翼型幾何的區間不確定性，問題(6)可描述為

(8)

式中：U為q維不確定向量，用q維區間向量UI表示。目標函數和約束函數均為不確定變量U的連續函數，故其可能取值也構成區間形式，而不是確定值，因此，上述問題無法通過傳統的確定性優化方法求解。

魯棒優化設計中，將不確定性變量導致的目標性能區間作為設計目標，因此在優化過程中，需比較不同樣本的目標性能區間的大小。本文中，利用≤Cw區間序關系[16]比較2個區間AI和BI的大小，對于最小化問題其定義如下：

(9)

P(AI≤BI)=

(10)

(11)

如2.1節所述，升力系數等式約束自動滿足，其他約束采用罰函數法處理，則問題(11)進一步轉換為無約束優化問題：

(12)

式中：

本文基于區間不確定性分析方法及自適應多目標粒子群算法，建立了考慮幾何區間不確定性的魯棒優化設計方法，求解上述優化問題，其流程如圖9所示。

2.3 自適應多目標粒子群算法

本文采用粒子群算法具有形式簡潔、收斂快速以及參數調節機制靈活等優點。在粒子群算法中，一個粒子i可由位置向量xi=[xi,1,xi,2,…,xi,d]T∈Rd和速度向量vi=[vi,1,vi,2,…,vi,d]T∈Rd表示，其中i=1,2,…,N，N為種群中的粒子個數，d為設計變量的維數。粒子在進化過程中根據式(13)更新速度和位置：

(13)

式中：t為迭代次數；ω≥0為慣性權重；c1、c2為加速系數；r1、r2為在[0,1]上均勻分布的隨機數；pbesti為第i個粒子的歷史個體最優解;gbest為群體的全局最優解。為避免陷入局部最優，可采用精英學習策略進行局部極值擾動，精英學習率為Lr。本文自適應多目標粒子群算法基于Pareto占優原則建立并維護外部最優解集，通過引入Pareto熵表示Pareto前緣的分布均勻性，間接表示種群的多樣性。優化迭代時，Pareto前緣發生變化，相應的Pareto熵也發生變化，故可以采用差熵表示Pareto前緣重新分布范圍的大小，通過差熵可以推測當前種群發現新解的能力，估計種群所處的進化狀態包括收斂狀態、多樣化狀態及停滯狀態，進而調整參數ω、c1、c2、Lr以設計適應當前狀態的勘探和開采策略。

圖9 考慮幾何不確定性的魯棒優化設計流程Fig.9 Robust optimization design process considering geometric uncertainties

優化算法的具體流程如圖10所示，文獻[22]采用此優化算法對不同測試算例進行了驗證，并與其他常用多目標優化算法進行了對比，結果顯示該算法總體上具有更好的多樣性和收斂性。

2.4 結果分析

采用2.2節魯棒優化設計方法對問題(12)優化。幾何不確定性變量仍采用DFFD方法參數化產生，對超臨界翼型的上半部分和下半部分分別進行參數化，各布置7個直接操作點，采用16×12階控制體，直接操作點及控制體建立如圖11所示。翼型前緣和后緣的直接操作點仍保持固定，故共10個不確定性變量，變化區間仍為[-0.0015 ,0.001 5]。為提高計算效率，對計算精度與計算時間進行折中，選取通過25個初始樣本構建Kriging模型的方式進行阻力系數及力矩系數的不確定性分析，直接優化的方式進行翼型面積的不確定分析。代理模型精度校驗得阻力系數預測的最大誤差為7.4%，力矩系數預測的最大誤差為0.63%。多目標優化所得Pareto解集如圖12所示。

由于多目標優化中粒子群容量及進化步數設置較小，且跨聲速問題機理復雜，所得的Pareto前緣不甚明顯，但工程中只需選取一個或幾個最優解即可，故選取圖中Pareto前緣的一點作為本次魯棒優化的最優解，并對最優翼型進行CFD分析，得其阻力系數CD；同時采用相同的參數化及代理模型建立方法對確定性優化得到的最優翼型及初始RAE2822翼型進行一次不確定性分析，其相關數值如表5所示。2種優化方法得到的最優翼型及對應壓力分布如圖13所示。對RAE2822翼型及最優翼型進行不確定性分析，得到阻力目標、力矩與面積約束區間上下界對應的翼型及壓力分布，如圖14所示。

圖11 直接操作點及FFD控制體Fig.11 Pilot points and FFD control body

圖12 多目標魯棒優化的Pareto解集Fig.12 Pareto set of multi-objective robust optimization

首先從優化翼型的CFD計算結果可看出，無論是確定性優化還是魯棒優化最優翼型的阻力系數較初始RAE2822翼型均有較大程度的下降。從翼型幾何和對應壓力分布情況來看，優化使得翼型上表面前緣下降后緣抬升，造成前緣吸力峰值更高，但減弱了激波強度，從而減小阻力。確定性優化結果較魯棒優化結果激波更平緩，阻力更小。從考慮阻力及其他約束的綜合不確定性分析的結果來看，原始RAE2822翼型以及確定性優化的區間半徑目標fw較魯棒優化結果差距很大，但從表5也可看出,區間中值目標fC與CFD計算結果CD相差較大，說明不確定分析結果fw和fC較大主要是由約束函數造成的。

圖15為不確定分析時各采樣樣本的目標函數和約束函數值比較，該結果一定程度上可顯示優化目標和各約束對不確定分析結果的影響。其中阻力系數結果顯示，確定性優化結果的整體阻力系數較低且變化區間較小，但從力矩和面積約束的結果可看出，確定性優化的結果均在一定程度上違反了約束，而魯棒優化結果則基本滿足約束。這是因為確定性優化時，由于目標和約束的沖突，所得最優結果一般分布在約束邊緣，當考慮不確定性因素時，極易違反約束。由此可見，本文魯棒優化流程考慮了不確定性對優化目標及約束的影響，相對確定性優化，得到了綜合魯棒性更強的結果。實際操作中，也可以根據決策者對各優化目標和約束的偏好，靈活設定約束罰函數因子，得到更理想的結果。

表5 優化結果比較Table 5 Optimization result comparison

圖13 優化所得最優翼型及壓力分布比較Fig.13 Optimized airfoils and pressure distribution comparison

圖14 目標及約束區間上下界對應的翼型及壓力分布Fig.14 Airfoils and pressure distribution corresponding to upper and lower bounds of objective and constraint intervals

圖15 不確定分析采樣樣本的阻力系數、力矩系數及面積對比Fig.15 Drag coefficient, moment coefficient and area comparison of samples for uncertainty analysis

3 結 論

1) 對NACA0012翼型的阻力系數性能進行了跨聲速條件下考慮幾何不確定性的快速非線性區間分析。通過建模方法與直接分析方法對比，最大誤差在5%以內，而分析效率可提高95%。

2) 以跨聲速條件下RAE2822翼型的阻力性能為目標，并考慮升力、力矩及面積約束，分別進行了翼型的確定性優化設計及魯棒氣動優化設計。確定性優化所得最優翼型較初始翼型阻力減小35.9%，且嚴格滿足各項約束，同等精度下優于其他文獻結果。魯棒優化所得最優翼型較初始翼型和確定優化最優翼型，其魯棒性得到了大幅提高。分析總體優化目標各項可知，本次優化結果主要通過降低違反約束的概率提高魯棒性，通過調節阻力目標與其他約束的罰函數因子，可得到偏好不同的結果。