基于生成拓撲映射的氣動外形反設計方法研究

宋超, 劉紅陽, 周鑄, 羅驍, 李偉斌

(中國空氣動力研究與發展中心 計算空氣動力研究所, 四川 綿陽 621000)

飛行器氣動外形直接決定其飛行性能與飛行品質。氣動外形與性能之間存在復雜關系，并受多種約束條件限制，使得氣動外形設計成為整個飛行器設計流程中最為困難的環節。氣動外形數字化設計技術充分利用計算流體力學、數值優化方法，在提高設計效率的同時，有助于獲得更優異的氣動性能。氣動外形數字化設計技術包括優化設計方法和反設計方法。優化設計方法魯棒性好，對設計經驗依賴程度低，但依賴于數值優化算法的特性，如全局優化算法需要大量的流場分析，計算效率低，而基于梯度的優化設計結果強烈依賴于初始構型在設計空間中的位置。

反設計方法通過改變氣動外形，使得其表面壓力分布與給定的目標壓力分布盡量吻合，該方法通常計算效率較高，在給定優秀目標壓力分布后能夠很快得到滿意的氣動外形。Lighthill[1]最早提出了基于保角變換的翼型反設計方法。之后，Takanashi[2]發展了余量修正法及其改進方法[3-4]，該方法是工程常用的反設計方法之一。此類經典的反設計方法要求給出合理的目標壓力分布，否則無法收斂。Bui-Thanh等[5]利用本征正交分解(proper orthogonal decomposition,POD)降階方法發展了Gappy POD翼型反設計方法。白俊強等[6]對Gappy POD進行了改進，引入了迭代校正法提高反設計精度。Gappy POD反設計方法對目標壓力分布合理性沒有要求，計算效率高，但是不能處理各種幾何和氣動約束。

近年來，機器學習技術不斷發展并被應用到翼型反設計領域。單志輝[7]基于高斯過程回歸模型進行了翼型反設計研究，壓力分布和外形數據用于模型訓練，以預測目標壓力分布對應的翼型。Sun等[8]基于人工神經網絡開展了翼型反設計研究。Sekar等[9]利用深度卷積神經網絡將翼型和壓力分布的圖像作為學習對象，實現了壓力分布圖像特征提取，從而計算獲得翼型外形。目前，基于神經網絡，尤其是深度神經網格的反設計方法還存在模型訓練復雜、訓練時間長等問題。

生成拓撲影射(generative topographic mapping,GTM)[10]是機器學習領域的一種非線性隱變量模型，能實現高維數據的非線性降維，廣泛應用于數據分析及可視化分析[11-13]。GTM能夠盡可能地保持數據在高維空間上原有的拓撲相對關系，并且建立了高維與低維空間的雙向映射。Viswanath等[14-15]利用氣動外形與氣動力組成的高維數據作為GTM模型訓練數據，進行了翼型的氣動力優化設計，是GTM模型在氣動設計領域的首次應用。

本文針對反設計方法目標壓力分布難以給定、對設計經驗依賴度高的問題，利用GTM模型的特點，建立了翼型設計變量與壓力分布的映射模型，并結合全局優化算法實現了翼型壓力分布高效反設計。

1 基于生成拓撲映射的氣動反設計方法

1.1 生成拓撲映射模型

GTM能夠實現高維數據到低維空間的映射，并保持原有的拓撲相對關系。設有D維數據集T={t1,t2,…,tN}，但數據本質是L維的(L

圖1 GTM原理示意圖

文獻[10]給出了GTM模型詳細的描述，采用中心位于y(x,W)且方差為β的高斯函數作為D維數據集T的分布

對x積分可得：

(2)

式中，p(x)為先驗分布

(3)

表示均勻分布在隱變量空間內的K個點上。

由此得到

(4)

通過最大對數似然函數可求得參數W和β

(5)

分布p(t|W,β)為中心固定的帶約束混合高斯模型。矩陣W是隱空間點向高維空間映射函數的權重參數矩陣，其維度為D×M,M為基函數個數。GTM算法可歸納為：已知隱變量空間中x的分布和數據空間的數據集合T，使用最大似然法求解參數W,β。常采用期望最大化(expectation maximization,EM)算法訓練GTM模型參數[16-17]。

1.2 基于GTM的氣動反設計方法

基于GTM的氣動反設計方法流程如圖 2所示。

圖2 基于GTM的目標壓力分布反設計流程圖

首先在根據氣動外形參數化方法確定的設計變量確定設計空間。然后在設計空間內對設計變量進行抽樣，得到樣本集S，并計算對應外形的流場，這一步驟為實驗設計。若設計變量個數為m，樣本點數目為n，則S為n×m的矩陣。利用流場求解器得到樣本點i對應的壓力分布，若繞氣動外形的壓力分布點數為k，則得到維度為n×k的壓力分布數據集Cp。將樣本集與壓力分布數據集合并，作為GTM模型的訓練高維數據Xtrain。Xtrain是維度為n×(m+k)的矩陣，每一行由一個樣本點及其對應的壓力分布數據組成。

GTM模型訓練完成后，算法建立了高維空間與隱空間的映射關系。隱空間的維度即為L維，給定隱空間的任意位置x，其中x是一個長度為L的向量，每個分量在±1范圍內。利用GTM反映射關系可以得到對應高維空間的數據，即設計變量與壓力分布組成的數據。因此，對于隱空間中的每一個點，可以通過GTM映射得到一個氣動外形及其對應的壓力分布。

定義表面壓力分布與目標壓力分布之間的差異為Ep

(6)

本文采用遺傳算法在隱空間進行尋優，目標函數為Ep。GTM模型充當類似代理模型的作用。遺傳算法的設計變量為隱空間變量，通過調用GTM模型，得到對應的壓力分布和氣動外形參數化變量。通過公式(6)可以得到隱空間變量對應的適應度值，供遺傳算法執行交叉變異等操作。通過遺傳算法得到最接近目標壓力分布的氣動外形參數化變量，進一步由參數化方法得到對應的氣動外形。

隱空間維度較低，因此在種群數目及迭代步數較少的情況下即可高效找到全局最優解。最后通過CFD校驗設計結果。若不滿足設計目標，可以進行設計樣本的更新及GTM模型的更新，或者調整目標壓力分布，進行新一輪的設計。

2 模型參數影響分析

GTM模型參數對模型精度有重要影響，進而影響氣動外形反設計結果，有必要對模型參數的影響規律進行研究。以NACA2412翼型為基本翼型，采用類函數/型函數參數化方法(class function/shape function transformation,CST)[19]，翼型上下表面型函數階數均取6，共14個設計變量，設計變量上下限取為基本翼型參數的±30%。計算狀態為：α=2.2,Ma=0.2,Re=4.5×106。采用拉丁超立方抽樣方法(Latin hypercube sampling,LHS)取50個樣本點，并利用Xfoil求解器得到翼型壓力分布。GTM模型隱空間維度、隱空間網格點數、基函數個數取值如表1所示。

表1 GTM參數與預測誤差的關系

在隱空間均勻取2 500個樣本點，利用GTM映射回高維數據空間，得到翼型的設計變量與模型預測的壓力分布。利用Xfoil求解器計算得到翼型的真實壓力分布，并根據(5)式定義計算預測壓力分布的誤差，結果在表1中列出。從表中可以看出，模型預測的誤差最小值、均值、標準差等總體變化趨勢為隨隱空間維度的增加而減小，隨隱空間網格點數與基函數個數的增加而減小。采用4維隱空間，204個隱空間網格點、10個基函數訓練GTM模型，選取預測誤差最大的一組壓力分布，如圖4所示。測試表明GTM模型能夠建立高精度的映射關系，滿足壓力分布反設計的要求。綜合考慮GTM精度、計算機內存、訓練時間等因素，本文選取的隱空間維度為4，網格點數為204，基函數個數為10。

圖3 GTM模型預測的壓力分布與Xfoil驗證結果比較

3 翼型反設計算例

3.1 低雷諾數翼型壓力反設計

本節根據1.2節的設計框架開展翼型壓力反設計研究。基準翼型為NACA2412，設計狀態為α=2.2,Ma=0.2,Re=4.5×106，壓力分布通過Xfoil求解得到。為驗證方法的穩健性，給出不同的目標壓力分布形態。遺傳算法的種群數目取為20，迭代步數取為30步。在算法得到最優解之后，利用求解器進行驗證。

圖4給出了基本翼型形狀及其壓力分布，以及輸入的目標壓力分布、求解器驗證得到的設計壓力分布，其中目標壓力分布共給出了4種不同的形式。本文方法不需要給出繞翼型一周的完整壓力分布，也不要求壓力分布具有實際物理意義。給出的目標壓力分布約在弦長40%～50%位置處。

圖4 給定不同目標壓力分布的反設計結果

圖4a～4b)給出的壓力分布較為合理，設計壓力分布與目標壓力分布幾乎重合，達到了設計預期。圖 4c)～4d)給出的壓力分布形態不合理，尤其是圖4d)“W”形的壓力分布形式顯然不具有物理意義，設計方法仍然能夠找到與目標壓力分布最為接近的結果，證明了設計方法的穩健性?？傊?，從翼型壓力分布反設計的結果來看，設計得到的壓力分布都很好地趨近于目標壓力分布。

3.2 跨聲速翼型壓力反設計

本節以NACA0010翼型為初始翼型，進行壓力分布反設計研究。設計狀態為α=0,Ma=0.76,Re=7.7×106。翼型壓力分布利用PMB3D求解器得到[20]，湍流模型為Spalart-Allmaras(SA)模型。計算網格如圖 5所示，網格單元數為21 285，物面第一層網格數滿足y+=1，遠場距離翼型30倍弦長。

圖5 計算網格

翼型采用CST參數化方法，設計變量數目為14。采用拉丁超方方法選取50個樣本點，并計算對應的壓力分布。給定的目標壓力分布參考了超臨界翼型的壓力分布，翼型上表面壓力分布平緩，翼型下表面前半部分壓力分布保持一定前加載，其值約束在-0.3左右，靠近后緣部分保持翼型后加載。第1輪反設計結果如圖6所示，設計壓力分布已經靠近目標壓力分布，但仍有一定差距。將第1輪的設計變量作為初始值，進行第2輪設計。第2輪的設計結果如圖7所示，設計壓力分布在翼型上表面、下表面前部基本與目標值重合，整體更加接近目標值。為更好逼近目標壓力分布，嘗試進行第3輪反設計，如圖8所示，設計壓力分布整體上更好地趨近于目標壓力分布，達到了理想設計結果。表2給出了3輪設計翼型的氣動特性，從表中看出，隨著設計迭代，翼型升力系數不斷提升，阻力系數不斷減小，翼型升阻力特性提升明顯。

圖6 跨聲速翼型壓力反 圖7 跨聲速翼型壓力反圖8 跨聲速翼型壓力反設計 設計結果(第1輪) 設計結果(第2輪) 結果(第3輪)

表2 3輪設計翼型的氣動特性

由于目標翼型與基本翼型外形差異較大，本次設計進行了3輪設計，不斷更新設計空間，取得了理想的設計效果。反設計問題中，無法確定事先給定的設計空間是否包含目標氣動外形。本次設計也很好反映了實際設計過程的這一特點。即使設計空間未能包含目標壓力分布對應的氣動外形，本方法依然表現出很好的穩健性，通過不斷調整設計空間，能夠達到理想外形與壓力分布。

3.3 三維層流短艙反設計

為了進一步驗證本文方法適應性，選取客機發動機短艙模型為研究對象，以提高短艙外表面層流區域為目標進行壓力分布反設計。選取的短艙模型長度5.54 m，最大直徑3.719 m，模型左右對稱，上下非對稱。計算網格采用結構網格，網格量為508萬。設計狀態為馬赫數Ma=0.85，雷諾數Re=3.1×106，自由來流湍流度Tu=0.1%，攻角α=2°，側滑角β=0°，湍流采用SST模型求解，轉捩模型采用γ-Reθ模型。采用FFD參數化方法單獨對短艙外形進行參數化。選取位于最左側的截面為設計剖面，如圖9所示，將圖中所示紅色圓點的沿短艙徑向坐標變化量作為設計變量，共9個設計變量。

圖9 短艙模型及剖面參數化方法

采用拉丁超方方法選取50個樣本點，獲取短艙流場信息，并提取短艙設計剖面外形和對應的壓力分布作為訓練樣本。為提高短艙外表面層流區域范圍，提高10%～30%弦長范圍內的壓力梯度，給定的目標壓力分布如圖10所示。

圖10 短艙剖面壓力分布反設計結果

圖中的設計壓力分布是經CFD驗證的，結果表明反設計的壓力分布很好地與目標壓力分布吻合，達到了設計意圖。設計后的剖面前緣半徑增大，后段相對厚度有所減小。圖11給出了基本外形和反設計得到外形外表面對應的摩阻分布。從圖中可知，層流區域延長了約弦長長度的4.4%。

圖11 短艙剖面摩阻分布反設計結果

圖12給出了基準短艙與設計短艙的摩阻云圖，從圖中可以明顯看出經過設計的剖面處層流區域有明顯延長。

圖12 基準短艙與設計短艙摩阻云圖

4 結 論

本文結合GTM機器學習方法與全局優化算法，發展了一種簡單、魯棒、高效的氣動外形目標壓力反設計方法。

1) GTM模型能夠建立氣動外形及其壓力分布數據在低維隱空間的準確映射關系，借助GTM模型，能夠準確預測壓力分布對應的氣動外形；

2) GTM模型將高維壓力分布數據降低到低維隱空間，有利于全局優化算法高效尋優。另外設計目標定義靈活，不要求壓力分布物理上存在，不要求給出繞氣動外形一周的全部壓力分布，極大提升了反設計方法的魯棒性，使得設計問題定義更加靈活，使用更加方便。

3) 本文方法需要樣本點較少，計算效率高，并且不局限于二維設計，在三維外形的設計中同樣具有很好的適應性。在下一步研究中，本方法將用于復雜三維氣動外形設計，進一步驗證方法的實用性。