火箭末級電動繩系離軌動力學(xué)建模與分析

翟光,宋國梁,楊少偉

(北京理工大學(xué)宇航學(xué)院,北京100081)

1 引言

近年來,空間碎片的主動移除技術(shù)逐漸受到各國航天機構(gòu)的關(guān)注[1-2],國內(nèi)外許多學(xué)者提出并研究了不同的空間碎片主動移除方法[3]。采用電動繩系 (Electro-Dynamic Tether)對空間碎片實施離軌操作,具有無需燃料消耗的特殊優(yōu)勢,因此受到了各國重視。針對電動繩系離軌過程的控制問題,國內(nèi)學(xué)者張健[4-5]等人研究了電動繩系軌道運動與繩系運動的耦合問題,并采用結(jié)構(gòu)設(shè)計方法,提出了解耦控制策略;針對電動繩系橫線振動問題,徐大富[6-7]通過建立電動繩系的偏微分方程和邊界條件,采用有限元方法分析了系統(tǒng)的橫向振動規(guī)律;電動繩系的自主離軌速率是較多學(xué)者關(guān)注的問題,Steven G[8]等人在控制電流變化的基礎(chǔ)上,設(shè)計了電動繩系軌道機動的變軌控制策略;R.Zhong[8]基于高斯攝動方程建立了電動繩系納星系統(tǒng)的軌道動力學(xué)模型,并分析了在不同軌道傾角條件下系統(tǒng)的離軌特性。Claudio[9]建立了電動繩系離軌時間的理論計算模型,并結(jié)合第十一代國際地磁場模型進行了仿真驗證。然而,由于受安培力長期作用,電動繩系的軌道和姿態(tài)運動存在嚴重耦合,安培力矩可通過改變姿態(tài)運動的方式,直接影響電動繩系的軌道變化速度。在目前已有文獻當(dāng)中,多數(shù)學(xué)者均在假設(shè)繩系方向與軌道矢徑平行的基礎(chǔ)上,開展相關(guān)建模與分析工作,但事實上,在不施加控制的條件下,電動繩系與軌道矢徑方向存在著夾角,且是動態(tài)變化的,因此,其離軌過程運動特性分析結(jié)果具有較大的誤差。

本文針對火箭末級等大型空間碎片離軌問題,在考慮安培力矩的基礎(chǔ)上,首先建立電動繩系釋放過程的非線性姿軌耦合動力學(xué)模型,充分反映了加電之后安培力對電動繩系姿態(tài)運動的影響,隨后結(jié)合第十一代國際地磁場模型,建立了電動繩系離軌過程的動力學(xué)數(shù)值仿真模型,并通過仿真分析了離軌過程中電動繩系姿態(tài)和軌道變化情況,由此獲得了不同初始軌道條件下電動繩系系統(tǒng)的離軌時間。

2 電動繩系姿軌耦合動力學(xué)建模

2.1 建模參考坐標(biāo)系

首先建立電動繩系的軌道坐標(biāo)系和本體固連坐標(biāo)系。如圖1所示,軌道坐標(biāo)系原點位于系統(tǒng)的質(zhì)心,xo軸和zo軸分別指向系統(tǒng)軌道速度方向和地心方向;本體固連坐標(biāo)系為最大主軸慣量坐標(biāo)系,原點固連于系統(tǒng)的質(zhì)心,zb軸由飛行平臺質(zhì)心指向末端質(zhì)量質(zhì)心,xb軸方向由式 (1)確定,yb軸根據(jù)右手螺旋法則確定。

建模作如下假設(shè):重力場為理想的中心力場;忽略地球扁率和空間環(huán)境攝動因素的影響;火箭末級和末端質(zhì)量均為質(zhì)點,電動繩系末端2個連接點都位于連接對象的質(zhì)心。

圖1 電動繩系動力學(xué)建模參照坐標(biāo)系Fig.1 The EDT reference frame for dynamic modeling

2.2 系統(tǒng)的動能

如圖2所示,根據(jù)質(zhì)心的定義有:

將上式中x替換為y、z亦成立。式 (2)中m為末端質(zhì)量大小,M為火箭末級質(zhì)量,ρ電動繩系的線密度,l?、l為電動繩系總長度和展開長度;x1和x2分別火箭末級和末端質(zhì)量為在軌道坐標(biāo)系內(nèi)的x軸坐標(biāo)分量。

圖2 電動繩系構(gòu)型示意圖Fig.2 The configuration of EDT system

根據(jù)質(zhì)心定義可得:

電動繩系在釋放過程當(dāng)中,整個系統(tǒng)的動能由3部分組成,其中包括系統(tǒng)質(zhì)心沿軌道運動所產(chǎn)生的動能T1;系統(tǒng)繞質(zhì)心轉(zhuǎn)動所產(chǎn)生的動能T2;連接電動繩系長度變化所產(chǎn)生的動能T3。為節(jié)約篇幅,推導(dǎo)過程不再贅述,這里直接給出表達式如下:

2.3 系統(tǒng)的勢能

對于質(zhì)量微元dm,其在重力場中的勢能為:

末端質(zhì)量和火箭末級的勢能可分別表示為:

μ為地球引力常數(shù),ρ1和ρ2兩段繩系在軌道坐標(biāo)系內(nèi)的矢量,r表示軌道坐標(biāo)系內(nèi)的軌道位置矢量,根據(jù)2質(zhì)點在軌道坐標(biāo)系中的位置關(guān)系,并忽略高階項,對式 (8)分母部分進行泰勒展開,則有:

式中,r=‖r‖,α和β分別表示電動繩系的面內(nèi)外擺角。將式 (9)和式 (10)代入到式(8)中,由此可以得到火箭末級和末端質(zhì)量的勢能函數(shù),其表達式可以表示為

采用同樣的方法可以得到電動繩系的勢能函數(shù)為:

此外,在實際電動繩系釋放過程當(dāng)中,電動繩系在一定的拉伸作用下發(fā)生彈性變形。引入電動繩系的應(yīng)變系數(shù)ξ,E為電動繩系的彈性模量,A為橫截面積,則系統(tǒng)的彈性應(yīng)變勢能可以表示為:

系統(tǒng)總勢能函數(shù)可由V=V1+V2+V3計算得到。

2.4 系統(tǒng)的拉格朗日函數(shù)

根據(jù)拉格朗日定理,系統(tǒng)的拉格朗日函數(shù)為:

將系統(tǒng)的動能函數(shù)和勢能函數(shù)代入上式可得系統(tǒng)的拉格朗日函數(shù)具有以下的形式 (見式14):

由于模型考慮了電動繩系的質(zhì)量,系統(tǒng)的拉格朗日函數(shù)十分復(fù)雜。在實際的工程當(dāng)中,由于電動繩系質(zhì)量遠遠小于系統(tǒng)其他部分質(zhì)量,因此可以忽略不計,即可認為ρ=0,此外,忽略系統(tǒng)的柔性,此時拉格朗日函數(shù)中的彈性勢能項可略去不計,進一步定義歸一化廣義外力:

式中,FB(t)、τB(t)和分別表示廣義外力和外力矩,其中包括安培力和力矩。根據(jù)拉格朗日定理,按式 (12)對拉格朗日函數(shù)中的廣義坐標(biāo)求導(dǎo),

式中,qj和Qj分別表示廣義坐標(biāo)和廣義外力,則可導(dǎo)出系統(tǒng)的剛性動力學(xué)模型為:

動力學(xué)方程是在連接電動繩系剛性假設(shè)的情況下推導(dǎo)建立的,而系統(tǒng)的剛性假設(shè)只有在連接電動繩系處于張緊狀態(tài)的時候才能夠成立,若電動繩系處于松弛狀態(tài)時,末端質(zhì)量和火箭上面級之間的動力學(xué)耦合關(guān)系不再成立,即式能夠描述系統(tǒng)動力學(xué)特性的前提條件是連接電動繩系處于張緊狀態(tài),即T^(t)>0。

2.5 模型的軌道面內(nèi)外運動解耦

觀察式(17),當(dāng)β和β·的初值為零時,其他3個變量方向上的運動并不會誘發(fā)面外擺角的運動,當(dāng)忽略面外擺角后,模型可以進一步簡化,其中剛性動力學(xué)模型可以簡化為:

在不考慮安培力的條件下,式 (18)中第1個方程可以表示為

方程中左側(cè)第2項為擾動項,這說明電動繩系長度或者姿態(tài)變化時,會帶來實際軌道相對于開普勒軌道的偏移。嚴格意義上來講,電動繩系系統(tǒng)質(zhì)心的軌道并不是精確的開普勒軌道,但是由于擾動項很小,在短時間運動過程當(dāng)中可以忽略不計。

3 系統(tǒng)的振蕩運動

3.1 基本模型

國際參考地磁場 (IGRF)是描述地球主磁場空間結(jié)構(gòu)和時間變化的通用模型。地磁場磁位勢滿足拉普拉斯方程,它在球坐標(biāo)系中用球諧函數(shù)表示:

地磁場強度B為磁勢函數(shù)V的負梯度,即

式中,re為地球參考球赤道半徑;r、λ、φ分別為 EDT的地心距、地心經(jīng)度和緯度;Pnm為schmidt函數(shù);為磁場球諧函數(shù)的系數(shù),IGRF模型提供相應(yīng)的高斯系數(shù)表;為n次m階schmidt準(zhǔn)歸一化締合勒讓德函數(shù)。

我國《國家安全法》第52條規(guī)定：“國家安全機關(guān)、公安機關(guān)、有關(guān)軍事機關(guān)根據(jù)職責(zé)分工，依法搜集涉及國家安全的情報信息。”特別在第四章第二節(jié)“情報信息”第51條～第54條中，分別規(guī)定了“情報工作制度”“各部門搜集上報情報信息職責(zé)”“情報信息工作運用現(xiàn)代科技手段和加強研判分析”“情報信息的報送要求”等內(nèi)容?！霸摲ǖ牧⒎Ｊ匠尸F(xiàn)‘原始型’（即分散式立法）向‘混合型Ⅰ’（即分散式＋專門式立法）再向‘混合型Ⅱ’（即分散式＋專門式＋綜合式立法）發(fā)展的趨勢”。《國家安全法》是我國國家安全領(lǐng)域的基礎(chǔ)性法律，也是統(tǒng)籌、引領(lǐng)國家安全領(lǐng)域立法工作的綜合性法律之一，該法涵蓋了國家安全多個領(lǐng)域。

3.2 地磁場傾斜偶極子簡化模型

n=1、m=1時稱傾斜偶極子模型。本文考慮了磁軸中心與地球中心不重合的狀況,采用了偏心傾斜偶極子模型,此時磁場強度在慣性系下為:

式中,Beq是偶極子模型中rE處的磁場密度;um是磁偶極子單位矢量;rd是EDT相對于磁偶極子中心的距離,ud是EDT相對于磁偶極子中心的單位矢量。磁偶極子密度Beq可以直接由高斯系數(shù)表示為:

若初始時刻磁偶極子軸在地心慣性坐標(biāo)系的x-z平面內(nèi),那么磁偶極子單位矢量為:

其中Ω是地球自旋角速度;β是磁偶極子軸相對于地球自轉(zhuǎn)軸的傾斜角,對于傾斜偶極子模型,其值可以由下式得到:

根據(jù)2015年的國際參考地磁場模型,磁偶極子中心在地心地固坐標(biāo)系的坐標(biāo) (單位是km)以及偶極矩向量的經(jīng)度為:

根據(jù)高斯定理,它們由下列各式求得:

磁偶極子中心在赤道慣性坐標(biāo)系內(nèi)可以表示為:

其中:

式中,rsc為EDT在慣性系內(nèi)的位置矢量。

4 電動繩系軌道運動分析

4.1 安培力作用下的軌道攝動

電動繩系系統(tǒng)軌道動力學(xué)方程為:

式中,r是位置矢徑,v=;fe為安培力攝動下的加速度。電動力纜繩穿過地磁場,可以產(chǎn)生的電流大小為:

纜繩系統(tǒng)的總電阻為R,電流在地磁場中的運動感應(yīng)產(chǎn)生作用于纜繩的洛侖茲力為:

將式 (32)式代入式 (31),最終得到電動繩系在安培力作用下的攝動方程:

4.2 安培力矩作用下的面內(nèi)姿態(tài)攝動

用Tio(t)表示慣性系到軌道坐標(biāo)系的矢量轉(zhuǎn)換矩陣,在軌道坐標(biāo)系內(nèi),則纜繩的長度矢量在軌道,面內(nèi)分量可表示為:

電動繩系所受安培力可進一步在軌道坐標(biāo)系內(nèi)表示為:

該矩陣是關(guān)于時間的函數(shù),可根據(jù)電動繩系的軌道運動進行計算獲得。則安培力在軌道平面的分量可在軌道坐標(biāo)系內(nèi)表示為:

則安培力在軌道面內(nèi)的攝動力矩可表示為:

結(jié)合式 (15)和式 (18),則電動繩系軌道面內(nèi)姿態(tài)攝動方程可以表示為:

觀察式 (35)可以看到,當(dāng)電動繩系完成釋放且達到平衡之后,有:

整理之后則有:

式 (39)說明,在無控的條件下,電動繩系和當(dāng)?shù)卮咕€方向存在夾角,且該夾角是動態(tài)變化的。

5 考慮纜繩姿態(tài)的離軌仿真驗證

本部分結(jié)合火箭末級離軌任務(wù),建立了電動繩系離軌過程的數(shù)值仿真模型,仿真模型將系統(tǒng)軌道動力學(xué)、地磁場模型和繩系動力學(xué)三者結(jié)合到一起,在綜合考慮纜繩的面內(nèi)姿態(tài)擺動的條件下,驗證電動繩系的離軌效果。需要說明的是,本仿真模型當(dāng)中暫未考慮地球形狀攝動、大氣阻力攝動等因素,以第十一代地磁場模型為參照,建立地磁傾斜偶極子數(shù)值模型。電動繩系的參數(shù)如表1所示。

表1 EDT仿真參數(shù)Tab.1 Parameters for EDT simulations

假設(shè)繩系以初始速度1m/s,在600km高圓軌道,沿軌道半徑方向釋放過程面內(nèi)構(gòu)型變化情況,釋放過程中電動繩系未加電,因此無安培力和力矩作用。根據(jù)式 (18),由于繩系長度發(fā)生變化,因此在釋放過程,面內(nèi)擺角開始逐漸增大,當(dāng)完成釋放之后,電動繩系仍然在面內(nèi)作周期性擺動,這充分說明繩系在釋放過程必然會偏離軌道半徑方向,因此在分析其離軌過程時,假設(shè)繩系始終和當(dāng)?shù)卮咕€方向重合必然導(dǎo)致明顯的分析誤差。為說明考慮電動繩系面內(nèi)姿態(tài)擺動的必要性,本文分別在考慮和不考慮面內(nèi)姿態(tài)變化的情況下,對其離軌效果做了仿真對比,仿真周期為1天。

圖3是1天之內(nèi)電動繩系所處地磁場強度的變化情況。圖4是1天內(nèi)電動繩系軌道高度變化情況。仿真結(jié)果說明,兩種情況下1天降軌高度相差大約180m左右。圖5是電動繩系從不同初始軌道高度降低至400km所需要的離軌時間,其中 ‘?’和 ‘×’分別表示考慮和不考慮面內(nèi)擺動的結(jié)果,可以看到,兩者離軌時間最大相差近5天。當(dāng)初始軌道高度增加時,離軌時間偏差也會進一步增大。

圖3 EDT 24h內(nèi)所處地磁場矢量的變化Fig.3 Earth magnetic field vector for EDT during 24 hours

上述兩種情況之所以降軌高度不同,是因為兩種情況下電動繩系生成的電壓不同,進而導(dǎo)致安培力不同,圖6和圖7分別是兩種情況下電壓和安培力加速度變化情況,可以看到,兩種不同的仿真條件對應(yīng)的仿真結(jié)果具有明顯差別。

圖4 EDT 24h內(nèi)軌道高度變化Fig.4 Orbit altitude change during 24 hours for EDT

圖5 不同初始軌道的離軌時間變化Fig.5 Deorbit time for different initial orbits

圖6 EDT 24h內(nèi)電壓變化情況Fig.6 EDT Voltage history during 24 hours

圖7 EDT 24h安培力加速度變化Fig.7 EDT acceleration by Ampere force during 24 hours

圖8和圖9是在600km高,軌道傾角為0°和45°情況下,2天內(nèi)電動繩系面內(nèi)擺角以及面內(nèi)擺角加速度大小的變化情況,其中面內(nèi)擺角加速度的大小說明了電動繩系所受安培力矩的大小。圖9和圖10仿真結(jié)果說明,電動繩系在軌運行過程當(dāng)中,在重力梯度力矩和安培力矩綜合作用下,面內(nèi)擺角和所受安培力矩都是動態(tài)變化的。

圖8 i=0°時面內(nèi)振蕩情況Fig.8 In-plane libration for i=0°

圖9 i=45°時面內(nèi)振蕩情況Fig.9 In-plane libration for i=45°

6 結(jié) 論

采用拉格朗日方法建立了電動繩系離軌過程姿態(tài)動力學(xué)模型,描述了其在軌運行過程當(dāng)中的面內(nèi)外姿態(tài)變化情況,分析了電動繩系動力學(xué)模型實現(xiàn)面內(nèi)外姿態(tài)運動解耦的條件,隨后結(jié)合軌道動力學(xué)模型第十一代地磁場模型,在考慮重力梯度力矩和安培力矩作用的條件下,從理論上綜合分析了電動繩系面內(nèi)運動特性。本文最后建立了電動繩系離軌數(shù)值仿真模型,分別在考慮和不考慮面內(nèi)擺角條件下,對電壓變化、降軌高度和面內(nèi)擺動等進行了對比分析,仿真結(jié)果說明,在考慮面內(nèi)擺角的條件下,其離軌過程分析結(jié)果更為合理。