基于Euler方程的返回艙氣動外形優化設計方法研究

李治宇，楊彥廣，袁先旭，唐志共

（中國空氣動力研究與發展中心，四川 綿陽 621000）

0 引 言

氣動外形的優化設計是飛行器設計的重要環節。由于飛行器方案論證階段氣動外形優化設計是通過數值優化方法和氣動性能計算的結合來實現的，其可靠性、精度和效率是關系到設計成敗的重要因素。而精度和效率主要由流場解算和優化算法決定。當前的工程算法求解速度快但精度較低，實用價值有限。在優化設計方面，氣動外形的優化設計一般采用已趨于成熟的梯度法等優化方法，但其最大的缺點在于設計結果容易陷入局部最優，并且對目標函數和約束的要求比較高。模擬退火算法等全局優化算法的計算量比較大，應用NS方程求解的時間長。而多目標優化設計經常選用的加權和方法，其權值不容易合理確定，并且不能很好地反應各優化目標之間的關系。

國內外研究者們針對各類飛行器外形進行了大量優化設計工作。印度的Rajesh Kumar Arora等人用牛頓近似的方法求解流場載荷，在法線方向引入壓力修正，通過多目標遺傳算法（Multi Objective Genetic Algorithm）優化設計了可重復使用的飛船返回艙外形，得到了滿足返回艙氣動性能并且可重復使用、成本較低的外形［1］。NASA的James L.Brown等人通過優化設計非對稱防熱底來提高升阻比、降低對流和輻射傳熱，其中氣動力計算部分調用了高保真度工程算法的CBAero軟件包，用修正牛頓理論方法計算表面壓力分布［2］。Alan Le Moigne和 Wai Sam Wong等人通過求解N－S方程計算翼型的氣動性能，并應用序列二次規劃（NLPQL）方法進行優化設計。該方法能夠快速得到優化解，精度較高，但卻容易陷入局部最優的境地［3－4］。唐偉、車競、夏露等人應用內伏牛頓理論及定常勢流格林函數法計算飛行器氣動性能，并且通過遺傳算法優化設計了新外形的氣動布局［5－9］。

本文基于iSIGHT優化軟件平臺，應用遺傳算法和數值優化算法相結合的混合算法，通過應用快速求解Euler方程的CART3D程序，發展了一種快速有效的氣動布局優化設計方法，并以CTV返回艙為例驗證該方法，將其應用于類神舟飛船返回艙的優化設計中得到可靠結果。一般說來，該方法可用于飛行器氣動布局的優化設計。

1 研究方法

1.1 研究方法簡介

通過SolidWorks軟件作圖得到基本氣動布局外形并完成參數化建模。在優化設計過程中用到的參數化建模方法適用于SolidWorks和CFD分析軟件。本文應用一個VBS腳本文件來完成外形的參數化，該VBS文件通過直接修改SolidWorks文件的特征尺寸來修改外形，并輸出parasolid格式的數模文件。在輸入錯誤尺寸數據時，如負值，該文件并不會中斷程序，而是在輸出原始外形后繼續優化過程，優化過程忽略錯誤數據。所以，在出現錯誤尺寸時該文件不會導致優化過程中斷，也不會影響到優化結果，具有很好的魯棒性。

應用快速氣動分析軟件CART3D計算外形的氣動性能。該軟件通過快速生成非結構笛卡爾網格、求解Euler方程的方法計算高超聲速氣動力特性。在第一次運算中記錄操作命令流文件，之后的氣動計算就可以直接通過所生成的命令流文件來運行CART3D程序，期間不需要人工干預操作。并且，iSIGHT可直接調用該文件并修改其中參數。該軟件精度高、求解速度快、網格生成高度自動化、可用于機翼、全機、導彈及噴管等復雜系統的氣動力計算分析，適用性很好［10］。

優化計算過程選擇iSIGHT優化平臺，調用合適的優化算法，完成優化設計部分的工作，最終得到優化結果。iSIGHT平臺將大量需要人工完成的工作實現自動化處理，并且易于處理復雜問題的優化設計過程，可采用多種方法探索，可縮短設計周期、降低設計成本，提高設計質量和可靠性。iSIGHT提供了多種商業軟件及自編程序接口，設計者通過圖形用戶界面GUI可實現設計問題的過程集成、問題表述、優化方案選擇及求解監控等，iSIGHT提供多種優化算法，通過組合或單獨使用可基本滿足設計者要求。

集成的整個優化過程包括四部分，分別為外形特征尺寸計算（calculation0）、參數化生成外形（SW）、網格生成及氣動力計算（CART）、目標函數計算（calculation1）。如圖1所示。

圖1 優化流程圖Fig.1 Flow chart of the optimization

1.2 流場求解方法驗證

1.2.1 解算精度驗證

CART3D求解器基于可壓縮Euler方程，采用有限體積法離散，流動變量位于網格中心，時間推進采用Runge－Kutta法得到穩態解，空間離散為迎風格式，選擇使用限制器，格式具有TVD性質，采用多重網格法加速迭代過程。

通過CART3D計算得到飛船返回艙的高超聲速氣動力系數，并與文獻數據對比，驗證計算方法的精度。計算來流馬赫數M＝6，攻角α＝20°，圖2為計算網格和流場壓力分布圖，表1為計算結果對比，表中數據誤差為CART3D計算結果與各參考文獻數據對比的相對誤差，結果顯示該方法計算精度是可以信賴的。其計算質量主要取決于所生成的笛卡爾網格的大小，可以通過調整網格加以改進。圖3為CART3D和FLY3D計算的飛船返回艙升阻比和對頂點的俯仰力矩系數隨攻角變化的曲線及相應實驗結果，可見其計算結果符合很好。

圖2 計算網格與流場壓力分布圖Fig.2 Cartesian grid and flow－fields pressure distribution

表1 計算結果與文獻數據對比Table 1 The numerical results vs.the reference

圖3 返回艙升阻比和俯仰力矩隨攻角變化的曲線Fig.3 Comparison of aerodynamic coefficient CL／CDand Cmvs.angles of attack

1.2.2 解算效率驗證

通過與FLY3D程序對比發現，FLY3D網格約190萬時，在PC機上運行一次的計算機時約6h；而當CART3D網格量為130萬時，在PC機上運行一次的時間小于10min。通過表2對比，機時與網格量有類似正比關系，而計算結果的精度在一定范圍內變化不大。這樣，我們就可以在保證計算精度的前提下適當減少網格，以縮短機時，使得CART3D的計算效率滿足優化設計多輪計算的要求。在優化過程中CART3D實際網格量約在20萬左右，在PC機上運行一次的時間約為70～90s，滿足優化設計要求。表2中計算條件與表1相同。另外，其網格生成高度自動化，僅需數秒時間，效率很高。

表2 不同網格量的CART3D運行時間與精度對比表Table 2 Comparison of run time and precision vs.grid quantity of CART3D

2 優化設計算例驗證

2.1 計算模型

［13］，通過完成歐洲航天局為其載人天地往返運輸系統計劃中宇航員輸送艙（Cosmonaut Transportation Vehicle，以下簡稱CTV）所提出的氣動構型CTV的優化設計來對比選擇合適的優化算法。CTV構型為鈍球雙錐外形，其基本構型特征尺寸及輪廓如圖4所示。

圖4 CTV返回艙輪廓線Fig.4 Contour line of CTV reentry capsule

2.2 計算目標及約束

為確定CTV的最優氣動布局，這里采用以下簡化模型。在底部直徑D和總體積V不變的情況下，滿足配平升阻比CL／CD不小于0.8的約束，設計目標為最大的體積利用率RV。設計變量為球頭半徑RN、前錐長LF、前錐半錐角θF?？杀硎緸椋?/p>

設計目標： minU＝1／RV

約束條件：D＝3.5m

V＝3m3

CL／CD＞0.8

設計變量： 0.5m＜RN＜1.0m，2.0m＜LF＜3.0m，10.0°＜θF＜30.0°

其中，體積利用率RV＝6／S3／2，V為體積，S為表面積。設計點來流馬赫數M＝6，配平攻角αT＝20°。

2.3 算法選擇

本文通過不同算法進行優化對比。選用算法包括數值優化算法：MMFD（Modified Method of Feasible Directions－修正可行方向算法）、NLPQL（序列二次規劃算法）、LSGRG（Generalized Reduced Gradient－廣義簡約下降梯度法）及全局優化算法MIGA（Multi－Island Genetic Algorithm－多島遺傳算法）。

數值優化算法通過梯度搜索，處理不等式和等式約束，快速得到最優解；MIGA通過選擇、交叉及變異等遺傳操作尋優，把種群分組提高全局搜索能力。

2.4 結果分析

所選用的不同算法所得到的優化解基本達到了優化目的，即在滿足約束的條件下，優化目標體積利用率RV達到最大值，并且除算法2因為在局部最優點中止優化外，其余算法得到的結果，都略優于文獻數據。對比結果如表3所示，表中長度單位為米，角度單位為度。

由表3中算法1與2對比可見，MMFD算法得到的是局部最優解。對比算法3、4、5發現，最優值接近。由優化數據可以得到，算法3中的NLPQL實際上沒有起到優化作用，即該算法的優化部分沒有得到更優解。算法4的MMFD及算法5的LSGRG的運行對結果進行了更進一步的局部優化?？梢哉f明的是MIGA在不同初始點開始計算，得到了近似的全局最優點。而從算法6中可以看到，將MIGA算法的進化代數由10改為12后，得到了“更優解”。所得最優解的氣動力系數如表4所示，與文獻［14］中配平攻角20°時最大升阻比為0.8左右的計算結論及FLY3D精細計算的結果相吻合。

表3 不同優化算法的優化結果對比Table 3 Comparison between different optimization techniques

表4 算法6得到最優解的氣動力系數Table 4 Aerodynamic coefficient of the optimization result for algorithms 6

通過以上分析，可以確定適合于氣動布局外形優化的優化算法。通過多島遺傳算法MIGA全局優化，得到全局最優解后再通過數值優化算法進行優化計算，得到局部更優解，而算法MIGA的參數設置將會影響到結果的優化程度。這樣的混合優化算法在避免陷入局部最優解的同時，數值優化算法的優化計算對全局優化算法得到的最優解進行“微調”，得到“更優解”。

3 類神舟飛船返回艙優化設計

3.1 單目標優化設計

類神舟飛船返回艙外形的主要特征參數包括最大截面直徑dm、球冠半徑Rn、大頭肩部半徑Rc、倒錐角θ、長度L、尾部圓弧半徑RB及尾部圓臺的底面半徑和高度。如圖5所示。

圖5 類神舟飛船返回艙輪廓線Fig.5 Contour line of Shenzhou Spaceship reentry capsule

將各特征參數，除倒錐角外，用dm無量綱化。球冠鈍度大頭肩部圓弧曲率半徑Rc／dm、尾部圓弧曲率半徑長細比L／dm。優化要求在體積不小于16m3、體積利用率不小于0.85時，通過變化最大截面直徑、球冠鈍度、大頭肩部曲率半徑及倒錐角四個變量得到升阻比最大的優化外形，計算來流M＝13.2，αT＝20.0°，優化計算模型表示為：

目標： maxCL／CD

約束：V≥16m3

Rv≥0.85

變量： 2.5m≤dm≤4.0m

選用MIGA和MMFD相結合的優化算法，其中的多島遺傳算法設置島數為10，每島10個個體，進化12代，染色體長度20，交叉概率0.9，變異概率0.01，遷徙率0.4，遷徙間隙3代。

3.2 雙目標優化設計

雙目標優化設計與以上單目標優化設計相比，約束條件及變量均相同，不同之處在于設計目標要求在升阻比最大的同時峰值熱流最小化，其中峰值熱流的計算采用Fay－Riddell公式。優化計算模型可表示為：

表5 單目標優化設計結果Table 5 Optimization result for single objective

目標： maxCL／CD，minqw

約束：V≥16m3

Rv≥0.85

變量： 2.5m≤dm≤4.0m

選用鄰域培植遺傳算法NCGA，種群大小20，進化150代，最后輸出pareto解集。兩個目標的比例因子分別取為1.0和1000.0，權系數均為1.0，則目標函數objective＝qw／1000.0－CL／CD／1.0。相應優化結果及優化外形如表6、圖6所示。

表6 雙目標優化設計結果Table 6 Optimization result for double objectives

圖6 類神舟飛船返回艙優化pareto最優解集及最優化外形圖Fig.6 Pareto plot and the optimization contour of Shenzhou Spaceship reentry capsule

4 小 結

（1）形成了基于iSIGHT軟件平臺的集成優化方法，用于返回艙優化設計能夠得到理想結果。其優化精度高、速度快，有較好的工程應用前景。

（2）所采用的流場解算方法適用于多種復雜系統的氣動力計算，適用性強，計算精度和效率較高。

（3）選擇MIGA與數值優化算法相結合的混合遺傳算法，可以滿足單目標多變量優化設計。NCGA同樣可以滿足多目標優化設計，影響優化結果的主要因素是遺傳算法選擇、交叉及變異等相關參數的選擇。下一步的研究工作將集中在優化算法的選擇應用及相應優化參數對優化結果的影響方面。

參考文獻:

［1］RAJESH KUMAR ARORA.Aerodynamic shape optimization of a Re－entry capsule［R］.AIAA 2003－5394.

［2］JAMES L BROWN，JOSEPH A GARCIA.An asymmetric capsule vehicle geometry study for CEV［R］.AIAA 2007－604.

［3］ALAN LE MOIGNE.A discrete Navier－Stokes adjoint method for aerodynamic optimisation of Blended Wing－Body configurationsp［D］.［PHD Thesis］.Cranfield U－niversity，2002.

［4］WAI SAM WONG.Mechanisms and optimisations of 3Dshock control bumps［D］.［PHD Thesis］.Sheffield University，2006.

［5］唐偉，桂業偉，張勇，等.未來升力再入飛船返回艙的氣動布局設計［C］.中國宇航學會首屆學術年會.北海，2005.

［6］唐偉，高曉成，李為吉，等.雙橢圓截面再入飛行器的氣動計算及布局優化設計［J］.空氣動力學學報，2004，22（2）：171－174.

［7］唐偉，桂業偉，方方.新型升力再入飛船返回艙氣動外形選型研究［J］.宇航學報，2008，29（1）：84－88.

［8］車競，唐碩，何開鋒，等.高超聲速飛行器氣動布局總體性能優化設計研究［J］.空氣動力學學報，2009，27（2）：214－219.

［9］夏露，高正紅，李天.飛行器外形多目標多學科綜合優化設計方法研究［J］.空氣動力學學報，2003，23（1）：275－280.

［10］JAMES S GREATHOUSE，BENJAMIN S KIRK，RANDOLPH P LILLARD.Crew exploration vehicle（CEV）crew module shape selection analysis and CEV aeroscience PROJECT overview［R］.AIAA 2007－603.

［11］袁先旭.非定常流動數值模擬及飛行器動態特性分析研究［J］.［博士學位論文］.中國空氣動力研究與發展中心，2002.

［12］紀楚群，周偉江.鐘形返回艙空氣動力特性數值模擬［J］.航天返回與遙感，2001，22（1）：32－37.

［13］安偉剛.多目標優化方法研究及其工程應用［D］.［博士學位論文］.西北工業大學，2005.

［14］唐偉，張勇，李為吉，等.歐洲返回艙CTV氣動特性綜述［J］.空氣動力學學報，2003，21（3）：362－366.