基于引力搜索算法的兩級串聯(lián)吻切錐乘波體設(shè)計

蔡致鵬,王新閣,張宴嘉

(空軍航空大學(xué)， 長春 130000)

1 引言

高超聲速武器是指飛行馬赫數(shù)達到5馬赫以上，能夠在大氣層或臨近空間飛行的武器[1]。因其具有高空高速、不易攔截等優(yōu)勢而深受追捧。未來基于先進的超音速轟炸機平臺，高超聲速巡航導(dǎo)彈勢必會得到更加廣泛的應(yīng)用。高超聲速巡航導(dǎo)彈相對于亞音速巡航導(dǎo)彈有反應(yīng)速度快、突防能力強、毀傷效果大等優(yōu)點[1,2]，是目前各國都在研究的主流武器之一。美在高超聲速武器研究中起步早，美國開發(fā)了多款高超聲速武器。X-51A驗證機進行了多次實驗，驗證機由B52轟炸機攜帶到1.5萬m高空投放，之后火箭助推器把驗證機加速到4.8Ma后脫落分離，機身內(nèi)的超燃沖壓發(fā)動機啟動，速度保持在5.1Ma，飛行時長約為4 min。先進高超聲速武器(advanced hypersonic weapons,AHW)計劃的目標(biāo)是發(fā)展一種射程6 000 km，飛行速度6馬赫，精度在10 m以內(nèi)的導(dǎo)彈，到2017年底前AWH進行過3次試飛，2次取得成功[3]。自美國退出《反彈道導(dǎo)彈條約》以來，俄羅斯為突破美國導(dǎo)彈防御系統(tǒng)，一直致力于研發(fā)高超聲速武器。匕首高超音速導(dǎo)彈由陸基伊斯坎德爾-M近程彈道導(dǎo)彈改裝而來，由米格-31K截擊機或圖-22M轟炸機攜帶投放，飛行速度10Ma，圓概率誤差7 m。俄羅斯的“鋯石”高超聲速巡航導(dǎo)彈經(jīng)過多次成功試射后已經(jīng)服役[3-4]。

常規(guī)構(gòu)型飛行器高超聲速飛行時，下表面高壓氣體會上溢至上表面，上表面壓力增大，導(dǎo)致升阻比降低。Kuchemann[5]根據(jù)風(fēng)洞試驗和試飛結(jié)果，提出了高超聲速“升阻比屏障”，傳統(tǒng)氣動構(gòu)型能達到的高超聲速升阻比極限為4。乘波體構(gòu)型是能突破“升阻比屏障”的飛行器外形之一，其設(shè)計思想是采用已知基準流場反設(shè)計飛行器外形。自1959年提出以來，發(fā)展出楔形乘波體、錐導(dǎo)乘波體、吻切錐乘波體和楔-錐乘波體4種外形[6]。面對日益復(fù)雜的周邊環(huán)境和世界軍事大國的領(lǐng)先地位，加緊對高超聲速巡航導(dǎo)彈的研究具有重要意義。本文中基于低-高馬赫數(shù)兩級串聯(lián)吻切錐乘波體對高超聲速巡航導(dǎo)彈氣動外形進行研究。

2 吻切錐乘波體模型建立

當(dāng)超聲速氣流繞流無限長零迎角圓錐時，且氣流馬赫數(shù)和圓錐半頂角在一定范圍內(nèi)，圓錐前會產(chǎn)生錐面激波。乘波體能夠?qū)⒓げê蟮母邏簹怏w限制在乘波體的下表面，阻止下表面氣體繞過前緣曲線與上表面形成連通，因而會得到更大的升力。

如圖1所示，將激波出口型線劃分為多個小段激波曲線，每小段曲線均可視為部分圓錐激波，產(chǎn)生該段激波的圓錐為吻切錐。預(yù)先選定的吻切錐乘波體流動捕捉曲線FCC(flow capture curve)通過自由流線法獲得流動捕捉管FCT(flow capture ture)。FCT與每個吻切錐所產(chǎn)生的交線即為吻切錐乘波體的前緣曲線。前緣曲線與FCT合圍部分為乘波體上表面。將前緣曲線按照激波出口型線劃分成對應(yīng)的幾段，每段前緣曲線都進行離散化，分別進行流線追蹤，最后將流線經(jīng)過曲面擬合生成乘波體下表面。兩級串聯(lián)吻切錐乘波體由一級乘波體、銜接段、二級乘波體組成。一、二級長度h相同，采用橋接曲面連接，銜接段長度l為全長的十分之一。

圖1 兩級串聯(lián)吻切錐乘波體截面示意圖

當(dāng)來流馬赫數(shù)M，內(nèi)圓錐半錐角δ，和流場長度L0固定時，乘波體的FCC曲線就成了改變乘波體外形的唯一變量。如圖2所示，本文中設(shè)定吻切錐乘波體FCC曲線為：

(1)

式中：a、b、c均為實數(shù)，保證曲線的指數(shù)段和水平段二階導(dǎo)數(shù)連續(xù)。式中的參數(shù)a、b、c、L1并不相互獨立。改變兩級串聯(lián)吻切錐乘波體外形僅有c、L1兩個變量。為滿足寬速域要求，兼顧投放加速起點，一級乘波體設(shè)計馬赫數(shù)為4，二級乘波體設(shè)計馬赫數(shù)為8。

圖2 乘波體底面截面圖

3 基于引力搜索算法的兩級串聯(lián)吻切錐乘波體設(shè)計

3.1 多目標(biāo)優(yōu)化

多目標(biāo)優(yōu)化問題(multi-objective optimization，MPO)一般由N個決策變量參數(shù)，M個目標(biāo)函數(shù)和J個約束條件組成，目標(biāo)函數(shù)、約束條件與決策變量之間呈函數(shù)關(guān)系。對于最大化多目標(biāo)優(yōu)化問題的描述有：

Maxf(x)=[fi(x),i=1,2,…,M]

gj(x)≤0,j=1,2,…,J

(2)

3.2 萬有引力搜索算法

引力搜索算法(gravitational search algorithm,GSA)是Rasshedi等人受牛頓的萬有引力定律啟發(fā)，于2009年提出的一種算法，引力算法的個體位置代表著問題的解，個體質(zhì)量是對結(jié)果優(yōu)劣的評價[9]。

設(shè)在D維空間中有N個粒子，則第i個粒子的位置為：

(3)

(4)

(5)

其中：fiti(t)為t時刻粒子i的適應(yīng)度；worst(t)表示最差適應(yīng)度；best(t)表示最佳適應(yīng)度。在求解最大化問題中：

best(t)=maxfitj(t)

(6)

worst(t)=minfitj(t)

(7)

t時刻粒子i受到粒子j的引力為：

(8)

(9)

其中：G(t)為引力常數(shù)；為了防止分母為零，添加一個大于零的小常數(shù)ε；Rij(t)為粒子之間的歐氏距離。粒子i受到的引力為粒子群內(nèi)除本身外其他粒子吸引力的合力。通過采用各粒子施加引力的隨機加權(quán)和代表粒子所受的合力，從而加強隨機性。

(10)

(11)

(12)

但在尋找最優(yōu)解的過程中，結(jié)果可能很快收斂在次優(yōu)解附近，并反復(fù)跳躍，難以得到全局最優(yōu)解。為此，任意粒子i上所受合外力不再是種群中除粒子i外全部粒子所提供的引力和，算法將選擇部分粒子發(fā)揮作用，該部分粒子的集合記為Kbest。調(diào)整后粒子i所受合力為：

(13)

Kbest內(nèi)的元素為種群中除粒子i外適應(yīng)度值排在前Kbest個的全部粒子，且集合內(nèi)的元素會隨著迭代步數(shù)的增加逐漸減少。從初始化的種群大小為N，經(jīng)過多次迭代，Kbest內(nèi)部的元素可能會減少到1個。為了算法能夠更加準確高效，GSA算法要求引力常數(shù)G隨時間t變化，調(diào)整后的萬有引力常數(shù)G(t)為：

(14)

其中：G0為萬有引力常數(shù)的初始值；β為引力常數(shù)縮減系數(shù)。整體算法流程如圖3所示[10]。

圖3 萬有引力搜索算法流程框圖

3.3 升阻比和容積率估算

高超聲速飛行器的設(shè)計應(yīng)同時滿足高升阻比、高升力系數(shù)和高容積率需求[11-12]。本文中選擇升阻比和容積率作為優(yōu)化目標(biāo)，容積率η的計算公式為[13]：

(15)

其中:V為乘波體飛行器的體積；Sw為乘波體飛行器的濕潤面積。

由于種群數(shù)量N較大，為精簡計算量，提高效率。采用正交實驗和半經(jīng)驗公式推導(dǎo)總結(jié)設(shè)計參數(shù)與性能之間關(guān)系，并建立二元非線性回歸模型，使用對升阻比和容積率的估算結(jié)果來代替模擬結(jié)果。通過改變c和L1這2個設(shè)計參數(shù)可以改變兩級串聯(lián)吻切錐乘波體外形，每個參數(shù)選取4個水平，故采用L8(42)正交表，共需要8次計算。在Fluent中進行8次計算得到結(jié)果匯總為表1。根據(jù)結(jié)果建立多元非線性回歸模型。參考升阻比估算和容積率的計算公式[14]，可以得到升阻比L/D和容積率η都與L1和c成二階關(guān)系。最后按照正交表的試驗結(jié)果，最終建立函數(shù)模型：

0.037 6c-0.039 9L1

3.4 升阻比和容積率優(yōu)化

容積率η和升阻比L/D是一對存在矛盾關(guān)系的參數(shù)。取維度d=2，粒子群數(shù)N=200，Pareto解集的最大規(guī)模為100，c的邊界條件為[1,2.1 838]，L1邊界條件為[2,4.867 6]，容積率η大于0.128，升阻比L/D大于3.2，設(shè)定迭代次數(shù)為200。得到容積率η與升阻比L/D的Pareto解集分布如圖4所示。選取點A(0.140 001，4.740 26)和點B(0.128 073，5.514 019)建立優(yōu)化模型。點A對應(yīng)的c和L1為(2.1，3.7)，點B對應(yīng)的c和L1為(1.8，3.3)。建立相應(yīng)的模型A和模型B，如圖5、圖6所示。

表1 正交實驗計算結(jié)果

圖4 多目標(biāo)優(yōu)化的Pareto解集曲線

圖5 模型A示意圖

圖6 模型B示意圖

使用Fluent對模型進行升阻比和容積率計算，得到如表2所示的結(jié)果。經(jīng)過計算，仿真計算結(jié)果與估算結(jié)果差值較小，證明估算結(jié)果具有可靠性。

表2 模型A、B的容積率、升阻比Table 2 Volume ratio and lift-drag ratio of model A and B

由此能夠提出一種反設(shè)計思路，可以在設(shè)計導(dǎo)彈之前優(yōu)先依據(jù)任務(wù)的需要，明確攜帶燃料和彈藥部質(zhì)量，由此選擇合適容積率，然后通過圖4中的解集，找到最大的升阻比。最后再依據(jù)選擇的升阻比與容積率約束兩級串聯(lián)吻切錐乘波體的設(shè)計參數(shù)。

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

選取模型B進行氣動分析。根據(jù)圖7(a)可以看出，在速度為6Ma、迎角為4°的條件下，氣流經(jīng)過模型B能夠生成激波且提供較大的升力。在考慮底阻的情況下，尾部離體后的區(qū)域會形成一段低壓區(qū)。可以通過在尾部安裝動力裝置的方式，利用尾部噴出的高壓氣體來抵消乘波體前后壓力差，從而能夠進一步增大升阻比。圖7(b)中，可以在多個截面上觀察到，在乘波體下表面會形成一個高壓封閉區(qū)，雖然邊緣位置有少量與上表面形成連通區(qū)，但在整體上的壓力分布上保持一個較好水平。

圖7 模型B的壓力云圖

如圖8(a)、(b)、(c)和(d)分別為在來流速度為4Ma時模型B的一級乘波體底面、二級乘波體底面的壓力云圖和來流速度為8Ma時模型B的銜接段、二級乘波體底面的壓力云圖。

圖8 乘波體截面壓力云圖

從4張圖都可以看出，壓力最大區(qū)域集中在下表面邊緣位置。因而在邊緣位置的封閉性就直接影響了升力的大小。對比圖8(a)圖8(b)可以發(fā)現(xiàn)，一級乘波體的下表面封閉性優(yōu)于二級乘波體。一級乘波體的設(shè)計馬赫數(shù)為4，二級乘波體的設(shè)計馬赫數(shù)為8，故在4Ma條件下，一級乘波體的壓力封閉性能更好。分析圖8(c)與圖8(d)，發(fā)現(xiàn)銜接段的封閉性要好于二級乘波體。將圖8(b)與圖8(d)進行對比，在8馬赫的條件下，壓力外泄的情況得到了明顯改善。上述結(jié)果都與預(yù)期設(shè)定吻合，證明達到了反設(shè)計的預(yù)期效果。

圖9為兩級串聯(lián)吻切錐乘波體在不同速度下的升阻比曲線。從圖9中看出，從4～8Ma部分升阻比隨著速度的增大會逐步增大，但是在馬赫數(shù)大于8部分，升阻比趨于穩(wěn)定且略有下降。根據(jù)馬赫數(shù)無關(guān)原理[15-16]，當(dāng)馬赫數(shù)增大到一定程度之后，飛行器表面的升力系數(shù)與阻力系數(shù)將不隨馬赫數(shù)發(fā)生變化，所以升阻比在8Ma后逐漸趨于穩(wěn)定。但是底部阻力會略有增大，從而導(dǎo)致升阻比略有下降。

圖9 模型B升阻比與馬赫數(shù)之間的關(guān)系曲線

4 結(jié)論

1) 采用正交實驗研究了設(shè)計參數(shù)與升阻比容積率的關(guān)系，得到升阻比L/D和容積率η都與L1和c成二階關(guān)系，并建立多元非線性回歸模型。分析發(fā)現(xiàn)容積率高的乘波體，其升阻比性能不佳，而升阻比高的乘波體，容積率又無法滿足要求，2個參數(shù)是矛盾關(guān)系。

2) 采用多目標(biāo)萬有引力搜索算法對兩級串聯(lián)吻切錐乘波體進行優(yōu)化，在滿足容積率需求的前提下得到最大升阻比的優(yōu)化模型，并通過CFD方法驗證，容積率和升阻比的誤差在1%以內(nèi)，證明了優(yōu)化方法的可靠性。

3) 對優(yōu)化后的乘波體模型仿真分析，一、二級乘波體在各自的設(shè)計馬赫數(shù)下都表現(xiàn)出較好的封閉性，在6Ma時整體的壓力分布較好，有較大升阻比，在馬赫數(shù)超過8后，升阻比略有下降后趨于穩(wěn)定。