面向飛行器概念設計的全速域氣動分析工具

呂凡熹, 李正洲, 鄧經樞, 肖天航, 余雄慶

(南京航空航天大學 航空宇航學院, 江蘇 南京 210016)

面向飛行器概念設計的全速域氣動分析工具

呂凡熹, 李正洲, 鄧經樞, 肖天航, 余雄慶*

(南京航空航天大學 航空宇航學院, 江蘇 南京 210016)

為滿足飛行器概念設計中快速氣動計算的需求，研究和開發了一種全自動化全速域氣動分析集成工具。針對不同的速度域，該工具由三種計算方法組成：1) 亞、跨聲速情況，采用基于自適應直角網格的非線性全速勢方程有限體積解法；2) 超聲速情況，采用基于自適應直角網格的歐拉方程有限體積解法；3) 高超聲速情況，采用基于當地表面斜度法的面元法。不同速度域的多種構型飛行器的算例結果表明，該集成工具的計算精度滿足要求，且自動化程度高、收斂速度快，可應用于飛行器概念設計的快速氣動分析。

飛行器；空氣動力學；概念設計；全速勢方程；歐拉方程；當地表面斜度法；自適應直角網格

0 引 言

概念設計是飛行器設計的早期階段，需要從眾多概念方案中優選出最佳方案。飛行器概念設計中需進行大量氣動特性分析，傳統的做法是采用工程估算方法快速地分析概念方案的氣動特性。目前飛行器概念設計的一個重要趨勢是氣動分析模型采用數值計算方法[1-2]，其優勢在于：能提高氣動分析模型的預測精度；能應用于非常規布局的氣動特性分析。

飛行器概念設計階段的氣動數值分析方法應滿足如下要求：1) 外形建模的參數化；2) 計算網格生成的自動化；3) 流場求解速度快、時間短；4) 能適用于不同飛行速度和高度范圍。雖然基于N-S方程的數值方法具有計算精度高，通用性強的優點，但所需計算時間仍然較長，目前還難以滿足概念設計階段快速響應要求。

國外已發展了幾種面向飛行器概念設計的氣動數值分析工具，如APAS[3]和CBAERO[4]。但上述兩個工具無法對跨聲速氣動特性進行較精確的估算。從計算精度和效率上看，在概念設計階段，基于全速勢方程的數值方法具有計算速度快、存儲空間占用小的特點，是眾多方法中最適用于亞/跨聲速計算的工具。目前，典型的基于全速勢方程的氣動分析程序有FLO[5]、BLWF[6]、TranAir[7]等。其中，FLO只適用于機翼，BLWF只適用于翼身組合體。TranAir采用自適應直角坐標網格和有限元方法，具有較好通用性，但TranAir的使用依賴附加的前處理模塊AGPS[8]，一定程度上也限制其自動化的程度。

網格生成是CFD中一項重要的內容。結構、非結構和直角網格是常用的幾種網格類型，與前兩種網格類型相比，直角網格有著空間填充能力強和復雜幾何外形適應性好的優點，較容易實現網格的自動生成和自適應加密。從滿足概念設計階段氣動計算對網格生成的要求來說，直角網格是值得研究和發展的方法[9-10]。

本文旨在研究和開發一種全速域氣動分析的集成工具。亞、跨聲速情況，采用全速勢方程有限體積解法；超聲速情況，采用歐拉方程有限體積解法；高超聲速情況，采用基于當地表面斜度法的工程面元法。計算網格采用自適應直角網格。通過集成參數化幾何模型、自適應直角網格方法、全速勢方法、歐拉方法、工程面元法的程序，形成面向飛行器概念設計的快速氣動分析工具。

以下首先闡述氣動分析工具的總體框架；然后說明幾何模型和計算網格的快速生成方法；之后簡述全速勢方法、歐拉方法、工程面元法的算法；最后通過算例對各氣動分析程序進行驗證。

1 氣動分析工具的框架

本文氣動分析工具由外形參數化建模、計算網格生成、流場計算組成，其中流場計算包含全速勢方法、歐拉方法、工程面元法。氣動分析工具的架構如圖1所示。通過集成各模塊，整個過程可自動進行。流程簡述如下：

1) 采用軟件Open VSP[11]進行飛行器幾何建模，快速生成的幾何模型文件(STL格式)。

2) 基于幾何外形，采用自適應直角網格方法自動生成全速勢和歐拉方法需要的網格。工程面元法則直接讀取STL文件作為氣動網格。

3) 針對不同來流速度選擇計算方法；分析結束后，根據結果的精度，全速勢和歐拉方法可自動執行網格自適應，直到得到足夠精度的結果。

4) 分析和處理計算結果，輸出分析報告和指定格式的數據文件存入飛行器氣動數據庫。

2 飛行器外形參數化建模

Open VSP是一種面向飛行器概念設計的參數化幾何建模開源軟件，特點如下：1) 與傳統CAD軟件相比，需要的計算機資源很少；2) 飛行器的幾何特征是由設計人員所熟悉的參數來定義的；3) 可快速生成概念方案的幾何模型；4) 提供了多種組件來定義飛行器的幾何特征，可涵蓋各類飛行器的外形建模。

Open VSP創建的飛機幾何模型的輸出圖形格式包括IGS、STL、STEP等，并可計算和輸出飛機各部件的外露面積、內部體積等幾何特性。同時，應用其腳本功能可實現從EXCEL或XML文件讀取給定的幾何參數，直接生成模型，并輸出指定格式的幾何文件。

3 計算網格的自動生成

高超工程面元法計算直接基于STL三角片數據進行。全速勢和歐拉方程的求解則在自動生成的直角網格上進行，生成過程如下。

3.1幾何模型和數據結構

采用外形建模模塊生成的STL格式的文件作為幾何模型。該格式的文件由包含外法向量信息的三角片組成。通過使用交替數字樹數據結構[12]來管理三角片，使得直角網格單元搜索相交三角片的操作能快速進行。

采用八叉樹結構定義網格的拓撲關系。完整的直角網格是從包含整個計算域的根節點逐步加密得到的。除根節點外，八叉樹中每個單元都是通過加密其父單元得到的。

3.2網格初始化和幾何自適應

從八叉樹的根節點開始生成初始均勻網格。將樹中的每個葉子節點加密，直至所有葉子節點層數加密至設定的初始層數。

顯然，初始均勻網格不能精確地描述幾何外形，故需要進行幾何自適應[9-10]。幾何自適應是根據物面邊界的復雜程度，通過判斷物面邊界單元是否需要加密來實現的。

根據幾何自適應和后續處理的需要，需將網格單元進行分類。直角網格單元主要分為三類：流場內非切割單元、固壁內非切割單元和物面邊界單元。使用交替數字樹算法搜索確定與物面相交的單元，即物面邊界單元。當所有物面邊界單元都判斷完畢，流場內外非切割單元便組成了各自的單連通域。非切割單元使用光線投射算法[10]進行類型判斷，并使用染色算法快速確定其單連通域內其他單元的類型。

3.3物面邊界單元的處理

因為直角網格并不貼體，所以需要對物面邊界單元進行處理才能正確的施加物面邊界條件。本文使用精度高且發展成熟的切割法[9-10,13]進行處理。由于物面邊界單元和物面邊界的相交情況十分復雜，因此八叉樹結構無法方便地存儲切割后的單元信息，見圖2。例如圖2(a)所示，當單元被分割成多個部分時，八叉樹結構就很難描述該拓撲關系了。本文采用“cell-to-face”[14]的數據結構來專門存儲物面邊界單元，通過網格面存儲的兩側單元編號直接提供單元間連接信息。這樣不同種類的復雜切割情況能夠得到正確的處理，其拓撲關系也能得到正確的存儲，如圖2(b)。

切割處理后會產生體積極小的單元，從而引起計算的不穩定。故采用網格融合[13]技術將極小單元與較大的鄰居單元合并，如圖2(c)所示。

4 數值方法

4.1全速勢方程及離散求解方法

4.1.1 控制方程

三維定常全速勢方程的守恒形式為:

式中:速度V是速度位φ的梯度。ρ是表達式為：

其中，γ是比熱比，Ma∞是自由來流馬赫數。

4.1.2 有限體積空間離散

1) 離散格式。將控制方程寫成積分形式，由散度定理得到:

對任意單元i，方程(4)用有限體積法離散為:

其中，ρj、Vj是單元i的網格面j面心處的密度和速度，nj、ΔAj是網格面j的外法矢和面積，Resi表示單元i的殘差。

2) 通量計算。采用文獻[15]中的K-exact Least Squares方法構建φ的梯度來計算方程(5)中的Vj。在面心進行梯度構建時，選擇面左右單元以及左右單元鄰居的值進行插值。每個參與單元的插值多項式為:

式中：xi、yi、zi為單元i中心的坐標。由于參與單元的個數不確定，采用最小二乘法求這四個系數。插值多項式在每個單元上的誤差為:

將各單元誤差的平方加權求和可得:

式中,ri是各單元中心相對面心距離的倒數。對F求偏導，得到方程組，該方程組的解分別是速度分量u、v、w。

式中,Δx、Δy和Δz分別是迎風單元中心到面心連線的向量分量的兩倍，ρx、ρy和ρz是迎風單元中心處的密度梯度。m表達式為:

Ma是迎風單元的馬赫數。Mac是截斷馬赫數，常取0.95。C用來調整人工粘性大小，常取1到2之間。

4.1.3 GMRES隱式迭代求解

隱式格式每步求解方程組:

使用牛頓線性化的操作得到線性方程組:

本文使用結合ILU預調制方法[15]的GMRES算法解方程(12)。

4.2歐拉方程及離散求解方法

三維積分形式的歐拉方程為：

其中：W=[ρρuρvρwρe]T為守恒變量，其中ρ為流體密度，u、v和w為笛卡爾坐標系3個方向的絕對速度分量，e為單位質量總能；F(W)為對流通量。歐拉方程在空間上采用迎風格式有限體積離散，在時間方向上運用全隱式推進求解，離散方程采用GMRES格式迭代求解。

4.3流場自適應的計算

由于計算資源有限，為準確呈現流場特性，本文實施流場自適應。歐拉方法選擇速度散度和旋度作為考察參數[9]。考慮到單元大小的影響，針對速度散度的參數的表達式為:

τdi=|

對于旋度的考察參數，同樣有

τci=|

4.4粘性阻力的計算

全速勢和歐拉方法包含無黏假設只能計算誘導阻力。為計算黏性阻力系數，采用基于附面層理論和部件形狀因子修正方法計算各部件的黏性阻力系數[3]，通過疊加各部件的黏性阻力系數，獲得全機的黏性阻力系數。

4.5基于當地表面斜度法的面元法

對于高超聲速計算問題，可以采用歐拉方法，但對于特定的問題計算速度低且穩定性差。因此采用基于當地表面斜度理論的面元法計算氣動力,并采用三角形面元描述幾何外形。面元按照撞擊角分為迎風面和背風面，根據不同馬赫數選取不同公式計算面元的壓力系數。對于迎風面，Ma∞為1.5～3.8時，采用活塞-牛頓撞擊理論[16]；Ma∞為3.8～10時，采用修正牛頓撞擊理論[17]；Ma∞>10時，采用活塞-牛頓撞擊理論。方法如下：

式中：

對背風面，則采用普朗特-邁耶膨脹波方法：

5 算例驗證

5.1全速勢方法

5.1.1 亞聲速算例驗證

亞聲速算例是針對DLR-F4進行數值分析，自由來流Ma∞=0.6，CL=0.5。4.3節中采用的流場自適應方法的效果在文獻[18]已經得到了驗證，并應用到了本文算例中。如圖3所示，計算網格進行了流場自適應，并且加密了流動參數變化劇烈的區域。圖4給出執行5次流場自適應后的收斂歷程圖。對同套網格，相比于課題組開發的歐拉方法求解器，該全速勢方法計算時間縮短75%。本算例網格數為295 412個，不考慮流場自適應的情況下，在Intel Xeon E5-2630V2 處理器上單核運行計算程序(下文所有算例均采用相同計算條件)，僅需12個迭代步和240 s的物理時間。

圖5給出了壓力系數云圖。圖6給出了本文全速勢和歐拉方法于7個展向站位的壓力系數和實驗數據[19]的對比。可見前緣處全速勢方法的精度更高，兩種方法的結果都與實驗數據相吻合。

5.1.2 跨聲速算例驗證

跨聲速算例針對ONERA M6進行數值分析，自由來流Ma∞=0.839,α=3.06°。采用流場自適應來更精確地捕捉激波位置，如圖7所示。圖8給出了收斂歷程圖，可見每次流場自適應后流場保持了良好的收斂性。本算例網格數為151 233個，不考慮流場自適應的情況下，計算時間約為15 min。

圖9展示了6個展向位置自適應前后壓力系數的對比，可見流場自適應能顯著提高計算精度。同時該圖還給出了實驗數據[20]，可以看到計算結果和實驗數據吻合的較好。圖10給出了機翼表面和對稱面上的馬赫數云圖，本文的全速勢方法成功捕捉了機翼上表面的Lambda 型激波。

5.1.3 翼身融合布局無人機亞聲速氣動特性

將該方法應用于翼身融合布局(BWB)無人機的亞聲速計算仿真，Ma∞=0.1763,α=-6.28°～9.42°。實驗數據來源于南京航空航天大學進行的低速風洞實驗。圖11中展示該算例的計算網格，圖12給出了α=-0.072°時壓力系數云圖。圖13給出了升阻力特性曲線，可見計算結果和實驗數據之間誤差很小，比如升力系數最大誤差小于4%。算例展示該方法在非常規布局飛行器研究中的應用價值。

5.2歐拉方法

針對超聲速流動問題，采用基于自適應直角網格的歐拉方法，相比于全速勢方法，歐拉方法雖然計算速度相對較慢，但計算精度更高，尤其是對于激波較強的問題有更好的穩定性和精度。本節的算例來源于HL-20的風洞實驗[21]。網格如圖14所示，可見本文的直角網格方法能夠適應更加復雜的外形。圖15給出了自由來流Ma∞=1.2，α=0°時，流場自適應計算結果的馬赫數云圖。圖16給出了Ma∞=1.2和1.6時，不同迎角下升力系數計算結果和實驗數據的對比，可見該方法和實驗數據吻合度較高。網格數為412 730個，不考慮流場自適應的情況下，約需830迭代步和51 min的物理時間。

5.3基于當地表面斜度法的面元法

高超聲速問題第一個算例是HL-20在Ma∞=6不同迎角下的升阻比的對比。圖17(a)為迎角30°工況下的表面壓力系數分布；圖17(b)為升阻比隨迎角的變化情況，可見計算結果與實驗數據吻合良好，最大誤差(發生在迎角2.5°)不超過7%。網格數為15 320個，約需要20 s的物理時間。

第二個算例是AS-202的數值分析[4,22]。對比數據來源于CBAERO和DPLR對AS-202的仿真分析[4]。DPLR是高超聲速流場精確求解CFD程序。在自由來流8.24 km/s、Ma∞=28.6、迎角18.2°的工況下，圖18為返回艙的表面壓力系數分布；圖19為本文方法的數據與CBAERO、DPLR仿真數據分析對比。從圖19可見，z<0時，本文壓強分布略小于CBAERO和DPLR；z>0時，本文方法的壓強分布介于CBAERO和DPLR之間。

6 結 論

本文研究和開發了一種面向飛行器概念設計的全速域快速氣動分析工具。若干不同速度域和不同構型的算例表明：

1) 由參數化幾何建模程序、自適應直角網格方法程序、不同速度域的氣動分析程序組成的氣動分析工具，能快速地分析各種構型飛行器的氣動特性，計算過程均可完全自動化。

2) 計算精度較好。基于自適應直角網格的全速勢方程有限體積解法用于分析亞/跨聲速氣動問題，能夠實現基于流場解的網格自適應，提高了激波捕捉的精度；基于自適應直角網格的歐拉方程有限體積解法用于分析有強激波流動問題，彌補了全速勢方法的不足；當地表面斜度法的面元法用于高超聲速聲速問題，拓展了整套工具的適用范圍，其計算速度非常快，計算精度較好。

綜上所述，該工具能適用于全速域不同構型飛行器的氣動特性分析，具有速度快、計算精度較好、所需計算資源少、自動化程度高的特點，適用于飛行器概念設計階段的氣動特性分析。

[1]朱自強, 王曉璐, 吳宗成, 等. 民機設計中的多學科優化和數值模擬[J]. 航空學報, 2007, 28(1): 1-13.

[2]De Weck O, Agte J, Sobieszczanski-Sobieski J, et al. State-of-the art and future trends in multidisciplinary design optimization[R]. AIAA 2007-1905, 2007.

[3]Bonner E, Clever W, Dunn K. Aerodynamic preliminary analysis system II: Part I - theory[R]. NASA CR-182076, 1991.

[4]Kinney D J, Garcia J A, Huynh L. Predicted convective and radiative aerothermodynamic environments for various reentry vehicles using CBAERO[C]//44th AIAA aerospace Sciences Meeting and Exhibit. Doi: 10.2514/62006-659.

[5]Holst T L. Transonic flow computations using nonlinear potential methods[J]. Progress in Aerospace Sciences, 2000, 36(1): 1-61.

[6]Bolsunovsky A L, Buzoverya N P, Karas O V, et al. An experience in aerodynamic design of transport aircraft [C]//28th International Congress of the Aeronautical Sciences. ICAS, 2012.

[7]Rubbert P E, Bussoletti J E, Johnson F T, et al. A new approach to the solution of boundary value problems involving complex configurations[J]. Computational Mechanics-advances & Trends Amd, 1986: 49-84.

[8]Johnson F T, Tinoco E N, Yu N J. Thirty years of development and application of cfd at boeing commercial airplanes, seattle[J]. Computers and Fluids, 2005, 34(10): 1115-1151.

[9]De Zeeuw D L. A quadtree-based adaptively-refined Cartesian-grid algorithm for solution of the Euler equations[D]. Ann Arbor(MI): University of Michigan, 1993.

[10]Aftosmis M J, Berger M J, Melton J E. Robust and efficient Cartesian mesh generation for component-based geometry[J]. AIAA Journal, 1998, 6(36): 952-960.

[11]Hahn A. Vehicle sketch pad: A parametric geometry modeler for conceptual aircraft design[R]. AIAA 2010-657, 2010.

[12]Bonet J, Peraire J. An alternating digital tree(ADT)algorithm for 3D geometric searching and intersection problems[J]. International Journal for Numerical Methods in Engineering, 1991, 31(1): 1-17.

[13]Yang G, Causon D M, Ingram D M. Calculation of compressible flows about complex moving geometries using a three-dimensional Cartesian cut cell method[J]. International Journal for Numerical Methods in Fluids, 2000, 33(8): 1121-1151.

[14]Aftosmis M J, Berger M J, Adomavicius G. A parallel multilevel method for adaptively refined Cartesian grids with embedded boundaries[C]//38th Aerospace Sciences Meeting and Exhibit. Reno, Nevada, USA; 2000.

[15]Neel R E. Advances in computational fluid dynamics: Turbulent separated flows and transonic potential flows[D]. Blacksburg(VA): Virginia Polytechnic Institute and State University, 1997.

[16]張魯明. 航天飛機空氣動力學分析[M]. 國防工業出版社, 2009.

[17]Dejarnette F R, Ford C P, Young D E. Calculation of pressures on bodies at low angles of attack in supersonic flow[J]. Journal of Spacecraft and Rockets, 1980, 17(6): 529-536.

[18]呂凡熹, 肖天航, 余雄慶. 基于自適應直角網格的二維全速勢方程有限體積解法[J]. 計算力學學報, 2016, 33(03): 424-430.

[19]Redeker G. DLR-F4 wing-body configuration[R]. Report No.: AGARD-AR-303, 1994.

[20]Schmitt V, Charpin F. Pressure distributions on the ONERA-M6-wing at transonic mach numbers, experimental data base for computer program assessment[R]. Paris: The Fluid Dynamics Panel Working Group 04, 1979.

[21]Ware G M, Cruz C I. Aerodynamic characteristics of the HL-20[J]. Journal of Spacecraft and Rockets, 1993, 30(5): 529-536.

[22]Wright M J, Prabhu D K, Martinez E R. Analysis of apollo command module afterbody heating part I: AS-202[J]. Journal of Thermophysics and Heat Transfer, 2006, 20(1): 16-30.

Anaerodynamicanalysistoolforaircraftconceptualdesignatfullspeedrange

LYU Fanxi, LI Zhengzhou, DENG Jingshu, XIAO Tianhang, YU Xiongqing*

(CollegeofAerospaceEngineering,NanjingUniversityofAeronauticsandAstronautics,Nanjing210016,China)

For fast and automatic aerodynamic computation in unconventional aircraft conceptual design, a full-automatic aerodynamic analysis tool for arbitrary Mach number flow was developed in this paper. This tool consists of three methods for three individual Mach number ranges: 1) A finite volume full potential method on adaptive Cartesian grids for subsonic and transonic flow; 2) A finite volume Euler method on adaptive Cartesian grids for supersonic flow; 3) An empirical panel method based on local surface inclination methods for hypersonic flow. The feasibility and accuracy of the proposed method are validated by several typical cases of different flows around ONERA M6 wing, DLR-F4 wing-body, unmanned aerial vehicle (UAV) with a blended wing body (BWB), HL-20 lifting-body, and AS-202 capsule. The validation cases demonstrate a fast computation convergence with fully automatic grid treatment, and the results confirm the capacity of the present method in applications for unconventional aircraft conceptual design in a wide range of Mach numbers.

aircraft; aerodynamic; conceptual design; full potential equation; Euler equation; local surface inclination method; adaptation Cartesian grids

V211.3

A

10.7638/kqdlxxb-2016.0152

0258-1825(2017)05-0625-08

2016-11-30;

2017-03-06

國家自然科學基金(11672133);中央高校基礎科研業務費(NZ2016101)；江蘇高校優勢學科建設工程資助項目

呂凡熹(1989-),男，江西九江人，博士研究生，研究方向：飛機總體設計，計算流體力學. E-mail：lvfanxi@nuaa.edu.cn

余雄慶*，博士，教授. E-mail：yxq@nuaa.edu.cn

呂凡熹, 李正洲, 鄧經樞, 等. 面向飛行器概念設計的全速域氣動分析工具[J]. 空氣動力學學報, 2017, 35(5): 625-632.

10.7638/kqdlxxb-2016.0152 LYU F X, LI Z Z, DENG J S, et al. An aerodynamic analysis tool for aircraft conceptual design at full speed range[J]. Acta Aerodynamica Sinica, 2017, 35(5): 625-632.