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

蔡致鵬,王新閣,張宴嘉

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

1 引言

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

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

2 吻切錐乘波體模型建立

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

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

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

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

(1)

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

圖2 乘波體底面截面圖

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

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

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

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

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

(2)

3.2 萬(wàn)有引力搜索算法

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

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

(3)

(4)

(5)

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

best(t)=maxfitj(t)

(6)

worst(t)=minfitj(t)

(7)

t時(shí)刻粒子i受到粒子j的引力為：

(8)

(9)

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

(10)

(11)

(12)

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

(13)

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

(14)

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

圖3 萬(wàn)有引力搜索算法流程框圖

3.3 升阻比和容積率估算

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

(15)

其中:V為乘波體飛行器的體積；Sw為乘波體飛行器的濕潤(rùn)面積。

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

0.037 6c-0.039 9L1

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

容積率η和升阻比L/D是一對(duì)存在矛盾關(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所示。選取點(diǎn)A(0.140 001，4.740 26)和點(diǎn)B(0.128 073，5.514 019)建立優(yōu)化模型。點(diǎn)A對(duì)應(yīng)的c和L1為(2.1，3.7)，點(diǎn)B對(duì)應(yīng)的c和L1為(1.8，3.3)。建立相應(yīng)的模型A和模型B，如圖5、圖6所示。

表1 正交實(shí)驗(yàn)計(jì)算結(jié)果

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

圖5 模型A示意圖

圖6 模型B示意圖

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

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

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

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

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

圖7 模型B的壓力云圖

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

圖8 乘波體截面壓力云圖

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

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

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

4 結(jié)論

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

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

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