梁欣欣，李世鵬，陳 陽(yáng)，張 艷，傅 瑜

（北京宇航系統(tǒng)工程研究所，北京100076）

在滿(mǎn)足一定的固體發(fā)動(dòng)機(jī)總體性能的設(shè)計(jì)約束下，內(nèi)彈道決定了發(fā)動(dòng)機(jī)推力大小的變化趨勢(shì)，飛行程序角則決定了推力方向，將二者聯(lián)合尋優(yōu)，可使固體運(yùn)載器的總體方案和飛行程序更加優(yōu)化［1－4］。在開(kāi)展固體運(yùn)載器總體設(shè)計(jì)的初步論證過(guò)程中，普遍將發(fā)動(dòng)機(jī)內(nèi)彈道以平均推力曲線(xiàn)形式來(lái)開(kāi)展發(fā)動(dòng)機(jī)的設(shè)計(jì)與總體方案論證。事實(shí)上，平均推力曲線(xiàn)形式的內(nèi)彈道設(shè)計(jì)對(duì)于實(shí)現(xiàn)射程指標(biāo)以及大氣飛行段姿控穩(wěn)定、載荷條件等并非最優(yōu)。本文旨在探尋固體運(yùn)載器具有較好性能時(shí)的發(fā)動(dòng)機(jī)內(nèi)彈道曲線(xiàn)形式，在實(shí)際工程研制過(guò)程中可通過(guò)發(fā)動(dòng)機(jī)藥柱設(shè)計(jì)來(lái)改變藥柱燃燒過(guò)程中的燃面，獲得有利于運(yùn)載器總體性能的曲線(xiàn)形式。

本文以射程最大為目標(biāo)，用離散化方法建立了內(nèi)彈道優(yōu)化設(shè)計(jì)模型，并開(kāi)展了內(nèi)彈道優(yōu)化仿真，在滿(mǎn)足飛行過(guò)程中的姿控和載荷約束的同時(shí)，提高固體運(yùn)載器的射程。

1 數(shù)學(xué)模型

1.1 飛行程序設(shè)計(jì)

固體運(yùn)載器飛行程序采用目前通常使用的工程設(shè)計(jì)方法［5］。在主動(dòng)段，運(yùn)載器垂直起飛一段時(shí)間后，利用攻角轉(zhuǎn)彎段調(diào)整推力的方向，此后進(jìn)入攻角為0的重力轉(zhuǎn)彎段，在減小氣動(dòng)損失的同時(shí)，達(dá)到彈道轉(zhuǎn)彎的目的，然后進(jìn)入定軸俯仰程序角飛行段。被動(dòng)段為無(wú)動(dòng)力橢圓軌道飛行段。固體運(yùn)載器主動(dòng)段俯仰飛行程序角曲線(xiàn)如圖1所示。

圖1 固體運(yùn)載器主動(dòng)段俯仰飛行程序角曲線(xiàn)

圖1 中，t1可由運(yùn)載器推重比確定，t2在運(yùn)載器飛行亞音速段結(jié)束，重力轉(zhuǎn)彎段可在大氣稀薄時(shí)的某一高度結(jié)束。通過(guò)以上程序角設(shè)計(jì)原則，推力方向的變化僅與攻角轉(zhuǎn)彎段的最大負(fù)攻角αM有關(guān)。

1.2 飛行動(dòng)力學(xué)建模

在方案論證階段，可將固體運(yùn)載器的運(yùn)動(dòng)視為質(zhì)點(diǎn)運(yùn)動(dòng)，在發(fā)射坐標(biāo)系內(nèi)建立質(zhì)點(diǎn)彈道模型［6－7］。固體運(yùn)載器運(yùn)動(dòng)學(xué)方程為

式中：x，y，z為運(yùn)載器位置在發(fā)射系的坐標(biāo)；vx，vy，vz為速度在發(fā)射系的分量。

動(dòng)力學(xué)方程為

式中：GB為彈體系到發(fā)射系的轉(zhuǎn)換矩陣；GV為速度系到發(fā)射系的轉(zhuǎn)換矩陣；Fx，F(xiàn)y，F(xiàn)z分別為氣動(dòng)阻力、升力和側(cè)力；gx、gy、gz為引力在發(fā)射系的分量；Fe為發(fā)動(dòng)機(jī)的有效推力，按下式計(jì)算：

式中：F0（t）為發(fā)動(dòng)機(jī)地面推力，Se為發(fā)動(dòng)機(jī)噴口面積，F(xiàn)v（t）為發(fā)動(dòng)機(jī)的真空推力，p為大氣壓強(qiáng)，p0為地面大氣壓強(qiáng)。

氣動(dòng)力按下式計(jì)算：

式中：Cx為阻力系數(shù)，Cαy為升力系數(shù)對(duì)攻角α的導(dǎo)數(shù)，Cβz為側(cè)力系數(shù)對(duì)側(cè)滑角β的導(dǎo)數(shù)，SM為固體運(yùn)載器特征面積，q為動(dòng)壓。

地球采用靜止圓球假設(shè)，因此引力可按下式計(jì)算：

式中：g沿地心矢徑方向，可根據(jù)角度關(guān)系獲得發(fā)射系分量gx，gy，gz；μ為地球引力系數(shù)；r為地心矢徑。

質(zhì)量按下式計(jì)算：

式中：m為運(yùn)載器實(shí)時(shí)質(zhì)量，m0為各級(jí)起飛質(zhì)量，為各級(jí)發(fā)動(dòng)機(jī)的質(zhì)量流量。

2 優(yōu)化模型及求解方法

2.1 優(yōu)化模型

本文建立推力曲線(xiàn)的假設(shè)條件：①發(fā)動(dòng)機(jī)裝藥量、總沖、工作時(shí)間保持不變；②內(nèi)彈道推力與壓強(qiáng)、流量成正比；③不考慮發(fā)動(dòng)機(jī)上升沿與下降沿的曲線(xiàn)變化；④不考慮內(nèi)彈道曲線(xiàn)的變化導(dǎo)致的發(fā)動(dòng)機(jī)總重的變化。

將內(nèi)彈道推力曲線(xiàn)離散化為多個(gè)結(jié)點(diǎn)，將每個(gè)結(jié)點(diǎn)上的推力值作為優(yōu)化變量，然后采用擬合方法獲得整個(gè)內(nèi)彈道推力曲線(xiàn)。

離散及擬合方法在軌跡優(yōu)化的各類(lèi)算法中應(yīng)用較為廣泛，本文采用等時(shí)間間隔離散、線(xiàn)性擬合方法。將發(fā)動(dòng)機(jī)工作時(shí)間為T(mén)的內(nèi)彈道推力曲線(xiàn)F（t）離散化為等時(shí)間間隔的N段，N＋1個(gè)時(shí)間節(jié)點(diǎn)上的推力值分別為F0，F(xiàn)1，…，F(xiàn)N，n＝1，2，…，N，則曲線(xiàn)F（t）可利用時(shí)間節(jié)點(diǎn)上的推力值線(xiàn)性擬合獲得，即：

推力結(jié)點(diǎn)值F0，F(xiàn)1，…，F(xiàn)N和最大負(fù)攻角αM構(gòu)成了優(yōu)化問(wèn)題的優(yōu)化變量。發(fā)動(dòng)機(jī)總沖為定值，且推力結(jié)點(diǎn)值受發(fā)動(dòng)機(jī)承壓影響需在一定范圍內(nèi)變化，即構(gòu)成了優(yōu)化問(wèn)題的約束條件。綜上，可建立以最大射程為優(yōu)化目標(biāo)的內(nèi)彈道優(yōu)化模型，如下所示：

式中：J（X）為目標(biāo)函數(shù)，X為優(yōu)化變量，L為射程，I為發(fā)動(dòng)機(jī)總沖，F(xiàn)為平均推力，δ為平均推力上下浮動(dòng)的百分比，可根據(jù)發(fā)動(dòng)機(jī)可承受壓強(qiáng)確定。流量的計(jì)算可根據(jù)與推力成正比的關(guān)系計(jì)算得到。

為了將上述帶約束的優(yōu)化問(wèn)題轉(zhuǎn)化為無(wú)約束的優(yōu)化問(wèn)題，將待優(yōu)化的各節(jié)點(diǎn)推力值進(jìn)行等比調(diào)節(jié)和值域調(diào)控的迭代處理，使得約束自動(dòng)滿(mǎn)足。無(wú)約束的優(yōu)化模型為

將帶約束的優(yōu)化問(wèn)題（8）轉(zhuǎn)換為無(wú)約束問(wèn)題（9）的流程，如圖2所示。

圖2 帶約束的優(yōu)化問(wèn)題轉(zhuǎn)換為無(wú)約束問(wèn)題的流程圖

2.2 優(yōu)化流程

針對(duì)優(yōu)化模型，可將攻角轉(zhuǎn)彎段的最大負(fù)攻角利用一維搜索獲得優(yōu)化值，其余優(yōu)化變量利用成熟的優(yōu)化算法進(jìn)行優(yōu)化。針對(duì)優(yōu)化模型（8）的仿真優(yōu)化流程如圖3所示。

圖3 離散法優(yōu)化流程

2.3 優(yōu)化算法

本文采用自適應(yīng)差分進(jìn)化法［8－9］進(jìn)行優(yōu)化求解。1）基本思想。

自適應(yīng)差分進(jìn)化法的基本思想為：針對(duì)種群中個(gè)體適應(yīng)度的分散程度進(jìn)行度量，并利用個(gè)體的適應(yīng)度及其分散程度進(jìn)行交叉和變異操作，用種群熵表征種群搜索過(guò)程中的收斂程度，并根據(jù)種群熵的變化自適應(yīng)地減小種群的搜索范圍，即在搜索初期快速減小種群的搜索范圍，提高早期收斂的快速性；隨著搜索的進(jìn)行，逐步降低種群搜索范圍減小的速度，以保證全局收斂性能，避免局部收斂。

2）基本操作。

設(shè)種群規(guī)模為m，設(shè)Fmin和Fmax是在搜索過(guò)程中第k代種群的最小和最大適應(yīng)度值，將適應(yīng)度值劃分為等距的m個(gè)Aj（j＝1，2，…，m）區(qū)域：

設(shè)種群中有mj個(gè)個(gè)體的適應(yīng)度值屬于Aj區(qū)域，定義個(gè)體出現(xiàn)在Aj區(qū)域的概率為Pj＝mj／m。由此，定義第k代種群的熵為設(shè)krc、S分別為初始種群或改變范圍后種群的代數(shù)和熵，kmax為最大的設(shè)定代數(shù)，定義SR為是否改變種群范圍的控制參數(shù)：SR＝0.98exp（－krc／kmax）。

設(shè)優(yōu)化變量xj∈ ［，］，j＝1，2，…，m，和分別為第k代種群中適應(yīng)度值前60%的個(gè)體的第j個(gè)變量的最大值和最小值。當(dāng)Sk／S＜SR時(shí)，改變種群搜索的范圍：

式中：w為松弛因子。在搜索初期，SR較大，保證算法快速縮小種群的搜索范圍，提高搜索效率。隨著搜索的進(jìn)行，SR逐漸變小，種群搜索范圍縮小的速度減緩，以保證算法的全局收斂性。

3）算法流程。

自適應(yīng)差分進(jìn)化法的流程如圖4所示。

圖4 自適應(yīng)差分進(jìn)化算法流程圖

3 優(yōu)化結(jié)果及分析

針對(duì)典型三級(jí)固體運(yùn)載器大氣飛行段開(kāi)展分析。對(duì)于載荷條件，由于該飛行器一級(jí)飛行段不是最?lèi)毫虞d荷環(huán)境，其最大動(dòng)壓提高30%依然滿(mǎn)足載荷條件；對(duì)于姿控條件，一二級(jí)分離的高度盡可能高，一二級(jí)分離動(dòng)壓盡可能小，則對(duì)姿控條件越有利。

利用自適應(yīng)差分進(jìn)化法開(kāi)展離散模型的優(yōu)化。將一級(jí)發(fā)動(dòng)機(jī)工作時(shí)間等間隔分為24段，以25個(gè)節(jié)點(diǎn)值作為優(yōu)化變量，最優(yōu)轉(zhuǎn)彎攻角αM在25個(gè)節(jié)點(diǎn)值確定后由一維搜索算法獲得，δ＝0.3，即F（t）內(nèi)彈道曲線(xiàn)各節(jié)點(diǎn)值上、下限分別設(shè)定為平均推力值的130%和70%。種群規(guī)模為40，迭代收斂最大代數(shù)為200代。

優(yōu)化迭代收斂曲線(xiàn)如圖5所示，得到的優(yōu)化內(nèi)彈道曲線(xiàn)如圖6所示，優(yōu)化結(jié)果與平均推力對(duì)比如表1所示。圖5中，Lmax為每代最大射程值，表1中以平均模型為基準(zhǔn)給出了各量?jī)?yōu)化后的變化百分比，L為射程，qmax為最大動(dòng)壓，Hs為一二級(jí)分離高度，qs為一二級(jí)分離動(dòng)壓。

圖5 優(yōu)化收斂曲線(xiàn)

圖6 優(yōu)化內(nèi)彈道曲線(xiàn)

表1 優(yōu)化結(jié)果與平均推力對(duì)比

由圖表可知，“前高后低”的雙推力內(nèi)彈道曲線(xiàn)對(duì)于提高射程來(lái)說(shuō)具有積極意義，可提高射程約3.1%；同時(shí)，由于跨音速段推力大，導(dǎo)致最大動(dòng)壓提高27.1%，小于可承受最大動(dòng)壓；一二級(jí)分離高度提高為115.6%，分離動(dòng)壓減小為37.8%，這對(duì)于姿控分離穩(wěn)定十分有利。因此，在固體運(yùn)載器方案論證時(shí)，可針對(duì)“前高后低”的雙推力內(nèi)彈道形式對(duì)大氣飛行段工作的發(fā)動(dòng)機(jī)開(kāi)展研究論證工作，保證在符合發(fā)動(dòng)機(jī)性能約束和運(yùn)載器飛行載荷環(huán)境約束的情況下，提高射程能力，改善一二級(jí)分離環(huán)境。

4 結(jié)論

本文以固體運(yùn)載器為研究對(duì)象，針對(duì)大氣飛行段內(nèi)彈道優(yōu)化設(shè)計(jì)建模，并以自適應(yīng)差分進(jìn)化法進(jìn)行仿真驗(yàn)證。通過(guò)對(duì)仿真結(jié)果的分析，有以下幾點(diǎn)結(jié)論：①將內(nèi)彈道及外彈道聯(lián)合優(yōu)化，可進(jìn)一步改善固體運(yùn)載器的整體性能和飛行環(huán)境；②內(nèi)彈道曲線(xiàn)優(yōu)化形式并非傳統(tǒng)認(rèn)為的平均推力形式；③“前高后低”的雙推力的內(nèi)彈道曲線(xiàn)形式對(duì)于提高射程和改善飛行環(huán)境具有積極意義。

研究成果可為固體運(yùn)載體總體方案論證提供借鑒意義。

［1］JAROSLAW S S.Sensitivity analysis and multidisciplinary optimization for aircraft design：recent advances and results［J］.Journal of Aircraft，1990，27（12）：993－1001.

［2］BROWN N F，OLDS J R.Evaluation of multidisciplinary optimization techniques applied to a reusable launch vehicle［J］.Journal of Spacecraft and Rockets，2006，43（6）：1 289－1 300.

［3］BRAUN R D，POWELL R W.Comparison of two multidisciplinary optimization strategies for launch－vehicle design［J］.Journal of Spacecraft and Rockets，1995，32（3）：404－410.

［4］楊軍，陳汝訓(xùn)，趙錫良.戰(zhàn)術(shù)火箭／固體火箭發(fā)動(dòng)機(jī)一體化優(yōu)化設(shè)計(jì)［J］.宇航學(xué)報(bào)，1999，20（1）：21－27.YANG Jun，CHEN Ru－xun，ZHAO Xi－liang.Integral optimization design for tactical rocket and it’s SRM by utilizing the synthesis of internal ballistics and trajectory［J］.Journal of Astronautics，1999，20（1）：21－27.（in Chinese）

［5］龍樂(lè)豪.總體設(shè)計(jì)［M］.北京：中國(guó)宇航出版社，1989.LONG Le－h(huán)ao.Overall design［M］.Beijing：China Astronautic Publishing House，1989.（in Chinese）

［6］賈沛然，陳克俊，何力.遠(yuǎn)程火箭彈道學(xué)［M］.長(zhǎng)沙：國(guó)防科技大學(xué)出版社，1994.JIA Pei－ran，CHEN Ke－jun，HE Li.The ballistics of long－range rocket［M］.Changsha：National University of Defense Technology Press，1994.（in Chinese）

［7］洪蓓，梁欣欣，辛萬(wàn)青.固體運(yùn)載火箭多約束彈道優(yōu)化［J］.導(dǎo)彈與航天運(yùn)載技術(shù)，2012，（3）：1－5.HONG Bei，LIANG Xin－xin，XIN Wan－qing.Multi－constraint trajectory optimization for solid launch vehicle［J］.Missiles and Space Vehicles，2012，（3）：1－5.（in Chinese）

［8］PRICE K.Differential evolution：a fast and simple numerical optimizer［C］／／Biennial Conf of the North American Fuzzy Information Processing Sociey.New York：IEEE 1996：524－527.

［9］PRICE K.Differential evolution vs.the functions of the 2nd ICEO［C］／／IEEE Int Conf on Evolutionary Computation.Indian Upolis：IEEE，1997：153－157.