系統參數對電動力繩系動力學的影響

楊育偉，蔡 洪

(國防科技大學空天科學學院，長沙 410073)

0 引 言

電動力繩系是一種由兩個衛(wèi)星和一條柔性可導電的系繩組成的空間飛行器系統。其中，衛(wèi)星連接于系繩的兩端[1]。在與地磁場及電離層的作用下，電動力繩系能夠感應產生洛倫茲力，從而用于軌道機動[2]。由于具有不消耗燃料、質量輕及功率大等優(yōu)點，電動力繩系技術在空間碎片清除和廢棄衛(wèi)星降軌等空間任務中具有非常廣闊的應用前景[3-8]。而在當前的理論研究中，電動力繩系的擺動動力學和軌道動力學研究是重點關注的問題。

在對電動力繩系的擺動動力學的研究中，周期擺動是關注最多的問題[9-14]。Peláez等[11- 12]首先提出了啞鈴模型且推導了在圓軌道和橢圓軌道上運行的電動力繩系的擺動動力學方程。該方程具有非常強的非線性，在一定的參數范圍內存在周期解。通過Poincaré映射法[11]和Legendre偽譜法[13]等數值方法可以求解得到方程的周期解。由于電動力繩系是處于電動力的周期激勵下，周期擺動從本質上是不穩(wěn)定的。通過計算周期解遷移矩陣的特征值，Peláez和Andrés[12]發(fā)現不穩(wěn)定性隨著電動力參數和軌道偏心率的增長而增長。

降軌動力學是研究電動力繩系軌道動力學的重要問題之一[15-19]。胡長偉等[15]研究了系統與空間等離子體之間的接觸電阻對離軌特性的影響。張健等[16]通過系繩結構的設計，實現了系繩穩(wěn)定與軌道機動控制的解耦，并研究了氣動阻力作用下系統的動力學及控制特性。文獻[17]中運用高斯擾動方程對繩系-納星的軌道動力學進行了建模。文獻[18-19]對保證快速、穩(wěn)定的飛行器降軌的電流控制進行了研究，通過比較不同軌道傾角下的降軌時間得出：電動力繩系的降軌效率與軌道傾角和系繩擺動相關。

在擺動動力學和軌道動力學方程中，存在系統參數：主星質量、子星質量、系繩質量和系繩中電流。在所有對電動力繩系動力學的研究中，上述系統參數對系統擺動動力學和軌道動力學的影響的研究較少。Wang等[20]研究了系統總質量與主星質量的比和子星質量與總質量的比對降軌穩(wěn)定性和效果的影響。然而，每個系統參數對軌道動力學的影響并沒有研究。因此，本文將研究：1)系統參數對擺動動力學的影響；2)系統參數對軌道動力學的影響。

1 系統參數對擺動動力學的影響

1.1 坐標系統

下面介紹文中用到的兩個坐標系。

1)地心慣性坐標系EXYZ：該坐標系以地心E為原點，EX軸指向第一升交點，EZ與地球自轉軸重合，EY符合右手定則。

2)軌道坐標系Oxyz：該坐標系以電動力繩系的質心O為原點，Ox軸沿當地垂線并指向天頂，Oz軸沿軌道動量矩的方向，Oy軸符合右手定則。

1.2 擺動動力學方程和周期解

文中僅研究在留位狀態(tài)階段系統參數對電動力繩系擺動動力學的影響。本階段中，系統的系繩處于張緊狀態(tài)，繩長保持不變。因此，將電動力繩系假設為啞鈴模型，該模型將系繩視為非彈性且無柔性的細桿[11-12]。雖然在空間運行的電動力繩系的系繩長度可達幾千米，但與軌道半徑相比仍然是小量。因此，在建模過程中，假定沿系繩周圍空間的地磁場是均勻的，且取系統質心處的值。為了計算作用于系繩上的電動力，需要對地磁場B建模。本文使用了最簡單的地磁場模型—非傾斜偶極子模型。在軌道坐標系Oxyz中，該模型的各分量(Bx,By,Bz)可表示為

(1)

式中：r為軌道半徑，μm為磁偶極子強度，i為軌道傾角，為近心點角距，f為真近點角。

定義系繩在軌道坐標系的平面Oxy內的投影與Ox軸的夾角為面內角θ，系繩與平面Oxy的夾角為面外角φ，用θ和φ來描述電動力繩系的擺動動力學。對于電動力繩系，主星及子星的尺寸與繩系的長度相比很小；另外，主星與繩系之間、子星與繩系之間的連接均為柔性連接。基于這兩個因素，系統的主星與子星的姿態(tài)對系統擺動的影響在本文中忽略不計。空間中電動力繩系模型和各坐標系如圖1所示。

基于啞鈴模型和非傾斜偶極子磁場模型的假設，運行于橢圓軌道上的電動力繩系的擺動動力學方程可表示為[14]

圖1 電動繩系啞鈴模型和坐標系示意圖Fig.1 Dumbbell model of EDT and coordinate systems

(2)

其中，e為軌道偏心率。在式(2)的推導過程中，將變量對時間的求導變換為了對真近點角f的求導，即式中求導符號(·)表示d( )/df。另外，在推導中，引入了電動力參數ε，其表達式為

(3)

在電動力繩系的運行過程中，雖然在洛倫茲力作用下，系統會產生軌道機動，但是與系統的擺動相比，軌道元素如偏心率、軌道傾角及近地點高度的變化非常緩慢。因此，本文中在研究擺動動力學時，假設系統的軌道保持不變。此時，式(2)中包含參數ε的激勵項具有周期性，在一定的參數范圍內，其具有周期解。但由于式(2)具有很強的非線性，只能采用數值方法來對其進行研究。

采用數值方法進行研究時，需令

(4)

將擺動動力學方程(2)表示為一階形式。

圖2 不同ε值對應的周期解Fig.2 Periodic solution of libration dynamic equation for different ε

擺動動力學方程的周期解，可以通過Poincaré映射等數值方法求解得到[11]。擺動動力學方程的周期解對電動力繩系具有非常重要的意義，其表示系統能量輸入為零的狀態(tài)，且可以作為空間操作的參考軌道。因此，本文研究系統參數對擺動動力學方程的周期解及其特性的影響。

在擺動動力學方程中，各系統參數通過ε來影響該方程。圖2為當i=40°，e=0.2，=0°時，不同ε值對應的周期解。為了不失一般性，算例中i和e的值是隨機選取的，并沒有特殊含義。

圖2顯示，對于不同的軌道參數，隨著ε的增長，擺動動力學方程周期解的擺動幅度也在增大。電動力繩系擺動動力學方程的周期解是不穩(wěn)定的[12]。周期解的不穩(wěn)定性越強，在無控制狀態(tài)下，當系統受到擾動時，周期擺動狀態(tài)會越快地發(fā)生偏離，變?yōu)檗D動運動。而將非周期擺動控制為周期擺動時，此時各參數對應的周期解不穩(wěn)定性越強，控制越難實現。非線性動力學方程周期解的不穩(wěn)定性可以通過周期解的遷移矩陣的特征值的模來表征。

1.3 穩(wěn)定性分析理論

令

(5)

為基于參量σ的n維自治系統的非線性微分方程。當σ=σ0時，如果系統(5)存在周期解，則周期解的局部穩(wěn)定性由該系統的線性化方程決定，即

(6)

式中：A(t)=?f(x,σ)/?x|x=x0(t),σ=σ0為雅閣比矩陣。由于x=x0(t)具有周期性，A(t)為周期時變矩陣。

此時，將非線性系統(5)的周期解穩(wěn)定性判斷近似轉化為線性系統(6)的周期解局部穩(wěn)定性判斷問題。而線性系統周期解的穩(wěn)定性可通過Floquet理論來判斷[21-22]。設n×n的非奇異矩陣Φ(t)為系統(6)的基解矩陣(Fundamental matrix)，即

(7)

另外，存在周期為T的n×n非奇異矩陣P(t)和n×n常值矩陣R滿足

Φ(t)=P(t)etR

(8)

而遷移矩陣(Monodromy matrix)Φ(T)為當t=T、初值條件Φ(0)=I(I為單位陣)時Φ(t)的值；Φ(T)的特征值稱為Floquet特征乘子。

參考文獻[21-22]中有以下結論：方程(6)的周期解漸近穩(wěn)定的必要條件是方程(6)的特征乘子的模都小于1(所有的特征乘子在復平面的單位圓內)。否則，周期解是不穩(wěn)定的。另外，模大于單位1的特征乘子的模越大，周期解的不穩(wěn)定性越強。

對于擺動動力學方程，周期解的遷移矩陣有兩對特征值。其中至少有一個特征值的模大于單位1[12]。本文用遷移矩陣特征值的模即周期解的不穩(wěn)定性來衡量系統參數對電動力繩系擺動動力學的影響。在仿真中采用了文獻[23]介紹的基于精細積分的遷移矩陣計算方法。

1.4 系統參數對擺動動力學的影響仿真

仿真中，分別計算了不同系統參數對應的周期解以及周期解遷移矩陣特征值的模隨系統參數的變化關系。圖3給出了周期解遷移矩陣的特征值的模隨各系統參數的變化關系。針對每個系統參數的仿真中，該系統參數的變化范圍如圖3中對應各子圖中的橫坐標所示，而其余參數的值保持恒定，且取表1中的值。本文假設系繩的長度與質量呈正比，因此圖3(c)中存在兩個橫坐標：系繩質量mt，以及對應的繩長L。另外，在所有的仿真中，軌道元素設為e=0.1，i=25°，=0°。

表1 參考系統的系統參數的值Table 1 Values of system parameters of the baseline system

由圖3可知，對于所有的系統參數，均存在模大于單位1的特征值λ3,4，說明對應的周期解均是不穩(wěn)定的。

2 系統參數對軌道動力學的影響

2.1 軌道動力學方程

為了避免當軌道偏心率或者軌道傾角非常小時造成的奇異性，航天器的軌道運動可以用春分點軌道元素的形式來描述。春分點軌道元素和經典軌道元素之間的關系可表示為

圖3 周期解遷移矩陣特征值的模與系統參數的關系Fig.3 Moduli of eigenvalues of monodromy matrix as functions of system parameters

(9)

其中，a為軌道長半軸，Ω為軌道升交點赤經，λ為平經度，M為平近點角；h，k，p和q分別為用來代替經典軌道元素e，i，Ω和的對應變量。則以春分點軌道元素描述的系統軌道動力學方程為

(10)

(11)

作用在系繩上的洛倫茲力FL可以表示為

(12)

其中，t為與系繩相切的單位向量，它在軌道坐標系中的表達式為

t=cosθcosφi+sinθcosφj+sinφk

(13)

本文僅考慮洛倫茲力和重力作用于系統上，則由洛倫茲力引起的擾動加速度(S,T,W)為

(14)

將式(1)、(14)和(13)代入式(15)，可得

(15)

在洛倫茲力的作用下，電動力繩系可進行降軌。本文通過系統參數對降軌過程的影響來研究系統參數對軌道動力學的影響。在降軌時，軌道長半軸的變化率可表征電動力繩系的效率。在式(10)中，da/dt依賴于加速度分量S和T。令Ta表示項(1+ecosf)T在一個軌道周期內的平均值。在一個軌道周期內，項esinfS的積分遠小于Ta的積分。因此，分量T為決定軌道長半軸變化率的最主要因素。

令

(16)

式(15)中計算T時，地磁場分量Bx所占的比重遠大于分量Bz。則T主要依賴于項-Bzcosθcosφ和ζ。從而，da/dt主要由-Bzcosθcosφ和ζ決定。本文降軌過程中系繩的擺動運動也被考慮進軌道動力學方程中，且系繩以第2節(jié)中的周期解為初始擺動狀態(tài)。第2節(jié)顯示，對于相同的軌道元素，擺動角θ和φ的幅度均隨著ε的增大而增大。因此，項-Bzcosθcosφ的絕對值隨著ε的增大而減小，從而造成T的幅值隨著ε的增大而減小。而T與ζ的關系為正比例關系。因此，Ta與ε呈反比例關系，與ζ呈正比例關系。其中，ε對Ta的影響表示了系繩的擺動對Ta的影響。而Ta的變化趨勢則可反映降軌時間的變化趨勢。下面通過仿真算例來研究系統參數對軌道動力學的影響。

2.2 系統參數對軌道動力學的影響仿真

首先，在考慮了系繩擺動情況下，對運行于橢圓軌道上的電動力繩系的降軌過程進行仿真。本算例中，選取的系統參數的值與表1中所示相同；初始近地點高度為800 km，目標近地點高度為500 km，初始軌道偏心率為0.1，初始軌道傾角為25°。系繩的初始擺動狀態(tài)為對應于系統參數和初始軌道元素的周期解，降軌過程中對應的狀態(tài)則由擺動動力學方程積分得到。圖4給出了降軌過程中軌道近地點高度、偏心率及軌道傾角隨時間的變化關系。圖5給出了降軌期間系繩面內和面外擺動角的變化。

圖4顯示，電動力繩系的降軌過程需要約9天的時間。在降軌過程中，軌道偏心率呈增大趨勢，軌道傾角呈減小趨勢。但是與近地點高度的明顯變化相比，這兩個軌道元素的變化幅度均很小。

圖5顯示，降軌過程中電動力繩系的擺動周期約為1.5～2.0 h。在仿真過程中，繩系的擺動考慮了軌道元素的變化，即圖5中擺動角的計算中計入了軌道元素的實時變化。在擺動動力學方程(2)中，軌道近點高度的變化對該方程沒有影響，只有軌道傾角和軌道偏心率的變化對周期解產生擾動。在一個擺動周期內，圖4中的軌道傾角變化約0.00375°，軌道偏心率變化約0.0000225；即在整個降軌過程中，軌道偏心率和軌道傾角的變化非常小且非常緩慢。雖然周期解是不穩(wěn)定的，但圖5的結果顯示，在降軌過程中軌道元素存在小的擾動下，繩系在長時間內仍然為近似周期擺動。這也反過來說明，在單獨研究系統的擺動動力學過程中將軌道元素假設為不變是合理的。

圖4 軌道近地點高度、偏心率及軌道傾角隨時間的變化關系Fig.4 Time histories of the orbital altitude of the perigee, eccentricity and inclination of EDT during deorbit

圖5 降軌期間擺動角的變化Fig.5 Time histories of the in-plane and out-of-plane libration angles during deorbit

然后，利用仿真算例來研究各系統參數對軌道動力學的影響。目前電動力繩系最主要的應用是對廢棄衛(wèi)星進行降軌，降軌效率是最重要、最需要考察的因素，軌道高度是最明顯的變化量。因此，本文中僅關注近地點高度的變化，用系統的近地點高度從800 km降到500 km所耗時間來衡量降軌效率。

圖6給出了降軌時間隨各系統參數的變化關系。各子圖中橫坐標為所研究的系統參數的變化范圍，而其余參數的值與表1中所示相同。所有的仿真中，初始的e和i設為e=0.1，i=25°；且在飛行過程中考慮繩系的擺動。初始的擺動狀態(tài)則為對應于初始e,i及ε的周期解。

圖6(a)顯示，降軌時間隨m1的增大而增大。m1的增大導致擺動動力學方程中的ε增大，在軌道動力學方程中，ζ的值隨m1的增大而減小，因此分別正比于ζ和反比于ε的Ta的值隨著m1的增大而減小；在減小的Ta作用下，降軌過程中的da/dt的絕對值減小，即降軌效率變慢。圖6(a)中結果說明，增長的m1對系統的軌道動力學具有消極的影響。

圖6 降軌時間隨各系統參數的變化關系Fig.6 Values of deorbit time as functions of system parameters

圖6(b)顯示，降軌時間隨m2的增大先減小后增大。這是因為隨著m2的增大，擺動動力學方程中ε的值減小；而軌道動力學方程中ζ的值則隨著m2的增大而減小；這兩個因素導致Ta的絕對值隨著m2的增大先減小后增大。Ta的絕對值越大，降軌效率越高，降軌時間越短；Ta的絕對值越小，降軌效率越低，降軌時間越長。圖6(b)中的結果說明，增長的m2對系統的軌道動力學的影響為先積極后消極。

圖6(c)顯示，降軌時間隨mt的增大而變短。這是因為隨著mt的增大，擺動動力學方程中的ε減小，而軌道動力學方程中的ζ則增大，從而導致Ta的絕對值隨著mt的增大而增大。圖6(c)結果說明，增長的mt對系統的軌道動力學的影響為積極的。本質上，這是因為mt的增大意味著系繩的增長，從而提高了降軌效率。

3 討 論

本節(jié)對系統參數對擺動動力學和軌道動力學的影響進行綜合討論，給出了關于每個系統參數的設計建議并提出了降軌策略。系統參數對擺動動力學和軌道動力學的影響總結于表2中。

表2 系統參數對擺動動力學和軌道動力學的影響Table 2 Summary of effects of the system parameters on libration dynamics and orbital dynamics

綜合仿真結果可知，增長的m1對系統的擺動動力學和軌道動力學均有消極的影響。另外，增大的m1意味著發(fā)射成本的增加。因此，主星質量m1設計得越小越好。

而m2較小的增長就會造成擺動動力學不穩(wěn)定性很大程度地降低，且在所有的系統參數中，增長的m2對于減小周期解的不穩(wěn)定性效果最為明顯。盡管增長的m2對于軌道動力學的影響為先積極后消極，但是最長的降軌時間與最短的降軌時間的差僅僅為3天，相較其他參數對降軌動力學的影響，是非常小的。相比其對擺動動力學的積極的影響，其對軌道動力學消極的影響可以忽略。另外，能夠對擺動動力學產生明顯效果的增長的m2也僅僅幾千克的量級，這對于發(fā)射成本的產生的影響也比較小。因此，綜合考慮，當進行系統設計時，子星質量m2越大越好。

mt的增長對系統的擺動動力學和軌道動力學均有積極的影響。其中，隨著mt的增長降軌時間明顯減少。mt的增長對降軌動力學產生的影響相對較小，在考慮發(fā)射成本時可以忽略。但是，mt的增長意味著系繩長度的增加，這樣就增加了系統的展開難度。因此，在不影響系統展開的前提下，系繩的質量mt設計得越大越好。

電動力繩系最主要的功能是進行降軌，因此有必要綜合考慮各種對降軌產生影響的因素，對應于降軌的整個系統進行設計，并提出合理的降軌策略。

在進行系統設計時，主星的質量設計要盡量小；子星的質量設計在滿足發(fā)射成本約束、質量小于主星質量約束以及其它約束的情況下，質量盡可能大；繩系長度的設計則是在滿足展開機構約束條件下，長度盡可能長。三個因素中，繩系長度的變化對降軌的影響最為明顯，且最具有可行性，因此，在進行系統設計時，繩系的長度為最主要的設計變量。

在降軌過程中，繩系中的電流和繩系的姿態(tài)擺動均對系統的軌道動力學具有影響。因此，降軌策略的制定需主要考慮這兩個因素。首先，由于常值電流情況下，系統存在規(guī)律的周期姿態(tài)擺動，因此降軌過程中需將繩系中電流控制為常值。而常值電流的大小則需要綜合考慮姿態(tài)擺動的不穩(wěn)定性以及電流產生能力來確定。在滿足穩(wěn)定性控制等約束范圍內，電流的值越大越好。當系統的各組件質量、繩系中的電流大小確定后，以對應的周期解為初始擺動狀態(tài)。在運行的過程中，當干擾力和軌道元素攝動等因素使得繩系的擺動偏離周期解達到一定的程度時，對繩系的擺動進行控制，使其回到此時對應的周期解上。繩系擺動的控制可通過控制繩系中的電流，或者安裝在子星上的小推力發(fā)動機等方法實現。當繩系擺動恢復到周期解上時，停止控制；而當繩系擺動再次偏離周期解達到一定程度時，再次施加控制，這樣直到降軌過程完成。

4 結 論

本文研究了主星質量、子星質量、系繩質量以及系繩中的電流對電動力繩系擺動動力學和軌道動力學的影響。通過各系統參數對周期解不穩(wěn)定性的影響研究了系統參數對擺動動力學的影響。通過各系統參數對降軌時間的影響研究了系統參數對軌道動力學的影響。分析結果表明：主星質量的增大對擺動動力學和軌道動力學均有消極的影響；子星質量的增大對擺動動力學具有積極的影響，而對軌道動力學的影響則為先積極后消極；系繩質量的增大對擺動動力學和軌道動力學均有積極的影響；增大的電流對擺動動力學具有消極的影響，而對軌道動力學則有積極的影響。綜合考慮各參數對系統擺動動力學和軌道動力學的影響程度，以及實際的發(fā)射、空間的展開等因素，可為各系統參數的設計以及降軌策略的制定提供理論參考和建議。