梁欣欣，李世鵬，陳 陽，張 艷，傅 瑜

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

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

本文以射程最大為目標，用離散化方法建立了內彈道優(yōu)化設計模型，并開展了內彈道優(yōu)化仿真，在滿足飛行過程中的姿控和載荷約束的同時，提高固體運載器的射程。

1 數(shù)學模型

1.1 飛行程序設計

固體運載器飛行程序采用目前通常使用的工程設計方法［5］。在主動段，運載器垂直起飛一段時間后，利用攻角轉彎段調整推力的方向，此后進入攻角為0的重力轉彎段，在減小氣動損失的同時，達到彈道轉彎的目的，然后進入定軸俯仰程序角飛行段。被動段為無動力橢圓軌道飛行段。固體運載器主動段俯仰飛行程序角曲線如圖1所示。

圖1 固體運載器主動段俯仰飛行程序角曲線

圖1 中，t1可由運載器推重比確定，t2在運載器飛行亞音速段結束，重力轉彎段可在大氣稀薄時的某一高度結束。通過以上程序角設計原則，推力方向的變化僅與攻角轉彎段的最大負攻角αM有關。

1.2 飛行動力學建模

在方案論證階段，可將固體運載器的運動視為質點運動，在發(fā)射坐標系內建立質點彈道模型［6－7］。固體運載器運動學方程為

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

動力學方程為

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

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

氣動力按下式計算：

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

地球采用靜止圓球假設，因此引力可按下式計算：

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

質量按下式計算：

式中：m為運載器實時質量，m0為各級起飛質量，為各級發(fā)動機的質量流量。

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

2.1 優(yōu)化模型

本文建立推力曲線的假設條件：①發(fā)動機裝藥量、總沖、工作時間保持不變；②內彈道推力與壓強、流量成正比；③不考慮發(fā)動機上升沿與下降沿的曲線變化；④不考慮內彈道曲線的變化導致的發(fā)動機總重的變化。

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

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

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

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

為了將上述帶約束的優(yōu)化問題轉化為無約束的優(yōu)化問題，將待優(yōu)化的各節(jié)點推力值進行等比調節(jié)和值域調控的迭代處理，使得約束自動滿足。無約束的優(yōu)化模型為

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

圖2 帶約束的優(yōu)化問題轉換為無約束問題的流程圖

2.2 優(yōu)化流程

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

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

2.3 優(yōu)化算法

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

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

2）基本操作。

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

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

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

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

3）算法流程。

自適應差分進化法的流程如圖4所示。

圖4 自適應差分進化算法流程圖

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

針對典型三級固體運載器大氣飛行段開展分析。對于載荷條件，由于該飛行器一級飛行段不是最惡劣載荷環(huán)境，其最大動壓提高30%依然滿足載荷條件；對于姿控條件，一二級分離的高度盡可能高，一二級分離動壓盡可能小，則對姿控條件越有利。

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

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

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

圖6 優(yōu)化內彈道曲線

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

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

4 結論

本文以固體運載器為研究對象，針對大氣飛行段內彈道優(yōu)化設計建模，并以自適應差分進化法進行仿真驗證。通過對仿真結果的分析，有以下幾點結論：①將內彈道及外彈道聯(lián)合優(yōu)化，可進一步改善固體運載器的整體性能和飛行環(huán)境；②內彈道曲線優(yōu)化形式并非傳統(tǒng)認為的平均推力形式；③“前高后低”的雙推力的內彈道曲線形式對于提高射程和改善飛行環(huá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］楊軍，陳汝訓，趙錫良.戰(zhàn)術火箭／固體火箭發(fā)動機一體化優(yōu)化設計［J］.宇航學報，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］龍樂豪.總體設計［M］.北京：中國宇航出版社，1989.LONG Le－h(huán)ao.Overall design［M］.Beijing：China Astronautic Publishing House，1989.（in Chinese）

［6］賈沛然，陳克俊，何力.遠程火箭彈道學［M］.長沙：國防科技大學出版社，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］洪蓓，梁欣欣，辛萬青.固體運載火箭多約束彈道優(yōu)化［J］.導彈與航天運載技術，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.