張意國，趙長見，趙俊鋒，賈生偉

（中國運載火箭技術(shù)研究院，北京，100076）

0 引 言

對于大部分運載火箭和彈道導彈而言，助推器工作結(jié)束后需要進行分離并被拋掉以降低結(jié)構(gòu)重量。對于一般飛行器而言，在子級助推器分離時刻飛行速度不足以達到第一宇宙速度，使得分離后的助推器在地球引力的作用下將落入地球表面。助推器再入過程由于沒有主動控制力的作用，在重力及非線性氣動力的作用下將出現(xiàn)復雜的姿態(tài)運動。對于助推器再入姿態(tài)運動模態(tài)（即無控狀態(tài)下姿態(tài)周期性運動規(guī)律）[1]進行預示，有助于對助推器回收系統(tǒng)進行設(shè)計，同時也有助于發(fā)射場安控工作的組織實施。王景國等[2]人采用了傳統(tǒng)的基于3DOF模型的落點預測方法對火箭整流罩的殘骸落點進行估計，其中通過大量的試驗數(shù)據(jù)進行聚類分析結(jié)合工程經(jīng)驗對阻力系數(shù)和參考面積這兩個關(guān)鍵彈道參數(shù)進行了計算。陳彬[3]等人對助推器再入運動的單通道運動特性進行了分析，分析了助推器再入過程姿態(tài)運動的穩(wěn)定性。

為全面準確地描述助推器再入自由運動過程的姿態(tài)運動特性，有必要開展其三通道運動模態(tài)的研究。對助推器三通道運動模態(tài)的研究建立在助推器再入三通道耦合非線性模型的基礎(chǔ)上，相較于單通道線性化模型[3]，三通道運動學模型能夠提供更高的模型精度，同時也存在不止一組平衡點，對助推器再入動力學特性的描述可以更準確。但也因為非線性模型強耦合以及強非線性的特點，使得基于小擾動理論的線性系統(tǒng)理論對一些系統(tǒng)特性無法展開準確地分析。

本文采用了基于演化搜索算法的改進郭濤算法，解決高階、強耦合非線性動力學模型的多平衡點求解問題[4～7]，并在此基礎(chǔ)上對助推器再入過程平衡點的穩(wěn)定性及演變規(guī)律進行了分析。

1 助推器再入非線性模型

在再入彈頭和自旋穩(wěn)定飛行器動力學的研究中，常使用基于經(jīng)典歐拉角描述的運動學模型，以便于更好地描述大姿態(tài)角運動下的運動規(guī)律[8,9]。對于再入助推器而言，其存在大范圍的姿態(tài)角運動，且姿態(tài)角運動存在高動態(tài)特性，適合于使用基于經(jīng)典歐拉角及其角速度描述的姿態(tài)動力學模型來進行動力學分析工作。

經(jīng)典歐拉角的定義如圖1所示，包括進動角、章動角及自旋角，是描述助推器體坐標系與準速度坐標系之間轉(zhuǎn)換關(guān)系的一組歐拉角，相應坐標系描述如下：

a）助推器體坐標系（Oxbybzb）。

坐標系原點O位于質(zhì)心，O xb軸沿助推器縱軸向前為正，Oyb在對稱平面內(nèi)指向上為正，Ozb軸與其余兩軸滿足右手定則。

b）準速度坐標系（Ox5y5z5）。

坐標系原點O位于質(zhì)心，O x5軸沿速度方向向前為正，O y5位于包含Ox5的豎直平面內(nèi)向上為正，O z5軸與其余兩軸滿足右手定則。

由準速度系轉(zhuǎn)到助推器體坐標系為1-3-1轉(zhuǎn)序：首先Ox5y5z5坐標系沿Ox5旋轉(zhuǎn)ζ（進動角）至Ox5y' z'坐標系；然后繞Oz'軸旋轉(zhuǎn)δ（章動角）至Oxby'' z'坐標系；最后繞Oxb軸旋轉(zhuǎn)λ（自轉(zhuǎn)角）至Oxbybzb坐標系。

圖1 經(jīng)典歐拉角定義Fig.1 Classical Euler Angle Definitionζ—進動角；δ—章動角；λ—自轉(zhuǎn)角

助推器再入飛行過程沒有控制力，只在重力和氣動力的作用下自由運動，而空氣動力一般在彈體系下描述，是與攻角α及側(cè)滑角β相關(guān)的變量，根據(jù)坐標轉(zhuǎn)換關(guān)系可得到α，β以ζ，δ，λ表示的表達式：

考慮到助推器姿態(tài)運動相對質(zhì)心運動為快時變運動，在分析姿態(tài)運動時可以認為助推器速度的大小及方向保持不變，準速度系可視為慣性坐標系，則再入助推器的角速度（xω，yω，zω）可以用δ、ζ及λ的角速度之和表示：

式中xω，yω，zω為彈體系下歐拉角的3個分量；ζ˙為進動角速度；λ˙為自旋角速度；δ˙為進動角速度。

對式（2）進行求導并進行轉(zhuǎn)換，可以得到以歐拉角表示的助推器非線性動力學模型：

式中xω˙，yω˙，zω˙為彈體系下歐拉角速度bω˙的3個分量；δ˙˙為進動角加速度；ζ˙˙為章動角加速度；λ˙˙為自旋角加速度。

可根據(jù)歐拉方程求得，矢量形式為

式中 I為彈體轉(zhuǎn)動慣量矢量；ω為彈體角速度矢量；M為合外力矩矢量。

2 非線性模型平衡點求解

對于助推器再入飛行無控姿態(tài)動力學系統(tǒng)，為了得到系統(tǒng)在固定參數(shù)下的平衡點，需要將模型中的所有狀態(tài)變量關(guān)于時間的導數(shù)設(shè)為0。以上標“P”表示平衡點處的各狀態(tài)變量，則得到平衡點處的運動方程：

式中 Ip為彈體平衡點轉(zhuǎn)動慣量矢量； ωp為彈體平衡點角速度矢量； Mp為平衡點合外力矩矢量。

式（6）即為平衡狀態(tài)的運動方程，展開可得到平衡方程的標量形式為

3 郭濤算法

針對高階強非線性方程組，在理論上不存在解析解，對于該問題的求解一般采用數(shù)值方法。比較常用的數(shù)值方法包括牛頓法和遺傳算法等，但這些算法對初值敏感，求解結(jié)果具有局部性，無法實現(xiàn)全局多解的求解，甚至在一些非線性較強的方程組求解過程中會出現(xiàn)直接發(fā)散的現(xiàn)象。所以有必要研究對于該類方程組的有效求解方法。

武漢大學的郭濤提出了一種演化搜索算法，其在求解非線性優(yōu)化問題及多解問題時有獨特的優(yōu)勢，能夠保證解的全局性，同時收斂速度較快，是一種適合解決本文問題的方法，其特點為：

a）提出了一種基于全局的群體搜索策略，理論上能夠保證解的全局性。

b）提出了種群迭代優(yōu)化策略，采用遺傳算法中父代雜交與隨機子空間搜索并行策略，保證了在種群優(yōu)化過程中兼顧搜索的遍歷性。

c）提出了末等淘汰策略，種群迭代中每次只淘汰最劣等個體，能夠較好地控制種群的整體性。

郭濤算法同時在理論上保證了全局性及遍歷性，這樣就保證了算法能夠依概率一收斂到全局最優(yōu)解。

該方法的主要求解策略及流程如圖2所示。

圖2 郭濤算法的求解過程Fig.2 Guo Tao Algorithm Solution Process

a）步驟一：問題的描述。

郭濤算法求解的第1步，將常見的非線性多維方程組多解問題抽象為有約束下的最優(yōu)化問題。

設(shè)非線性方程組為

假定該方程組有解，其求解問題可以轉(zhuǎn)化為

其中，S為搜索空間， S ?Rn，當fitness(x):S→R稱為目標函數(shù)，fitness(x)= 0 時，x為方程組的精確解，通常S是一個n維長方體，即由所確定；Φ為可行點集或可行區(qū)域。

b）步驟二：確定編碼方式。

確定問題的數(shù)學模型后，需要制定研究對象的編碼方式，對于復雜的非線性問題，采用浮點數(shù)編碼能夠提高算法求解的準確程度。

c）步驟三：種群初始化。

d）步驟四：個體評價。

本文的優(yōu)化問題可歸納為連續(xù)優(yōu)化問題，因為約束及解空間在實數(shù)空間上是處處連續(xù)的，因此針對該問題的求解過程能夠保證較好的局域性及精確程度。

用邏輯函數(shù) b etter(x1, x2)表示測試點 x1優(yōu)于測試點：

定義：

用定義的邏輯運算函數(shù)，求解種群中的最優(yōu)個體xbest和 最 差 值xworst， 使 得 ： ?x∈P，better(xbest, x) ,better(x, xworst)；如果 b etter(xworst,xbest)，轉(zhuǎn)到步驟八。

e）步驟五：子代優(yōu)化（父體雜交）。

以相同概率選取上述兩種方法中的一種決定搜索空間V，并在V中隨機選取一點x。

f）步驟六：劣汰策略。

如果better(x, xworst)，則xworst= x 。

g）步驟七：轉(zhuǎn)到步驟四。

h）步驟八：輸出 xbest。

i）步驟九：結(jié)束。

4 助推器再入三通道動力學特性分析

4.1 平衡點求解

求解平衡狀態(tài)非線性方程組時，氣動力系數(shù)通過攻角α和側(cè)滑角β在（-180°，+180°）范圍內(nèi)插值得到，系統(tǒng)具有定義域范圍大和參數(shù)強非線性的特點（各特征馬赫數(shù)下法向力系數(shù)NC和俯仰力矩系數(shù)MZC的變化規(guī)律見圖 3、圖 4）。傳統(tǒng)解非線性方程組的方法如半解析法和牛頓法只能解決小定義域和參數(shù)定常的問題，不適用于本文的模型。郭濤算法采用變空間搜索及劣勢淘汰策略，對方程組解的搜索具有全局性，不受定義域大小的限制。同時由于郭濤算法基于演化搜索的思想，只要求模型參數(shù)可求，不要求其定?；蚩晌ⅲ沟霉鶟惴ň哂懈鼜姷倪m應性，適用于本文的參數(shù)強非線性模型。

圖3 NC 隨攻角變化Fig.3 CN’s Relationship to the Change in Angle of Attack

圖4 MZC 隨攻角變化Fig.4 CMZ’s Relationship to the Change in Angle of Attack

根據(jù)非線性微分方程組的特性和助推器氣動軸對稱的特點，可知若{ζ ˙ ,δ, λ}為方程組的一個解，那么{ζ ˙ ,? δ,π+λ}，{? ζ ˙, δ, ?λ}，{?ζ˙,? δ, π ?λ}也是方程組的解，且在這 4個解處，系統(tǒng)具有相似的姿態(tài)運動規(guī)律，可以認為它們代表1個解。因此在下文中僅給出一種解的形式作為相同動態(tài)特性解的代表，且如{0 ,π,0} 等顯而易見的解也不再記錄。最后通過郭濤算法仿真計算得到助推器在不同特征馬赫數(shù)下對應的非線性姿態(tài)運動的平衡點，見表1。

表1 郭濤算法解不同特征Ma下對應非線性模型平衡點Tab.1 Guo Tao Algorithm’s Solutions Depend on Diferrent Ma

續(xù)表1

由表1可知，利用郭濤算法求解出每個特征馬赫數(shù)下助推器姿態(tài)運動非線性方程組存在4組解，每組解對應的適應值函數(shù)均接近于零，表明算法精度較高，結(jié)果可信。

通過經(jīng)典歐拉角的定義可知：ζ˙表征助推器的旋轉(zhuǎn)速度；δ表征助推器縱軸與速度軸夾角，等價于總攻角；λ表征助推器縱平面與總攻角平面夾角。通過表 1對不同特征馬赫數(shù)下的平衡點分析可以得到平衡點隨馬赫數(shù)變化的規(guī)律，見圖5、圖6。由圖5以及圖6的分析結(jié)果知，不同于線性系統(tǒng)，助推器非線性系統(tǒng)存在多組平衡點。不考慮對稱分布的相似平衡點，在再入馬赫數(shù) M a∈ (1 ,9)的情況下分布著 4組平衡點：其中平衡點1對應δ=0°；平衡點2對應δ≈90°，隨著馬赫數(shù)的增加章動角變化不大，而相應的進動角速度隨馬赫數(shù)增加而增加；平衡點3、4對應δ隨馬赫數(shù)的增加而增加，由δ≈110°向δ≈160°變化。

圖5 平衡點對應章動角隨馬赫數(shù)變化規(guī)律Fig.5 Pδ’s Relationship to the Change in Angle of Attack

圖6 平衡點對應進動角速度隨馬赫數(shù)變化規(guī)律Fig.6 ζ ˙p ’s Relationship to the Change in Angle of Attack

4.2 平衡點穩(wěn)定性分析

助推器再入姿態(tài)動力學是強耦合高維非線性系統(tǒng)，很難直接根據(jù)李雅普諾夫穩(wěn)定性的定義來判斷其穩(wěn)定性。對于此類非線性系統(tǒng)，李雅普諾夫第一法（或稱間接法）在穩(wěn)定性的分析中具有著廣泛的應用。

李雅普諾夫第一法（間接法）：

若x=0是非線性系統(tǒng) x˙ =f( x)的一個平衡點，其中f →Rn是連續(xù)可微的，設(shè)雅克比矩陣為

如果D的所有特征值都滿足 R eλi< 0，則原點是漸近穩(wěn)定的；

如果D至少有一個特征值滿足 R eλi> 0，則原點是不穩(wěn)定的；

如果有某一解的特征根實部在零點而其余特征根的實部為負值，此時則無法判斷該平衡點的穩(wěn)定性。

運用Lyapunov的第1種方法求解本文再入非線性系統(tǒng)平衡點的穩(wěn)定性，可分析得到平衡點的穩(wěn)定性隨馬赫數(shù)變化規(guī)律，可得到再入助推器姿態(tài)動力學平衡點穩(wěn)定性的變化規(guī)律，見圖7。

圖7 δ解曲線的穩(wěn)定性Fig.7 Stability of the δ’s Solution

由圖7可知，在助推器再入飛行絕大部分馬赫數(shù)域內(nèi)，非線性系統(tǒng)具有一組穩(wěn)定平衡點和3組不穩(wěn)定平衡點，這意味著在合適的初始條件下系統(tǒng)會被鎖定在這一平衡狀態(tài)下。而在Ma∈（0.76,0.84）這一馬赫數(shù)域內(nèi)，系統(tǒng)的2組不穩(wěn)定解曲線發(fā)生了性態(tài)變化，演變成穩(wěn)定平衡解曲線，此時系統(tǒng)存在3組穩(wěn)定平衡點，意味著在合適的初始條件下系統(tǒng)可能被鎖定在其中任意一種平衡狀態(tài)下。此處以Ma=0.8的工況為例，將不同初始狀態(tài)對應的相軌跡在δ-δ˙平面上進行投影，所得結(jié)果如圖8所示。

a）相軌跡在δ-δ˙平面上進行投影

續(xù)圖8

圖8 中不同曲線代表的是不同初始條件對應的相圖在δ-δ˙平面的投影，從圖 8中可以看到在仿真工況下的吸引子包括2個漸進穩(wěn)定的平衡點和一個穩(wěn)定的極限環(huán)，不同的初始條件會使得助推器的姿態(tài)收斂到不同的吸引子，驗證了基于演化搜索算法求解的平衡狀態(tài)的正確性。

將不同初始狀態(tài)對應的相軌跡在t-ζ˙平面上進行投影，所得結(jié)果如圖9所示。

圖9 解軌線在t-ζ˙平面投影Fig.9 Solve for the Projection onto t-ζ˙ Plane

由圖9可知，在Ma=0.8的工況下，在進動角速度存在初始偏差的情況下，均能收斂到平衡點對應的穩(wěn)定狀態(tài)（平衡點2和4）或在平衡點對應穩(wěn)定狀態(tài)附近震蕩（平衡點3），其中平衡點2對應的穩(wěn)定進動角速度ζ˙為7.006 21，平衡點3和4對應ζ˙為0。再次驗證了平衡點解算的正確性。

5 結(jié) 論

本文采用基于演化搜索策略的郭濤算法，解決了助推器再入飛行模態(tài)的預示問題。首先推導了以經(jīng)典歐拉角描述的姿態(tài)運動模型，使研究結(jié)論更直觀，模型更簡潔。利用該模型，進一步推導了五維姿態(tài)動力學平衡點方程組，并采用基于演化搜索算法的改進郭濤算法解決了非線性方程組的多解問題。進一步，利用Lyapunov第一法對平衡點穩(wěn)定性進行了分析。最后通過了時域分析的驗證。結(jié)果表明在再入過程中的特定馬赫數(shù)域內(nèi)（Ma=（0.76，0.84））存在3組飛行模態(tài)，一組模態(tài)對應高速旋轉(zhuǎn)運動，其余兩組模態(tài)對應姿態(tài)穩(wěn)定運動。進入該飛行速域時的不同初始狀態(tài)將使助推器姿態(tài)運動收斂到不同的飛行模態(tài)。

本文提出了助推器無控飛行模態(tài)的預示方法，獲得了再入飛行速域內(nèi)助推器可能出現(xiàn)的多個穩(wěn)定飛行模態(tài)，為進一步設(shè)計助推器回收系統(tǒng)或精細化設(shè)計安控區(qū)提供了有力支撐。但本文并未對各模態(tài)的穩(wěn)定域進行分析，如何解決高維非線性模型穩(wěn)定域的確定將是下一步研究的重點。