夏陳超，邵 純，姜婷婷，陳偉芳

(浙江大學 航空航天學院,杭州 310027)

?

基于FFD方法的高超聲速升力體氣動優(yōu)化

夏陳超，邵 純，姜婷婷，陳偉芳

(浙江大學 航空航天學院,杭州 310027)

將FFD(Free Form Deformation)自由變形法與無限插值動網格方法相結合，發(fā)展了一種飛行器參數(shù)化建模和網格生成方法。二維和三維的實例顯示自由變形之后得到的飛行器幾何外形及其對應的網格能保持平滑光順，驗證了方法的有效性。在此基礎上，結合徑向基函數(shù)代理模型和CFD技術發(fā)展了一套優(yōu)化設計方法并對高超聲速升力體外形進行了氣動優(yōu)化?；谧赃m應模擬退火算法的單目標優(yōu)化表明，在保持原有外形體積不減小的情況下，升阻比提高了1.28%；基于NSGA-II的多目標優(yōu)化得到了飛行器升阻比和體積的最優(yōu)解集，典型優(yōu)化外形的升阻比和體積分別提高了2.93%和2.49%。升力體的優(yōu)化結果表明了FFD方法的有效性和優(yōu)化設計方法的實用性。

自由變形；無限插值法；高超聲速；升力體；優(yōu)化設計

0 引言

隨著計算流體力學(CFD, Computational Fluid Dynamics)技術的發(fā)展和計算機水平的提高，基于CFD的高保真度優(yōu)化在飛行器設計中日漸增多。在優(yōu)化設計過程中，參數(shù)化建模和網格生成是兩個基礎而又十分重要的環(huán)節(jié)。高效的優(yōu)化設計系統(tǒng)應該包含快速、便捷和適用性強的參數(shù)化建模和網格生成模塊。常用的飛行器參數(shù)化建模方法有基于CAD法、解析法、多項式/樣條曲線法、自由變形法(FFD, Free Form Deformation)等[1]。相比于傳統(tǒng)的參數(shù)化方法，F(xiàn)FD方法可用較少的參數(shù)實現(xiàn)對復雜幾何外形的控制，操作便捷，變形能力強，同時能夠保持原有外形的光滑和連續(xù)，是較為理想的參數(shù)化方法。Taira等[2]基于FFD技術，建立了可視化的飛行器布局設計平臺；Lyu等[3]采用FFD方法對翼身融合體進行了極多設計變量的氣動外形優(yōu)化，降低了飛行器阻力。國內，王丹等[4]、黃江濤等[5]和王元元等[6]采用FFD技術分別對翼身組合體整流包、翼稍小翼和運輸機后體進行了優(yōu)化設計，取得了良好的效果。除了參數(shù)化建模，為提高高保真度優(yōu)化設計系統(tǒng)的效率，需要在幾何外形發(fā)生變化時，迅速形成新的網格用于數(shù)值計算。新的網格可以在優(yōu)化過程中重新生成[7]，更多的是采用動網格技術生成，原因是復雜外形飛行器的網格重新自動生成十分困難和耗時。常用的動網格方法有無限插值法[8-9]、彈簧法[10]和Delaunay圖映射法[11]等，其中無限插值法是一種簡單高效的結構動網格技術，在飛行器優(yōu)化設計、非定常數(shù)值計算中得到了廣泛應用。升力體是典型的高超聲速滑翔飛行器外形，具有較好的升阻特性，Kinney[12]和張珍銘等[13]采用工程估算方法分別對HL-20和類X-33升力體進行了優(yōu)化設計，馬洋等[14]采用CFD方法進行了滑翔式升力體外形優(yōu)化，均在一定程度上提高了飛行器升阻比。在采用CFD方法進行飛行器優(yōu)化時，計算量往往十分巨大，一種折中措施是采用代理模型的近似方法，Kriging模型[15]、徑向基函數(shù)模型(RBF, Radial Basis Function)[16]和響應面模型(RSM, Response Surface Model)[17]等是常用的代理模型，在飛行器設計中應用較多。

目前FFD方法主要應用于亞聲速流動，在高超聲速流動中的應用還很少。本文將FFD自由變形法與無限插值動網格方法相結合，建立了一種參數(shù)化建模和網格生成方法，接著基于徑向基函數(shù)代理模型、CFD方法和單目標/多目標優(yōu)化算法，對高超聲速升力體進行了氣動外形優(yōu)化設計。

1 FFD自由變形和動網格方法

FFD自由變形法由Sederberg和Parry[18]于1986年提出，其基本思想是將幾何模型嵌入到一個控制體中，建立數(shù)學關系后通過改變控制頂點來改變模型的形狀。本文采用Bernstein多項式來定義控制體，若在原點X0建立一個局部坐標系，則控制體中任意一點的全局坐標可寫為

X=X0+sS+tT+uU

(1)

其中，s、t、u為該點在局部坐標系中沿S、T、U方向上的值，可表示為

(2)

若局部坐標系3個軸矢量方向上分別有l(wèi)、m、n個控制頂點，則控制體將含(l+1)(m+1)(n+1)個頂點，其全局坐標可表示為

(i=0,1…,l；j=0,1,…,m；k=0,1,…,n)

(3)

對于待變形模型上任意一點，在求得其局部坐標(s,t,u)后，變形后新的全局坐標可由Bernstein多項式計算得到

(4)

實施FFD自由變形的主要步驟為：根據待變形幾何模型建立控制體和局部坐標系；計算幾何模型的局部坐標；改變控制頂點的全局坐標；根據式(4)計算幾何模型新的全局坐標。

在完成幾何變形之后，即可采用TFI方法進行計算網格的變形。TFI無限插值法[19]基于原有的多塊對接結構網格，通過歸一化的弧長參數(shù)將變形量傳遞至流場網格。以三維的體網格插值為例，其無限插值過程為

ΔSi,j,k=U+V+W-UV-UW-VW+UVW

(5)

其中

U=(1-αi,j,k)ΔS(1,j,k)+αi,j,kΔS(NI,j,k)

V=(1-βi,j,k)ΔS(i,1,k)+βi,j,kΔS(i,NJ,k)

W=(1-γi,j,k)ΔS(i,j,1)+γi,j,kΔS(i,j,NK)

UV= (1-αi,j,k)(1-βi,j,k)ΔS(1,1,k)+

αi,j,k(1-βi,j,k)ΔS(NI,1,k)+

(1-αi,j,k)βi,j,kΔS(1,NJ,k)+

αi,j,kβi,j,kΔS(NI,NJ,k)

UW= (1-αi,j,k)(1-γi,j,k)ΔS(1,j,1)+

αi,j,k(1-γi,j,k)ΔS(NI,j,1)+

(1-αi,j,k)βi,j,kΔS(1,j,NK)+

αi,j,kβi,j,kΔS(NI,j,NK)

VW= (1-βi,j,k)(1-γi,j,k)ΔS(i,1,1)+

βi,j,k(1-γi,j,k)ΔS(i,NJ,1)+

(1-βi,j,k)γi,j,kΔS(i,1,NK)+

βi,j,kγi,j,kΔS(i,NJ,NK)

UVW= (1-αi,j,k)(1-βi,j,k)(1-γi,j,k)ΔS(1,1,1)+

αi,j,k(1-βi,j,k)(1-γi,j,k)ΔS(NI,1,1)+

(1-αi,j,k)βi,j,k(1-γi,j,k)ΔS(1,NJ,1)+

(1-αi,j,k)(1-βi,j,k)γi,j,kΔS(1,1,NK)+

(1-αi,j,k)βi,j,kγi,j,kΔS(1,NJ,NK)+

αi,j,k(1-βi,j,k)γi,j,kΔS(NI,1,NK)+

αi,j,kβi,j,k(1-γi,j,k)ΔS(NI,NJ,1)+

αi,j,kβi,j,kγi,j,kΔS(NI,NJ,NK)

式中 ΔS為位移；α、β、γ為歸一化的弧長參數(shù)；NI、NJ、NK為網格塊3個方向編號的最大值。

2 自由變形和動網格實例

2.1 二維RAE2822翼型

RAE2822是經典的二維跨聲速翼型，常作為湍流模型的驗證算例。圖1為原始的翼型網格和FFD控制體，網格在近壁面進行了加密，圓點為控制頂點。圖2為控制頂點、翼型的變化以及對應的網格。可看出，翼型和網格在變形之后均可保持連續(xù)和光順。

圖1 原始網格和控制體Fig.1 Initial grid and control volume

圖2 變形后網格和控制體Fig.2 Grid and control volume after deformation

2.2 三維NACA64A010機翼

以三維NACA64A010機翼[20]進行驗證，圖3為機翼翼尖扭轉15°后的網格和拓撲，圖4為機翼向上翹曲的網格示意圖，其中圓點為控制頂點，控制體內為變化后外形?？梢钥闯觯現(xiàn)FD方法的變形能力較強，且變形后的機翼網格仍能保證較高的質量。

圖3 翼尖扭轉15°前后的網格Fig.3 Comparison of grid before and after 15°twist of wing tip

圖4 機翼翹曲前后的網格Fig.4 Comparison of grid before and after the warp of wing

2.3 類X-33升力體外形

以類X-33升力體[21]外形為例，檢驗自由變形和動網格方法對復雜飛行器外形的適用性。圖5為飛行器半模網格，網格量約140萬，共含29個網格塊。圖6為飛行器頭部迎風面局部變形示意圖，控制頂點的移動改變了飛行器表面局部形狀?？梢钥闯?，飛行器幾何外形在經過FFD變形后仍能保持很好的光滑性。FFD技術的這一特點使其十分適合于飛行器的局部修形和精細優(yōu)化。

圖5 類X-33網格Fig.5 Grid of the X-33-type vehicle

圖6 類X-33頭部迎風面FFD變形Fig.6 Free form deformation of the windward side of the X-33-type vehicle

3 升力體外形優(yōu)化設計

基于發(fā)展的自由變形和動網格技術，構建優(yōu)化系統(tǒng)并針對一類簡化的高超聲速升力體外形進行優(yōu)化設計。升力體初始外形長2 m，底部為長、短半軸之比為3∶1的橢圓形。初始網格如圖7所示，網格量約60萬，劃分為18塊，網格在壁面進行了加密保證y+值小于1。包圍升力體的FFD控制體如圖8所示，放大圖為控制底部截面形狀的控制頂點及其編號。設計變量為1～6號控制點(考慮了左右對稱性)，保持飛行器長寬不變，控制點僅在Y方向移動。

圖7 升力體網格示意圖Fig.7 Grid of the lifting body

圖8 控制體和設計變量示意圖Fig.8 Control volume and design variables

3.1 數(shù)值計算方法

飛行器氣動特性的獲取通過求解可壓縮N-S方程，對流項采用二階TVD格式進行離散，時間項采用隱式推進方法，湍流模型為k-εSST[22]兩方程模型。進口為來流條件，出口為外推條件，壁面為等溫條件。

3.2 試驗設計和代理模型

構建合適的代理模型是解決高精度優(yōu)化中計算量過大的常用措施，其過程包括試驗設計(DOE, Design of Experiment)、樣本點計算以及代理模型的建立和驗證。本文升力體的設計變量如表1所示?；诖耍捎美〕⒎?LHS, Latin Hypercube Sampling)試驗設計方法形成40個樣本點，結合適用于非線性全局優(yōu)化問題的徑向基函數(shù)方法構造代理模型。

表1 設計變量范圍Table 1 Design variables and ranges

圖9給出了飛行器體積和升阻比的分析值與預測值之間的對比。其中，體積的分析值通過飛行器體積分得到，升阻比的分析值通過CFD計算得到，預測值均通過代理模型獲得。通過交叉驗證方式計算得到代理模型對飛行器體積和升阻比的均方根誤差(RMSE, Root Mean Square Error)分別為2.6%和3.3%，決定系數(shù)(R-Squared)分別為98.5%和98%?？梢钥闯?，構建的代理模型具有較高的精度，滿足優(yōu)化設計的需求。

(a)體積 (b)升阻比

3.3 優(yōu)化流程和算法

本文基于代理模型方法，構建了優(yōu)化設計系統(tǒng)，基本流程如圖10所示。在建立滿足精度要求的代理模型后，采用自適應模擬退火算法(ASA, Adaptive Stimulated Annealing)[23]進行升阻比的單目標優(yōu)化，采用非支配排序遺傳算法(NSGA-II, Non-dominated Sorting Genetic Algorithm)[24]進行升阻比和體積的多目標優(yōu)化。

圖10 優(yōu)化設計流程Fig.10 Flowchart of optimization

3.4 優(yōu)化結果與分析

針對本文的升力體外形，在馬赫數(shù)10，攻角10°，高度30 km條件下進行優(yōu)化，初始外形升阻比2.186，體積0.281m3。首先進行單目標優(yōu)化，約束條件為飛行器體積不小于初始外形。圖11為采用自適應模擬退火算法優(yōu)化后升阻比的收斂過程，最優(yōu)外形對應的升阻比為2.214，相比于初始外形提高了1.28%。通過對最優(yōu)外形進行CFD計算得到的升阻比為2.203，與預測值相對誤差僅為0.5%。從圖12的底部截面形狀對比可看出，最優(yōu)外形在飛行器對稱面上的厚度略有增加，在對稱面與翼前緣之間的厚度有一定減小。

用NSGA-II優(yōu)化算法進行體積和升阻比的多目標優(yōu)化，種群數(shù)量為400，進化代數(shù)為100。由于體積和升阻比相互制約，因而多目標優(yōu)化得到的是關于體積和升阻比的諸多可行解及最優(yōu)解集(Pareto前沿)，如圖13。圖中三角形點為初始外形，空心原點為Pareto前沿上的4個典型結果。對于外形1結果，其體積為0.676 m3，但升阻比只有0.95，對于外形4結果則相反，其升阻比雖高達2.83，但體積只有0.103 m3。

圖13中外形3的升阻比和體積分別為2.250和0.288 m3，采用CFD數(shù)值計算得到的升阻比和體積分別為2.199、0.290 m3，誤差分別僅為2.32%和0.69%。相比于初始外形的升阻比和體積分別提高了2.93%和2.49%。圖14為外形3底部截面與初始外形的對比，與單目標優(yōu)化結果類似，此時飛行器迎風面的對稱面位置明顯增厚，翼前緣與對稱面之間減薄，即圖8中序號為2的設計變量向下移動量較大，序號為6的設計變量向上移動量較大。從圖15可看出外形3升阻比比初始外形大的原因：外形3翼前緣附近的馬赫數(shù)等值線(實線)比初始外形(虛線)更加貼近前緣，即其下表面的高壓氣流更少地溢流到上表面，有類似于乘波體的效果，因而升阻比較大。綜合單目標和多目標優(yōu)化，對本文10°攻角的優(yōu)化狀態(tài)，迎風面形狀對升阻比影響很大。實際飛行器設計中應根據需求，按照不同權重對升阻比和體積折中選擇。

圖11 升阻比收斂過程Fig.11 Convergence history of lift-drag ratio

圖12 優(yōu)化外形和初始外形的對比Fig.12 Comparison of initial and optimized shapes

圖13 多目標Pareto前沿Fig.13 Pareto front of multi-objectiveo ptimization

圖14 外形3和初始外形的對比Fig.14 Comparison of initial and optimized shapes of configuration 3

圖15 翼前緣附近馬赫數(shù)等值線對比Fig.15 Comparison of Mach contours near leading edge

4 結論

(1)基于FFD自由變形法和TFI無限插值技術，發(fā)展了參數(shù)化建模和網格生成方法，通過對翼型、機翼和復雜外形飛行器的建模與動網格生成驗證了方法的有效性。FFD方法操作便捷，變形能力強，能保持幾何外形的光順，TFI方法簡易高效，二者的結合可為優(yōu)化設計提供良好的基礎，提高優(yōu)化效率。

(2)基于徑向基函數(shù)代理模型和CFD技術，發(fā)展了高保真度優(yōu)化設計方法，針對高超聲速升力體外形進行了單目標/多目標優(yōu)化設計。在保持初始外形體積不變的條件下，基于ASA算法的單目標優(yōu)化將升阻比提高了1.28%。采用NSGA-II算法得到了升力體升阻比和體積的最優(yōu)解前沿，最優(yōu)解典型外形的升阻比和體積較初始外形分別提高了2.93%和2.49%，相比于初始外形其迎風面腹部厚度增大，腹部兩側厚度減小。

(3)本文發(fā)展的優(yōu)化設計方法基本滿足飛行器優(yōu)化設計的需求，具有較好的工程實用性。在此基礎上將進一步開展飛行器氣動布局優(yōu)化設計工作。

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

[2] Taira T,Jeong S,Obayashi S,et al.GUI-based geometry deformation tool for modification of aircraft configurations[R].AIAA 2013-0619.

[3] Lyu Z,Martins J.RANS-based aerodynamic shape optimization of a blended-wing-body aircraft[R].AIAA 2013-2586.

[4] 王丹,白俊強,黃江濤.FFD方法在氣動優(yōu)化設計中的應用[J].中國科學:物理學力學天文學,2014,44(3):267-277.

[5] 黃江濤,高正紅,白俊強,等.基于任意空間屬性FFD技術的融合式翼稍小翼穩(wěn)健型氣動優(yōu)化設計[J].航空學報,2013,34(1):37-45.

[6] 王元元,張彬乾,郭兆電,等.基于FFD技術的大型運輸機上翹后體氣動優(yōu)化設計[J].航空學報,2013,34(8):1806-1814.

[7] 夏陳超,趙文文,陳偉芳,等.類HTV-2升力體參數(shù)化建模與網格自動生成研究[C]//第十五屆全國激波與激波管學術交流會文集.杭州,2012:503-508.

[8] Byun C,Guruswamy G P.A parallel,multi-block,moving grid method for aeroelastic applications on full aircraft[R].AIAA 98-4782.

[9] Tsai H M,F.Wong A S,Cai J,et al.Unsteady flow calculations with a parallel multiblock moving mesh algorithm[J].AIAA Journal，2001,39(6):1021-1029.

[10] Zeng D,Ethier C R.A semi-torsional spring analogy model for updating unstructured meshes in 3D moving domains[J].Finite Elements in Analysis and Design,2005,41(11):1118-1139.

[11] Liu X,Qin N,Xia H.Fast dynamic grid deformation based on Delaunay graph mapping[J].Journal of Computational Physics,2006,211(2):405-423.

[12] Kinney D.Aerodynamic shape optimization of hypersonic vehicles[R].AIAA 2006-239.

[13] 張珍銘,丁運亮,劉毅.升力體外形設計的代理模型優(yōu)化方法[J].宇航學報,2011,32(7):1435-1444.

[14] 馬洋,楊濤,張青斌.高超聲速滑翔式升力體外形設計與優(yōu)化[J].國防科技大學學報,2014,36(2):34-40.

[15] 楊希祥,周張,彭科.基于Kriging方法的整流罩氣動外形設計優(yōu)化[J].固體火箭技術， 2014,37(2):167-171.

[16] Leary S J,Bhaskar A,Keane A J.Global approximation and optimization using adjoint computational fluid dynamics codes[J].AIAA Journal,2004,42(3):631-641.

[17] 成沉,鮑福廷,劉旸,等.基于響應面法的喉栓式噴管型面優(yōu)化設計[J].固體火箭技術,2014,37(1):47-51.

[18] Sederberg T W,Parry S R.Free-form deformation of solid geometric models[J].Computer Graphics,1986,20(4):151-160.

[19] Thompson J F,Soni B K,Weatherill N P.Handbook of grid generation[M].CRC Press,1998.

[20] Chaderjian N M,Guruswamy G P.Transonic Navier-Stokes computations for an oscillating wing using zonal grids[J].Journal of Aircraft,1992,29(3):326-335.

[21] Murphy K J,Nowak R J,Thompson R A,et al.X-33 hypersonic aerodynamic characteristics[J].Journal of Spacecraft and Rockets,2001,38(5):670-683.

[22] Menter F R,Kuntz M,Langtry R.Ten years of industrial experience with the SST turbulence model[C]//Proceedings of the 4th International Symposium on Turbulence,Heat and Mass Transfer.West Redding,2003:625-632.

[23] Ingber L.Adaptive simulated annealing (ASA):Lessons learned[J].Control and Cybernetics,1996,25(1):33-54.

[24] Deb K,Pratap A,Agarwal S,et al.A fast and elitist multiobjective genetic algorithm:NSGA-II[J].IEEE Transactions on Evolutionary Computation,2002,6(2):182-197.

(編輯：呂耀輝)

Aerodynamic optimization of hypersonic lifting body based on FFD method

XIA Chen-chao,SHAO Chun,JIANG Ting-ting,CHEN Wei-fang

(School of Aeronautics and Astronautics,Zhejiang University,Hangzhou 310027,China)

A general method of parametric modeling and grid generation was developed by combining the free form deformation(FFD) and transfinite interpolation techniques.Smoothness of geometry shape and satisfactory quality of grid were achieved from two-and three-dimensional cases.An optimization framework was then constructed by adopting a surrogate model and optimization algorithms.The aerodynamic shape optimization of a hypersonic lifting body was utilized as a demonstration case for the developed framework.The result of single objective optimization using adaptive simulated annealing algorithm shows that the lift-drag ratio has an increase of 1.28% with a constraint on the volume.Besides,multi-objective optimization is carried out by using NSGA-II algorithm,and the Pareto solutions about volume and lift-drag ratio are obtained.One of the typical solutions on the Pareto front shows that the lift-drag ratio has an increase of 2.93% and the volume has an increase of 2.49%.Optimization of the lifting body in this work shows the usability of the free form deformation technique and the good practicality of the optimization framework.

free form deformation;transfinite interpolation;hypersonic;lifting body;design optimization

2014-08-15；

：2014-11-17。

國家重點基礎研究發(fā)展計劃(2014CB340201)。

夏陳超(1989—)，男，博士生，研究方向為飛行器優(yōu)化設計。E-mail：aeroxia@zju.edu.cn

V441

A

1006-2793(2015)06-0751-06

10.7673/j.issn.1006-2793.2015.06.001