采用虛擬引力場(chǎng)的太陽(yáng)帆/電推進(jìn)混合推力機(jī)動(dòng)軌道設(shè)計(jì)

孫沖, 袁建平, 方群, 崔堯, 汪明曉, 劉易斯

(1．西北工業(yè)大學(xué)航天飛行動(dòng)力學(xué)技術(shù)國(guó)家級(jí)重點(diǎn)實(shí)驗(yàn)室，陜西西安710072;2．西北工業(yè)大學(xué)深圳研究生院，深圳518057;3．上海宇航系統(tǒng)工程研究所，上海201109; 4．國(guó)家航天局空間碎片監(jiān)測(cè)與應(yīng)用中心，北京100012)

連續(xù)小推力軌道設(shè)計(jì)在空間機(jī)動(dòng)軌道設(shè)計(jì)方面有廣泛的應(yīng)用[1-2]。與傳統(tǒng)脈沖推力或大推力相比，連續(xù)小推力軌道機(jī)動(dòng)容易控制，機(jī)動(dòng)能力強(qiáng)；另外其比沖小，實(shí)現(xiàn)軌道機(jī)動(dòng)所需要的燃料很少，因而成為研究熱點(diǎn)。文獻(xiàn)[3]中，Mcinnes首次提出了使用太陽(yáng)帆來(lái)改變引力場(chǎng)的大小，設(shè)計(jì)了懸浮軌道，并推導(dǎo)了實(shí)現(xiàn)懸浮軌道所需推力的解析解。Malcolm采用太陽(yáng)帆產(chǎn)生切向固定推力，設(shè)計(jì)了從地球軌道到太陽(yáng)極軌的轉(zhuǎn)移軌道，解決了軌道優(yōu)化問(wèn)題[4]。然而，由于單獨(dú)采用電推進(jìn)或者太陽(yáng)帆產(chǎn)生的推力較小，軌道轉(zhuǎn)移常常需要很長(zhǎng)時(shí)間。針對(duì)這個(gè)問(wèn)題，Mengali等在文獻(xiàn)[5]中提出了太陽(yáng)帆和太陽(yáng)帆/電推進(jìn)結(jié)合的混合推進(jìn)方法，并采用間接優(yōu)化方法進(jìn)行了軌道優(yōu)化。在其文章中，太陽(yáng)帆產(chǎn)生徑向推力，用于改變引力場(chǎng)的引力常數(shù)，而與此同時(shí)電推進(jìn)產(chǎn)生固定切向推力。在假定太陽(yáng)帆能夠產(chǎn)生50%的中心引力(μ2=0.5μ1)的條件下，混合推力能夠減少燃料的消耗并且縮短空間機(jī)動(dòng)任務(wù)時(shí)間。然而該方法是作為一種間接優(yōu)化方法，并沒(méi)有考慮推力約束問(wèn)題，很難應(yīng)用于實(shí)際空間機(jī)動(dòng)軌道設(shè)計(jì)中。

本文針對(duì)連續(xù)小推力軌道設(shè)計(jì)問(wèn)題，提出了一種基于虛擬引力場(chǎng)的混合連續(xù)小推力軌道設(shè)計(jì)方法[6-7]。該方法將虛擬引力場(chǎng)與平均算法結(jié)合，通過(guò)混合連續(xù)小推力形成虛擬引力場(chǎng)，進(jìn)而采用平均算法求解在虛擬引力場(chǎng)下固定大小切向力作用下的機(jī)動(dòng)軌道解析解。進(jìn)一步，本文采用四階龍格庫(kù)塔積分方法分析解析解的精度，并提出了誤差因子參數(shù)，從而保證解析解的精度。本文采用虛擬引力場(chǎng)方法，將航天器機(jī)動(dòng)軌道設(shè)計(jì)優(yōu)化問(wèn)題轉(zhuǎn)化成參數(shù)優(yōu)化問(wèn)題。該方法主要分為3步：①設(shè)計(jì)并優(yōu)化滿足任務(wù)約束(包括邊界約束和推力約束)的虛擬引力場(chǎng)，以及固定切向推力的大小；②求解實(shí)現(xiàn)虛擬引力場(chǎng)的推力大小；③考慮太陽(yáng)帆的動(dòng)力學(xué)模型，求解采用太陽(yáng)帆來(lái)實(shí)現(xiàn)虛擬引力場(chǎng)的控制律。

1 模 型

1.1 虛擬中心引力場(chǎng)的定義

虛擬中心引力場(chǎng)是指引力場(chǎng)中心不在地心上，且其引力常數(shù)不等于地心引力常數(shù)的中心引力場(chǎng)。該引力場(chǎng)可用參數(shù)(r0,μ2)來(lái)定義,其中r0為虛擬引力場(chǎng)中心與地心相對(duì)位置矢量,稱為虛擬引力場(chǎng)位置參數(shù);μ2為虛擬引力場(chǎng)引力常數(shù)[6-7]。

1.2 虛擬中心引力場(chǎng)下虛擬圓錐曲線軌道

根據(jù)虛擬中心引力場(chǎng)參數(shù)r0維數(shù)的不同,虛擬中心引力場(chǎng)可以分為2種,如圖1所示;

圖1 虛擬中心引力場(chǎng)示意圖

設(shè)定地心引力場(chǎng)引力常數(shù)為μ1,航天器在地心慣性坐標(biāo)系下位置矢量為r1,假定在航天器上施加推力加速度為at,則有:

(1)

假設(shè)地心慣性坐標(biāo)系為OXYZ,虛擬中心引力場(chǎng)坐標(biāo)系為O′X′Y′Z。推力與地心引力場(chǎng)的合力形成的虛擬中心引力場(chǎng)參數(shù)為(r0,μ2),則在虛擬中心引力場(chǎng)坐標(biāo)系下有:

(2)

(3)

由(2)式、(3)式可知,在虛擬中心引力場(chǎng)坐標(biāo)系下,航天器運(yùn)動(dòng)方程可表示為如下形式:

(4)

實(shí)現(xiàn)虛擬中心引力場(chǎng)所需推力為

(5)

將(3)式帶入(5)式, 可得

(6)

由上述分析可知,通過(guò)調(diào)整推力,可形成一種虛擬中心引力場(chǎng)。在該虛擬中心引力場(chǎng)中,航天器運(yùn)動(dòng)軌跡為圓錐曲線。若該虛擬圓錐曲線軌道滿足軌道轉(zhuǎn)移約束,則可實(shí)現(xiàn)軌道機(jī)動(dòng)。

1.3 固定切向力作用下軌道攝動(dòng)方程

在切向推力作用下,平面軌道高斯攝動(dòng)方程為

(7)

(8)

徑向推力以及周向推力和切向推力之間的關(guān)系為

(9)

而真近點(diǎn)角和偏近點(diǎn)角的關(guān)系為

(10)

因此,將公式(7)中自變量轉(zhuǎn)換為偏近點(diǎn)角,則軌道半長(zhǎng)軸和偏心率的變化[9]為

(11)

式中，E為偏近點(diǎn)角,ft為連續(xù)切向推力。在給定軌道交會(huì)的初始點(diǎn)A和終點(diǎn)B的條件下,采用混合推力方法設(shè)計(jì)機(jī)動(dòng)軌道就是求解引力場(chǎng)參數(shù)(μ2,rvg)以及切向推力ft。若選擇合適參數(shù)X=(μ2,rvg,fsep)時(shí),在虛擬引力場(chǎng)中和切向推力的同時(shí)作用下,航天器在虛擬引力場(chǎng)可從A點(diǎn)轉(zhuǎn)移到目標(biāo)點(diǎn)B點(diǎn),從而滿足軌道交會(huì)的邊界約束。這就是采用基于虛擬引力場(chǎng)方法的混合連續(xù)小推力軌道設(shè)計(jì)的基本思路。

采用虛擬引力場(chǎng)的混合連續(xù)小推力軌道設(shè)計(jì)方法來(lái)設(shè)計(jì)軌道,實(shí)際上是設(shè)計(jì)在虛擬引力場(chǎng)中的常值切向推力轉(zhuǎn)移軌道。轉(zhuǎn)移軌道需要滿足兩點(diǎn)邊界值約束,以及推力約束。當(dāng)初始軌道的偏心率很小,推力較小的條件下[9],在虛擬引力場(chǎng)中的切向推力軌道近似解析解為

(12)

當(dāng)初始點(diǎn)和目標(biāo)點(diǎn)已知時(shí),上式中的參數(shù)E已知。邊界約束可以轉(zhuǎn)化為

(13)

由此可見(jiàn),采用虛擬引力場(chǎng)下的連續(xù)推力軌道設(shè)計(jì)方法,可將所有滿足邊界約束的軌道用一組參數(shù)X=(μvg,rvg,fsep)表示,即轉(zhuǎn)移軌道參數(shù)化。這樣轉(zhuǎn)移軌道設(shè)計(jì)與優(yōu)化問(wèn)題就可以轉(zhuǎn)化成為參數(shù)設(shè)計(jì)與優(yōu)化問(wèn)題。

由上分析可知,采用虛擬引力場(chǎng)下的混合連續(xù)小推力軌道設(shè)計(jì)方法,需要給航天器施加推力,一部分推力用于形成虛擬引力場(chǎng),另一部分推力用于形成虛擬引力場(chǎng)中的切向力。具體推力分配:將形成虛擬引力場(chǎng)所需要的連續(xù)推力與在切向力的合力分解為徑向力和切向力,然后采用太陽(yáng)帆和電推進(jìn)分別產(chǎn)生分解的徑向力和切向力。由公式(6)可知,實(shí)現(xiàn)虛擬引力場(chǎng)所需要的推力是有虛擬引力場(chǎng)參數(shù)決定的如下式所示,

Trvg=F(μvg,rvg,fsep)

1.4 太陽(yáng)帆/電推進(jìn)混合最近最優(yōu)軌道設(shè)計(jì)

1.4.1 太陽(yáng)帆推進(jìn)模型

由公式(5)可知,實(shí)現(xiàn)虛擬引力場(chǎng)需要的力隨著軌道半徑的變化而變化。因而采用太陽(yáng)帆來(lái)實(shí)現(xiàn)虛擬引力場(chǎng),其推力大小和方向也都是實(shí)時(shí)變化的。在此采用太陽(yáng)帆推進(jìn)的理想模型[4]

Fsail=2PAcos2αn

(14)

式中，P=(S0/c)(r0/r)2,為太陽(yáng)光壓率。

采用理想推進(jìn)模型,太陽(yáng)帆產(chǎn)生的推力表示如下

(15)

式中，a1,a2,a3均為常數(shù)。而太陽(yáng)帆產(chǎn)生的力在軌道的徑向和切向的分量為

(16)

由此可知太陽(yáng)帆推進(jìn)可以實(shí)現(xiàn)部分徑向和周向的正推力,而剩余部分需要用電推進(jìn)來(lái)實(shí)現(xiàn)。

1.4.2 太陽(yáng)帆電推進(jìn)模型

本文選擇文獻(xiàn)[5]中的電推進(jìn)器,其產(chǎn)生的小推力和電推進(jìn)器功率參數(shù)相關(guān),而需要的燃料消耗也由點(diǎn)推進(jìn)器的功率決定。點(diǎn)推進(jìn)器的功率是可變的,如公式(17)所示

(17)

則其產(chǎn)生的推力加速度在徑向周向分量為:

(18)

1.4.3 推力分配

采用虛擬引力場(chǎng)混合小推力軌道設(shè)計(jì)方法設(shè)計(jì)軌道,需要推力來(lái)形成虛擬引力場(chǎng)以及在虛擬引力場(chǎng)中的切向力。本文采用太陽(yáng)帆/電推進(jìn)混合推力來(lái)實(shí)現(xiàn)虛擬引力場(chǎng),并提供在虛擬引力場(chǎng)中的常值推力。由于形成虛擬引力場(chǎng)所需的徑向力和周向力是固定的,而太陽(yáng)帆產(chǎn)生的推力在徑向和周向是耦合的,因此僅僅采用太陽(yáng)帆無(wú)法直接提供形成虛擬引力場(chǎng)需要的推力。在此,本文在選擇太陽(yáng)帆的同時(shí),另外選擇安裝2個(gè)電推進(jìn)器,分別用于產(chǎn)生徑向和周向的推力,使其與太陽(yáng)帆光壓推進(jìn)產(chǎn)生的徑向力和切向力相互匹配,共同作用產(chǎn)生形成引力勢(shì)場(chǎng),并提供常值推力。在此假定形成引力勢(shì)場(chǎng)所需要的徑向和周向的推力加速度為T(mén)r,Tc,而太陽(yáng)帆產(chǎn)生的徑向和周向的推力加速度為assr,assc。電推進(jìn)產(chǎn)生的徑向和周向的推力加速度為asepr,asepc,則有如下公式

(19)

由太陽(yáng)帆的數(shù)學(xué)模型可知,其產(chǎn)生的徑向力和軸向力與太陽(yáng)帆的位置和太陽(yáng)帆姿態(tài)有關(guān)。采用虛擬引力場(chǎng)和常值推力結(jié)合的方法,其太陽(yáng)帆的在每一時(shí)刻的位置是一定的,因此太陽(yáng)帆產(chǎn)生的推力,實(shí)際上是由太陽(yáng)帆的姿態(tài)角來(lái)確定的。為了能夠最大程度利用太陽(yáng)帆產(chǎn)生的推力,從而節(jié)省燃料,需要對(duì)其姿態(tài)角進(jìn)行優(yōu)化,在此假定太陽(yáng)帆的姿態(tài)角α和軌道機(jī)動(dòng)的轉(zhuǎn)移角θ有如下的關(guān)系:

α=a2·θ2+a1·θ+a0

(20)

式中，a0,a1,a2為待優(yōu)化的系數(shù)項(xiàng)。

當(dāng)給定虛擬引力場(chǎng)后,軌道轉(zhuǎn)移過(guò)程中需要電推進(jìn)器產(chǎn)生推力的大小為:

(21)

對(duì)比本文與文獻(xiàn)[5]可知,本文的推力分配與其不同。文獻(xiàn)[5]中軌道機(jī)動(dòng)的徑向推力完全由太陽(yáng)帆提供,而切向力由電推進(jìn)提供。本文采用虛擬引力場(chǎng)方法設(shè)計(jì)機(jī)動(dòng)軌道,太陽(yáng)帆推力一部分提供徑向力一部分提供切向力,這樣可以減少電推進(jìn)所需提供的推力。因而本文的推力分配方法比較靈活,可節(jié)省軌道機(jī)動(dòng)所需燃料。

1.4.4 混合小推力軌道設(shè)計(jì)與優(yōu)化

由上述分析可知,采用基于虛擬引力場(chǎng)的混合連續(xù)小推力軌道設(shè)計(jì)方法,可求出大量的滿足邊界約束的轉(zhuǎn)移軌道,如圖3所示。假定虛擬引力場(chǎng)參數(shù)以及切向力參數(shù)的參數(shù)范圍為

圖2 混合小推力軌道交會(huì)

圖3 混合小推力轉(zhuǎn)移軌道優(yōu)化示意圖

則軌道優(yōu)化問(wèn)題即可轉(zhuǎn)化為參數(shù)優(yōu)化問(wèn)題,優(yōu)化的參數(shù)為虛擬引力場(chǎng)參數(shù)和切向推力加速度,以及太陽(yáng)帆光壓推進(jìn)的推力系數(shù),而優(yōu)化的目標(biāo)函數(shù)可以為軌道轉(zhuǎn)移過(guò)程中,電推進(jìn)器消耗燃料質(zhì)量:

J=min(mfuel)

(22)

式中，mfuel為任務(wù)結(jié)束時(shí)航天器的質(zhì)量。

或是轉(zhuǎn)移軌道的任務(wù)時(shí)間:

J=min(Tf)

(23)

或是同時(shí)考慮燃料消耗和飛行時(shí)間的權(quán)重函數(shù)

J=min(k1·mfuel+k2·Tf)

(24)

式中，k1,k2為權(quán)重,k1+k2=1。

結(jié)合上述分析可知,整體軌道設(shè)計(jì)流程包括以下幾個(gè)步驟:

a) 給出初始軌道和目標(biāo)軌道;

b) 給出滿足推力約束的虛擬引力場(chǎng)參數(shù)范圍和切向推力范圍;

c) 在b)的范圍中,選擇虛擬引力場(chǎng)參數(shù)。

d) 將初始軌道狀態(tài)(ri,vi)轉(zhuǎn)化為虛擬引力

場(chǎng)中的軌道狀態(tài)(rvg,vvg);

e) 優(yōu)化切向力參數(shù),求解滿足邊界約束條件的轉(zhuǎn)移軌道;

f) 重復(fù)a)～e),得到大量的滿足約束的轉(zhuǎn)移軌道;

g) 采用Matlab軟件中的Fmincon函數(shù)優(yōu)化轉(zhuǎn)移軌道。

2 誤差分析

由第1節(jié)分析可知,采用虛擬引力場(chǎng)下混合連續(xù)小推力軌道設(shè)計(jì)方法,可求解出滿足推力約束和邊界值約束的解析解。然而在此需強(qiáng)調(diào),公式(10)中的解析解是存在誤差的。在此為了保證該解析方法的精度,本節(jié)采用四階龍格庫(kù)塔數(shù)值積分方法求解了軌道積分?jǐn)?shù)值解,求出了解析解的誤差。在此基礎(chǔ)上,分析誤差產(chǎn)生的原因并給出解析解的應(yīng)用范圍。具體的仿真條件見(jiàn)表1。

在此,將解析解的誤差分為位置誤差和速度誤差,并分別用Δεr,Δεv表示

Δεr=ranslytical-rnumerical=f(x)

Δεv=vanalytical-vnumerical=g(x)

對(duì)比圖4圖5,以及公式(10)、公式(11)可知,解析解與數(shù)值解的誤差和飛行時(shí)間、引力場(chǎng)參數(shù)、以及軌道參數(shù)有關(guān)。切向推力越大、飛行時(shí)間越長(zhǎng),該解析解的誤差會(huì)越大;初始軌道(虛擬引力場(chǎng)中)的偏心率越大、半長(zhǎng)軸越大其誤差越大。這一結(jié)論也與文獻(xiàn)[1]得出的結(jié)論一致。為描述解析解與數(shù)值解的誤差,在此定義一個(gè)誤差因子參數(shù),如(25)式所示

(25)

對(duì)比公式(11)及(12)可知,解析解的誤差因子會(huì)隨著切向推力、軌道半長(zhǎng)軸、軌道轉(zhuǎn)移時(shí)間(偏近點(diǎn)角)的增大而增大,當(dāng)(e<0.7)時(shí),偏心率越大,誤差參數(shù)越大。如圖5所示,采用解析方法可以得到大量的轉(zhuǎn)移軌道,同時(shí)采用四階龍格庫(kù)塔數(shù)值方法,就可以求解并分析這些解析解的速度誤差、位置誤差以及誤差因子之間的關(guān)系。

以地球-火星轉(zhuǎn)移軌道為例,采用虛擬引力場(chǎng)方法可以得到大量的轉(zhuǎn)移軌道,在其中任意選擇n條轉(zhuǎn)移軌道(n=157)。采用數(shù)值積分方法,可以得出這些解析解的速度誤差和位置誤差與誤差因子的關(guān)系,如圖6所示。

圖6 速度、位置誤差隨誤差因子的變化

圖6中,藍(lán)色的線條為位置誤差;紅色的線條為速度誤差;綠色的線條為誤差因子。由圖6可知,速度誤差、位置誤差和誤差因子是同步的,當(dāng)誤差因子很小時(shí),其對(duì)應(yīng)的速度誤差和位置誤差也很小;相反,當(dāng)誤差因子比較大時(shí),其速度誤差和位置誤差都很大。由圖6可知，速度誤差和位置誤差都隨著誤差因子的增大而增大。當(dāng)保證誤差因子小于0.002 5時(shí),速度誤差小于0.001%,而位置誤差小于0.002%。由以上分析可知,在采用虛擬引力場(chǎng)常值切向力解析方法設(shè)計(jì)轉(zhuǎn)移軌道時(shí),只需要保證誤差因子k

3 數(shù)值仿真

表1 地球火星軌道轉(zhuǎn)移虛擬引力場(chǎng)參數(shù)和常值推力

表2 太陽(yáng)帆姿態(tài)角系數(shù)

圖7 虛擬引力場(chǎng)混合推力方法的轉(zhuǎn)移軌道

圖8 切向推力方法的轉(zhuǎn)移軌道(相同時(shí)間)

圖7為采用虛擬引力場(chǎng)混合推力方法設(shè)計(jì)的轉(zhuǎn)移軌道；圖8為在相同轉(zhuǎn)移時(shí)間條件下，采用常值切向力方法設(shè)計(jì)的轉(zhuǎn)移軌道。圖9表述轉(zhuǎn)移軌道所需的周向推力和徑向推力以及太陽(yáng)帆提供的周向推力和徑向推力。圖10是轉(zhuǎn)移軌道太陽(yáng)帆姿態(tài)角。圖11是為相同情況下(切向推力大小相等)得到的轉(zhuǎn)移軌道。圖12為實(shí)現(xiàn)虛擬引力場(chǎng)所需要的推力加速度。圖13和圖14表示推力秒流量和推力方向角。由表5可知，相比于常值切向推力，虛擬引力場(chǎng)混合推力方法能夠很大程度上節(jié)省燃料，縮短時(shí)間。假定軌道轉(zhuǎn)移時(shí)間相等，虛擬引力場(chǎng)混合推力方法所需要的燃料僅為常值推力方法的 。采用虛擬引力場(chǎng)混合連續(xù)小推力方法，所需要的電推進(jìn)推力大大減小了，因而虛擬引力場(chǎng)混合推力方法能夠擴(kuò)大小推力軌道設(shè)計(jì)方法的應(yīng)用范圍。另一方面，在同樣大小推力條件下，采用常值推力方法，需要的軌道轉(zhuǎn)移時(shí)間和燃料消耗是相當(dāng)大的。

圖9 需要的推力和太陽(yáng)帆提供的推力 圖10 太陽(yáng)帆姿態(tài)角圖11 切向推力轉(zhuǎn)移軌道(與電推進(jìn)相同推力大小)

圖12 電推進(jìn)產(chǎn)生的推力加速度 圖13 電推進(jìn)需要的秒流量圖14 電推進(jìn)方向角

半長(zhǎng)軸(AU)偏心率真近點(diǎn)角/rad真近點(diǎn)角/rad 初始值1000 目標(biāo)值1.501.484(第一組)3.108(第二組)

表4 虛擬引力場(chǎng)參數(shù)和切向推力參數(shù)

表5 混合推力方法和切向推力力方法對(duì)比

4 結(jié) 論

本文提出了一種基于虛擬引力場(chǎng)的混合小推力設(shè)計(jì)方法。該方法能夠以解析的形式求出大量的滿足軌道邊界約束的轉(zhuǎn)移軌道，以及實(shí)現(xiàn)軌道轉(zhuǎn)移所需的推力。同時(shí)，在給定目標(biāo)函數(shù)的條件下，可對(duì)滿足推力約束的轉(zhuǎn)移軌道進(jìn)行優(yōu)化。與其它的小推力軌道設(shè)計(jì)方法相比，該方法有以下的特點(diǎn)：

1) 該軌道設(shè)計(jì)方法是一種解析方法，在保證精度的條件下，能夠減小計(jì)算量，因而該方法可以為數(shù)值優(yōu)化方法(直接法或間接法)提供軌道優(yōu)化初值；

2) 與常值切向推力相比，該方法能夠減小實(shí)現(xiàn)軌道轉(zhuǎn)移所需的推力，擴(kuò)大了小推力軌道應(yīng)用的范圍；

3)與常值切向推力相比，該方法不僅能夠減小軌道轉(zhuǎn)移所消耗的能量，而且能夠減小軌道轉(zhuǎn)移時(shí)間。

綜上所述，相比于傳統(tǒng)的小推力設(shè)計(jì)方法，虛擬引力場(chǎng)混合小推力方法，不僅能夠節(jié)省燃料消耗，而且能減小任務(wù)時(shí)間，并擴(kuò)大了小推力方法的應(yīng)用范圍。為小推力軌道設(shè)計(jì)與優(yōu)化領(lǐng)域提供了新的思路。