基于Euler方程的返回艙氣動(dòng)外形優(yōu)化設(shè)計(jì)方法研究

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

（中國空氣動(dòng)力研究與發(fā)展中心，四川 綿陽 621000）

0 引 言

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

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

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

1 研究方法

1.1 研究方法簡介

通過SolidWorks軟件作圖得到基本氣動(dòng)布局外形并完成參數(shù)化建模。在優(yōu)化設(shè)計(jì)過程中用到的參數(shù)化建模方法適用于SolidWorks和CFD分析軟件。本文應(yīng)用一個(gè)VBS腳本文件來完成外形的參數(shù)化，該VBS文件通過直接修改SolidWorks文件的特征尺寸來修改外形，并輸出parasolid格式的數(shù)模文件。在輸入錯(cuò)誤尺寸數(shù)據(jù)時(shí)，如負(fù)值，該文件并不會(huì)中斷程序，而是在輸出原始外形后繼續(xù)優(yōu)化過程，優(yōu)化過程忽略錯(cuò)誤數(shù)據(jù)。所以，在出現(xiàn)錯(cuò)誤尺寸時(shí)該文件不會(huì)導(dǎo)致優(yōu)化過程中斷，也不會(huì)影響到優(yōu)化結(jié)果，具有很好的魯棒性。

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

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

集成的整個(gè)優(yōu)化過程包括四部分，分別為外形特征尺寸計(jì)算（calculation0）、參數(shù)化生成外形（SW）、網(wǎng)格生成及氣動(dòng)力計(jì)算（CART）、目標(biāo)函數(shù)計(jì)算（calculation1）。如圖1所示。

圖1 優(yōu)化流程圖Fig.1 Flow chart of the optimization

1.2 流場求解方法驗(yàn)證

1.2.1 解算精度驗(yàn)證

CART3D求解器基于可壓縮Euler方程，采用有限體積法離散，流動(dòng)變量位于網(wǎng)格中心，時(shí)間推進(jìn)采用Runge－Kutta法得到穩(wěn)態(tài)解，空間離散為迎風(fēng)格式，選擇使用限制器，格式具有TVD性質(zhì)，采用多重網(wǎng)格法加速迭代過程。

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

圖2 計(jì)算網(wǎng)格與流場壓力分布圖Fig.2 Cartesian grid and flow－fields pressure distribution

表1 計(jì)算結(jié)果與文獻(xiàn)數(shù)據(jù)對比Table 1 The numerical results vs.the reference

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

1.2.2 解算效率驗(yàn)證

通過與FLY3D程序?qū)Ρ劝l(fā)現(xiàn)，F(xiàn)LY3D網(wǎng)格約190萬時(shí)，在PC機(jī)上運(yùn)行一次的計(jì)算機(jī)時(shí)約6h；而當(dāng)CART3D網(wǎng)格量為130萬時(shí)，在PC機(jī)上運(yùn)行一次的時(shí)間小于10min。通過表2對比，機(jī)時(shí)與網(wǎng)格量有類似正比關(guān)系，而計(jì)算結(jié)果的精度在一定范圍內(nèi)變化不大。這樣，我們就可以在保證計(jì)算精度的前提下適當(dāng)減少網(wǎng)格，以縮短機(jī)時(shí)，使得CART3D的計(jì)算效率滿足優(yōu)化設(shè)計(jì)多輪計(jì)算的要求。在優(yōu)化過程中CART3D實(shí)際網(wǎng)格量約在20萬左右，在PC機(jī)上運(yùn)行一次的時(shí)間約為70～90s，滿足優(yōu)化設(shè)計(jì)要求。表2中計(jì)算條件與表1相同。另外，其網(wǎng)格生成高度自動(dòng)化，僅需數(shù)秒時(shí)間，效率很高。

表2 不同網(wǎng)格量的CART3D運(yùn)行時(shí)間與精度對比表Table 2 Comparison of run time and precision vs.grid quantity of CART3D

2 優(yōu)化設(shè)計(jì)算例驗(yàn)證

2.1 計(jì)算模型

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

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

2.2 計(jì)算目標(biāo)及約束

為確定CTV的最優(yōu)氣動(dòng)布局，這里采用以下簡化模型。在底部直徑D和總體積V不變的情況下，滿足配平升阻比CL／CD不小于0.8的約束，設(shè)計(jì)目標(biāo)為最大的體積利用率RV。設(shè)計(jì)變量為球頭半徑RN、前錐長LF、前錐半錐角θF。可表示為：

設(shè)計(jì)目標(biāo)： minU＝1／RV

約束條件：D＝3.5m

V＝3m3

CL／CD＞0.8

設(shè)計(jì)變量： 0.5m＜RN＜1.0m，2.0m＜LF＜3.0m，10.0°＜θF＜30.0°

其中，體積利用率RV＝6／S3／2，V為體積，S為表面積。設(shè)計(jì)點(diǎn)來流馬赫數(shù)M＝6，配平攻角αT＝20°。

2.3 算法選擇

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

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

2.4 結(jié)果分析

所選用的不同算法所得到的優(yōu)化解基本達(dá)到了優(yōu)化目的，即在滿足約束的條件下，優(yōu)化目標(biāo)體積利用率RV達(dá)到最大值，并且除算法2因?yàn)樵诰植孔顑?yōu)點(diǎn)中止優(yōu)化外，其余算法得到的結(jié)果，都略優(yōu)于文獻(xiàn)數(shù)據(jù)。對比結(jié)果如表3所示，表中長度單位為米，角度單位為度。

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

表3 不同優(yōu)化算法的優(yōu)化結(jié)果對比Table 3 Comparison between different optimization techniques

表4 算法6得到最優(yōu)解的氣動(dòng)力系數(shù)Table 4 Aerodynamic coefficient of the optimization result for algorithms 6

通過以上分析，可以確定適合于氣動(dòng)布局外形優(yōu)化的優(yōu)化算法。通過多島遺傳算法MIGA全局優(yōu)化，得到全局最優(yōu)解后再通過數(shù)值優(yōu)化算法進(jìn)行優(yōu)化計(jì)算，得到局部更優(yōu)解，而算法MIGA的參數(shù)設(shè)置將會(huì)影響到結(jié)果的優(yōu)化程度。這樣的混合優(yōu)化算法在避免陷入局部最優(yōu)解的同時(shí)，數(shù)值優(yōu)化算法的優(yōu)化計(jì)算對全局優(yōu)化算法得到的最優(yōu)解進(jìn)行“微調(diào)”，得到“更優(yōu)解”。

3 類神舟飛船返回艙優(yōu)化設(shè)計(jì)

3.1 單目標(biāo)優(yōu)化設(shè)計(jì)

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

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

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

目標(biāo)： maxCL／CD

約束：V≥16m3

Rv≥0.85

變量： 2.5m≤dm≤4.0m

選用MIGA和MMFD相結(jié)合的優(yōu)化算法，其中的多島遺傳算法設(shè)置島數(shù)為10，每島10個(gè)個(gè)體，進(jìn)化12代，染色體長度20，交叉概率0.9，變異概率0.01，遷徙率0.4，遷徙間隙3代。

3.2 雙目標(biāo)優(yōu)化設(shè)計(jì)

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

表5 單目標(biāo)優(yōu)化設(shè)計(jì)結(jié)果Table 5 Optimization result for single objective

目標(biāo)： maxCL／CD，minqw

約束：V≥16m3

Rv≥0.85

變量： 2.5m≤dm≤4.0m

選用鄰域培植遺傳算法NCGA，種群大小20，進(jìn)化150代，最后輸出pareto解集。兩個(gè)目標(biāo)的比例因子分別取為1.0和1000.0，權(quán)系數(shù)均為1.0，則目標(biāo)函數(shù)objective＝qw／1000.0－CL／CD／1.0。相應(yīng)優(yōu)化結(jié)果及優(yōu)化外形如表6、圖6所示。

表6 雙目標(biāo)優(yōu)化設(shè)計(jì)結(jié)果Table 6 Optimization result for double objectives

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

4 小 結(jié)

（1）形成了基于iSIGHT軟件平臺的集成優(yōu)化方法，用于返回艙優(yōu)化設(shè)計(jì)能夠得到理想結(jié)果。其優(yōu)化精度高、速度快，有較好的工程應(yīng)用前景。

（2）所采用的流場解算方法適用于多種復(fù)雜系統(tǒng)的氣動(dòng)力計(jì)算，適用性強(qiáng)，計(jì)算精度和效率較高。

（3）選擇MIGA與數(shù)值優(yōu)化算法相結(jié)合的混合遺傳算法，可以滿足單目標(biāo)多變量優(yōu)化設(shè)計(jì)。NCGA同樣可以滿足多目標(biāo)優(yōu)化設(shè)計(jì)，影響優(yōu)化結(jié)果的主要因素是遺傳算法選擇、交叉及變異等相關(guān)參數(shù)的選擇。下一步的研究工作將集中在優(yōu)化算法的選擇應(yīng)用及相應(yīng)優(yōu)化參數(shù)對優(yōu)化結(jié)果的影響方面。

參考文獻(xiàn):

［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］唐偉，桂業(yè)偉，張勇，等.未來升力再入飛船返回艙的氣動(dòng)布局設(shè)計(jì)［C］.中國宇航學(xué)會(huì)首屆學(xué)術(shù)年會(huì).北海，2005.

［6］唐偉，高曉成，李為吉，等.雙橢圓截面再入飛行器的氣動(dòng)計(jì)算及布局優(yōu)化設(shè)計(jì)［J］.空氣動(dòng)力學(xué)學(xué)報(bào)，2004，22（2）：171－174.

［7］唐偉，桂業(yè)偉，方方.新型升力再入飛船返回艙氣動(dòng)外形選型研究［J］.宇航學(xué)報(bào)，2008，29（1）：84－88.

［8］車競，唐碩，何開鋒，等.高超聲速飛行器氣動(dòng)布局總體性能優(yōu)化設(shè)計(jì)研究［J］.空氣動(dòng)力學(xué)學(xué)報(bào)，2009，27（2）：214－219.

［9］夏露，高正紅，李天.飛行器外形多目標(biāo)多學(xué)科綜合優(yōu)化設(shè)計(jì)方法研究［J］.空氣動(dòng)力學(xué)學(xué)報(bào)，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］袁先旭.非定常流動(dòng)數(shù)值模擬及飛行器動(dòng)態(tài)特性分析研究［J］.［博士學(xué)位論文］.中國空氣動(dòng)力研究與發(fā)展中心，2002.

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

［13］安偉剛.多目標(biāo)優(yōu)化方法研究及其工程應(yīng)用［D］.［博士學(xué)位論文］.西北工業(yè)大學(xué)，2005.

［14］唐偉，張勇，李為吉，等.歐洲返回艙CTV氣動(dòng)特性綜述［J］.空氣動(dòng)力學(xué)學(xué)報(bào)，2003，21（3）：362－366.