馬洋 榮偉 江長虹 張青斌 霍霖

（1國防科學技術(shù)大學航天科學與工程學院，長沙 410073）

（2北京空間機電研究所，北京 100094）

1 引言

采用球冠倒錐外形的返回艙在再入過程中將經(jīng)歷十分嚴酷而復雜的氣動力熱環(huán)境。球冠表面受氣流強烈壓縮，承受著高溫、高壓氣流，氣流在球冠和倒錐聯(lián)結(jié)的肩部附近發(fā)生劇烈膨脹，在后體錐面上形成高速低壓的流動。精心設(shè)計返回艙外形可以有效地減少嚴酷力熱環(huán)境對其再入過程產(chǎn)生的不利影響，也是保證成功完成再入任務的關(guān)鍵點之一。返回艙在大氣層內(nèi)利用空氣阻力減阻，要獲得盡可能大的阻力以盡快減小飛行速度可以通過外形設(shè)計達到目的，但阻力的增大會受到諸多因素的制約。首先由于阻力與升阻比存在著天然的“矛盾關(guān)系”，而升阻比又是控制再入彈道、改進落點精度最關(guān)鍵的因素，因而在增大阻力時必須考慮到升阻比的要求。另外在高速再入過程中產(chǎn)生的氣動加熱也會對返回艙外形設(shè)計提出特殊要求。

氣動特性的精確預示是進行返回艙外形優(yōu)化設(shè)計的基礎(chǔ)，陳河梧采用試驗手段研究了 Ma=5~8，攻角在–27°~2°條件下，類神舟飛船返回艙的高超聲速氣動特性，認為返回艙球冠主導小攻角的氣動力變化，大攻角時倒錐的氣動力貢獻加大，升力和俯仰力矩出現(xiàn)比較明顯的非線性增量[1]。紀楚群和呂俊明采用計算流體力學(Computational Fluid Dynam ics, CFD)數(shù)值仿真手段對鐘形返回艙和火星再入返回艙進行數(shù)值模擬，對比分析了基于Euler方程和Navier-Stokes(NS)方程的CFD模擬結(jié)果，并分析了化學非平衡效應的影響，研究表明，在高超聲速情況下，由于物體迎風面的壓力遠遠大于背風面壓力，背風面壓力計算誤差對整體氣動力特性影響很小，同時氣流在物面的分離情況不嚴重，因此用 Euler方程可足夠準確地預測返回艙高超聲速情況的流場及氣動特性[2-3]?；瘜W非平衡效應對返回艙氣動特性影響較小，但它會對激波形狀和駐點位置產(chǎn)生影響。返回艙外形優(yōu)化設(shè)計研究方面，李治宇采用Euler方法在ISIGHT集成軟件環(huán)境下對返回艙外形展開優(yōu)化設(shè)計，得到了升阻比和駐點熱流的優(yōu)化前緣[4]。John E Theisinger對火星再入返回艙外形進行了詳細的氣動力熱分析，通過對外形進行優(yōu)化分析，對阻力、配平狀態(tài)下的穩(wěn)定性、表面熱流、駐點熱流、容積率、容積、升阻比等各種相互制約的性能進行協(xié)調(diào)[5-6]。

本文氣動特性分析采用三維 NS方程，不考慮真實氣體的化學非平衡效應，引入代理模型以減少計算開銷，采用Fay-Riddell公式計算駐點熱流，并利用成熟的改進非劣分類遺傳算法(Nondom inated Sorting Genetic A lgorithm-II, NSGA-II)[7]，對參數(shù)化設(shè)計的返回艙外形進行多目標優(yōu)化設(shè)計，該優(yōu)化設(shè)計方法在返回艙方案設(shè)計階段具有一定的應用價值。

2 研究方法

在氣動外形優(yōu)化設(shè)計中，設(shè)計結(jié)果的精度和計算效率往往相互沖突。結(jié)果的精度主要由氣動特性計算精度和尋優(yōu)算法精度決定，本文選擇經(jīng)典遺傳算法尋優(yōu)，尋優(yōu)算法精度和效率不在本文的研究之列。氣動特性計算方法主要有工程估算方法、CFD數(shù)值仿真以及兩種方法相結(jié)合的混合方法等。工程估算效率最高，其精度也最差，較少直接用于優(yōu)化設(shè)計。CFD方法計算精度很高，許多情況下可以替代試驗研究，但在飛行器氣動外形優(yōu)化設(shè)計中，需要反復修改氣動外形，多次進行氣動特性計算，采用CFD方法的計算成本即使是對最先進的計算資源來說仍顯得力不從心?；旌戏椒ㄓ捎诰哂泄こ坦浪愀咝Ш虲FD仿真高精度的潛質(zhì)，在優(yōu)化設(shè)計中正受到一些研究人員的重視[8-9]，但其計算精度受CFD結(jié)果和所采取的擬合算法的限制，其應用還需要進一步的研究。

本文引入代理模型技術(shù)以減小CFD計算的開銷，并采用改進的EI（Expected Improvement）函數(shù)加點策略[10]，充分發(fā)揮代理模型的高效性，并以高效的加點策略保證結(jié)果的精度，很好地解決了氣動外形優(yōu)化設(shè)計中精度與效率的問題，具體優(yōu)化設(shè)計方法如圖1所示。首先采用試驗設(shè)計方法產(chǎn)生構(gòu)建代理模型的樣本點，然后采用CFD方法計算樣本點的氣動特性，并基于計算結(jié)果構(gòu)建代理模型，再判斷代理模型計算結(jié)果與CFD計算結(jié)果的誤差是否滿足精度要求，如果不滿足則按加點策略加點重新進行CFD計算并構(gòu)建代理模型，直至滿足收斂要求，最后用構(gòu)建好的代理模型代替CFD仿真進行外形優(yōu)化設(shè)計。

圖1 優(yōu)化設(shè)計流程Fig.1 Optimization design process

2.1 氣動特性分析方法及驗證

氣動性能計算采用成熟的Fluent軟件。選取可實現(xiàn)的k-ε兩方程湍流模型，并配合使用非平衡壁面函數(shù)，這樣更適合于模擬返回艙附近擾流流場中出現(xiàn)的大分離和大漩渦特性，對流項離散采用二階AUSM格式，返回艙表面滿足無滑移邊界條件，進口取來流參數(shù)，出口數(shù)值邊界條件采用外推方式獲得。

采用結(jié)構(gòu)網(wǎng)格離散流場，由于流場的對稱性，只對一半流場進行網(wǎng)格劃分，在靠近返回艙表面對網(wǎng)格進行加密，返回艙和用于數(shù)值精度驗證的計算網(wǎng)格如圖2和圖3所示。

圖2 返回艙計算網(wǎng)格Fig.2 Computing grid of reentry capsul

圖3 數(shù)值驗證算例計算網(wǎng)格Fig.3 Computing grid of CFD validation

為驗證本文CFD計算模型的準確性，對文獻[11]提供的火星探測實驗飛行器外形進行數(shù)值模擬，并與文獻[3]和[11]計算結(jié)果進行對比（如表1所示，其中α、CL、CD、L/D和CM分別表示攻角，升、阻力系數(shù)，升阻比和俯仰力矩系數(shù)），可以看出，本文 CFD計算得到的氣動特性與文獻結(jié)果的誤差約為2~9%，數(shù)值仿真模型具有較高的可信度。

表1 火星探測實驗飛行器氣動特性計算驗證Tab.1 Aerodynam ic computing accuracy verification of Mars Science Laboratory

2.2 返回艙駐點熱流計算

由于熱流數(shù)值仿真需要比氣動力數(shù)值模擬更多的計算資源，而工程經(jīng)驗公式能在趨勢上給出氣動熱與其它性能的大致協(xié)調(diào)關(guān)系，因此在方案設(shè)計階段采用經(jīng)驗公式計算駐點熱流不失為一種折衷的選擇。返回艙駐點熱流采用簡化的Fay-Riddell公式[12]：

2.3 代理模型構(gòu)建

代理模型技術(shù)主要包括多項式響應面(Polynom ial Response Surface)、徑向基函數(shù)(Radial Basis Function)、神經(jīng)網(wǎng)絡(artificial neural networks)和Kriging方法等，前兩種方法在處理非線性強和高維數(shù)問題時表現(xiàn)不好，神經(jīng)網(wǎng)絡方法雖然較適合處理非線性強和高維數(shù)問題，但其預測精度和魯棒性在一定程度上不及Kriging方法[13]。

本文采用Kriging方法構(gòu)建替代CFD分析的代理模型[14]，為了增強代理模型的尋優(yōu)精度和效率，借鑒經(jīng)典的EI函數(shù)加點方法[10]并加以改進，提出在現(xiàn)有樣本點基礎(chǔ)上，同時加入EI值最大的點和Kriging模型預測方差最大的點來改善模型預測精度，EI值和 Kriging模型預測方差的具體計算方法參考文獻[10,14]，這樣的加點尋優(yōu)方法一方面具有EI函數(shù)方法兼顧局部和全局尋優(yōu)的特點，另一方面在模型預測方差最大處加點更注重模型的全局精度，并且多點加點方式也適合于并行計算，有利于提高優(yōu)化設(shè)計的效率。

2.4 尋優(yōu)算法

遺傳算法模擬生物在自然環(huán)境中的遺傳和進化過程，是一種自適應、全局優(yōu)化、概率搜索算法，它具有良好的魯棒性、全局性和并行性，在當今工程技術(shù)領(lǐng)域的多目標優(yōu)化問題中應用非常廣泛。在眾多的多目標遺傳算法中，NSGA采用簡潔明晰的非優(yōu)超排序機制，使算法具有逼近Pareto最優(yōu)解的能力，采用排擠機制保證得到的Pareto最優(yōu)解具有良好的分布，Deb在NSGA基礎(chǔ)上改進得到NSGA-II[15]，引入精英保留策略對產(chǎn)生的非劣解進行非劣快速分類和密度估計操作，比較其擁擠度，確定是否可接受為Pareto最優(yōu)解。NSGA2算法提出了新的基于分級的快速非劣排序算法，大大降低了計算復雜度；同時為使Pareto最優(yōu)解散布范圍較大，盡可能均勻遍布，引入了擁擠距離的概念，并采用擁擠比較算子代替需要計算共享參數(shù)的適應度共享方法；最后引入精英保留機制，不但擴大了采樣空間，而且有利于保持優(yōu)良個體，迅速提高了種群的整體水平。由于NSGA-II算法的上述優(yōu)點和可靠性，本文以其作為優(yōu)化算法。

3 多目標優(yōu)化與分析

3.1 返回艙外形描述

返回艙外形如圖4所示，可以采用最大直徑dm、最大高度H、端框直徑dD和高度HD、大底半徑RN、前體半徑RF、大底倒圓半徑 Rc和后體傾角θ來表示。

3.2 問題描述

在返回艙再入大氣的過程中，需要依靠大氣阻力使其減速才能完成再入返回的任務，同時返回艙也需要盡可能大的升阻比，這樣能實現(xiàn)橫向和縱向的機動，得到更寬的再入走廊，從而改善落點的精度。但是阻力和升阻比從來就存在著“矛盾”，要想得到滿足再入返回要求的返回艙外形必須同時考慮這兩個目標，進行多目標優(yōu)化。另外在返回艙高速再入過程中，嚴重的氣動加熱是必須考慮的問題，同時為滿足載人以及設(shè)備存放的要求，返回艙內(nèi)的空間也是設(shè)計需要關(guān)注的地方，本文將返回艙表面最大熱流（駐點熱流）和容積率作為約束條件加入到優(yōu)化設(shè)計過程中，以體現(xiàn)多種約束條件下的設(shè)計。

在圖4所示的8個外形參數(shù)中，最大高度H和最大直徑dm受返回艙總體規(guī)模的限制變化較小，端框直徑dD和高度HD由于處在背風區(qū)對返回艙的氣動力和氣動熱特性影響較小，對容積特性的影響也較小，它們都不適宜作為本文優(yōu)化問題的設(shè)計變量。后體傾角θ雖然可變范圍較大且對返回艙性能有較大影響，但在大底半徑RN、前體半徑RF和大底倒圓半徑Rc確定的前提下，θ是確定的，即θ、RN、RF和Rc只有是3個參數(shù)是獨立的。本文選取RN、RF和Rc作為設(shè)計變量，它們的取值范圍如表2。

圖4 返回艙外形參數(shù)Fig.4 Configuration parameters of reentry capsule

表2 設(shè)計變量及其變化范圍Tab.2 Design variables and their variation range

在某一典型的再入條件下[1]，返回艙外形優(yōu)化設(shè)計問題可以這樣描述：

3.3 優(yōu)化結(jié)果

多目標優(yōu)化結(jié)果如圖5所示。圖中空心圓圈為優(yōu)化得到的滿足約束條件的前緣點，它們分布均勻，滿足約束條件，并體現(xiàn)了兩個優(yōu)化目標間的沖突和妥協(xié)，設(shè)計者可以根據(jù)返回艙再入飛行的實際需求在前緣點上選擇相應的氣動外形；實心圓點為同時滿足 CD和 L/D最大的理想狀態(tài)（圖中標注“Ideal”），由于優(yōu)化目標相互沖突，該理想狀態(tài)是達不到的；“＋”標記點表示用于構(gòu)建代理模型的樣本點；前緣上首尾標注“Max CD”和“Max L/D”的點為滿足約束前提下，分別以CD和L/D為目標的單目標優(yōu)化結(jié)果；圖中標注“TPF”的點對應于某一典型的優(yōu)化前緣(Typical Pareto Front，TPF)。圖5還給出了樣本點和前緣點上典型狀態(tài)對應的返回艙外形。

圖5 多目標優(yōu)化結(jié)果Fig.5 Multi-objective optim ization results

3.4 優(yōu)化設(shè)計精度與效率分析

在構(gòu)建代理模型時初始樣本點為30個，經(jīng)過3次加點迭代和模型修正后，Kriging模型最大預測方差滿足收斂要求，然后隨機產(chǎn)生4個測試點對建立的代理模型進行驗證，升阻比、阻力系數(shù)和容積率的最大誤差分別為6.52%、3.58%和1.24%。另外在優(yōu)化前緣上選取3個標注代號為“Max CD”、“Max L/D”和“TPF”的外形，采用CFD方法計算對應的氣動特性，與采用代理模型預測的結(jié)果比較，表3為阻力系數(shù)和升阻比的比較結(jié)果，可見，CFD結(jié)果與代理模型響應差別最大不超過7%。由此可見基于本文介紹的方法構(gòu)建的代理模型具有較好的精度，可以替代真實物理模型進行優(yōu)化設(shè)計。

表3 典型前緣點特性驗證Tab.3 Aerodynam ic character verification of typical optim ization front

本文在優(yōu)化過程中，總共調(diào)用代理模型20 051次，這意味著如果氣動分析直接采用CFD方法的話，總共需要進行20 051次CFD數(shù)值計算，而采用代理模型替代數(shù)值仿真后，代理模型的計算時間與CFD計算時間相比幾乎可以忽略不計，整個優(yōu)化過程的計算時間主要體現(xiàn)在構(gòu)建代理模型時的40次CFD計算時間上，可見引入代理模型后，大約可以減少20 000次CFD分析時間，這樣高效的分析方法十分適合返回艙總體優(yōu)化設(shè)計。

4 結(jié)束語

通過對返回艙進行氣動外形優(yōu)化設(shè)計，得到了如下結(jié)論：

1）引入代理模型能極大提高氣動特性的計算效率，同時采用改進的EI加點策略能將樣本點的數(shù)目大大減少而不降低代理模型的精度，優(yōu)化計算結(jié)果表明，代理模型計算結(jié)果與CFD計算結(jié)果誤差可以控制在7%以內(nèi)，基于代理模型技術(shù)的優(yōu)化過程可以節(jié)省絕大部分（95%以上）的計算開銷；

2）以返回艙升阻比和阻力系數(shù)為目標，以駐點熱流和容積率為約束條件進行多目標多約束優(yōu)化，得到了優(yōu)化前緣，為返回艙總體設(shè)計提供了有益的借鑒。

References)

[1]陳河梧. 飛船返回艙高超聲速氣動特性的風洞實驗分析 [J]. 航天器工程, 2008, 17(5): 77-81.CHEN Hewu. Wind-tunnel Test Analysis on Hypersonic Aerodynam ic Characteristics of Returnable Module [J]. Spacecraft Engineering, 2008, 17(5): 77-81. (in Chinese)

[2]紀楚群, 周偉江. 鐘形返回艙空氣動力特性數(shù)值模擬 [J]. 航天返回與遙感, 2001, 22(1): 33-37.JI Chuqun, ZHOU Weijiang. A Numerical Simulation of Aerodynamic Characters for Bell Capsule [J]. Spacecraft Recovery and Remote Sensing, 2001, 22(1): 33-37. (in Chinese)

[3]呂俊明, 程曉麗, 王強. 火星科學實驗室氣動特性數(shù)值分析 [J]. 力學與實踐, 2013, 35(1): 31-35.LV Junm ing, CHENG Xiaoli, WANG Qiang. Numerical Aerodynam ic Analysis of Mars Science Laboratory [J]. Mechanics in Enginering, 2013, 35(1): 31-35. (in Chinese)

[4]李治宇, 楊彥廣, 袁先旭. 基于 Euler方程的返回艙氣動外形優(yōu)化設(shè)計方法研究 [J]. 空氣動力學學報, 2012, 30(5):653-657.LI Zhiyu, YANG Yanguang, YUAN Xianxu. The Study of the Reentry Capsule Shape Optimization Method Based on the Solving of the Euler Equations. [J]. Acta Aerodynam ic Sinica, 2012, 30(5): 653-657. (in Chinese)

[5]John E T, Robert D B. Multi-Objective Hypersonic Entry Aeroshell Shape Optim ization [J]. Journal of Spacecraft and Rockets, 2009, 46(5): 957-966.

[6]John E T, Robert D B. Aerothermodynam ic Shape Optimization of Hypersonic Entry Aeroshells [C]. 13th AIAA/ISSMO Multidisciplinary Analysis Optimization Conference, Fort Worth, Texas, 13-15 September, 2010.

[7]Deb K. Multi-Objective Optimization Using Evolutionary Algorithm [M]. Chichester: John W: ley & Sons, 2002.

[8]Slawomir K, Leifur L. Transonic Airfoil Shape Optim ization Using Variable-Resolution Models and Pressure Distribution Alignment [C]. 29th AIAA Applied Aerodynam ics Conference, Honolulu, Hawaii, 2011.

[9]Li H, Lin P, Yang Z. Modeling and Design Optimization of a Common Aero Vehicle w ith Parameterized Configuration [C].AIAA Modeling and Simulation Technologies Conference, M inneapolis, M innesota, 2012.

[10]Donald R J, Matthias S, William J W. Efficient Global Optim ization of Expensive Black-Box Functions [J]. Journal of Global Optimization, 1998, 13: 455-492.

[11]Artem A D, Edquist K, Shoenenberger M. Influence of the Angle of Attack on the Aerothermodynam ic Environment of the Mars Science Laboratory [J]. AIAA paper, 2006, 2006-3889.

[12]車競, 唐碩, 何開鋒. 類乘波體飛行器氣動加熱的工程計算方法 [J]. 彈道學報, 2006, 18(4): 93-96.CHE Jing, TANG Shuo, HE Kaifeng. Aerothemo Engineering Method for Quasi-Waverider Vehicle [J]. Journal of Ballistics,2006, 18(4): 93-96. (in Chinese)

[13]Ricardo M P, André R D C, Curran C. Comparison of Surrogate Models in a Multidisciplinary Optimization Framework for Wing Design [J]. AIAA Journal, 2010, 48(5): 995-1006.

[14]Lophaven S N, Nielsen H B, Sondergaard J. DACE a Matlab Kriging Toolbox [R]. Technical University of Denmark, 2002.

[15]Deb K, Pratap A, Agrawal S. A Fast and Elitist Multi-objective Genetic Algorithm: NSGA-II [R]. IEEE Transactions on Evolutionary Computation, 2000, 6(2): 182-197.