宋 明， 和興鎖， 閆業(yè)毫， 和東生

(西北工業(yè)大學(xué) 力學(xué)與土木建筑學(xué)院，陜西 西安 710029)

0 引言

太陽帆是利用巨大的、像鏡子一樣的輕質(zhì)帆面反射太陽光來產(chǎn)生推動力的.由日本航天局研制的“IKAROS”太陽帆成功發(fā)射后[1]，基于太陽帆的一系列空間任務(wù)被相繼提出，如：NanoSail-D, Cubesail, DeorbitSail, Sunjammer，LightSail等.毫無疑問，太陽帆將會成為本世紀的熱門研究課題.

限制性三體問題模型一直被視為基本的動力學(xué)模型[2]，而希爾型三體問題是限制性三體問題的一種近似[3].通常情況下人們把天體視為質(zhì)量均勻分布的理想球體，不考慮質(zhì)量為無窮小量的太陽帆對天體運動的影響.但是，宇宙中純粹意義的球體是不存在的，大部分星體都是扁的，也有一部分是不規(guī)則形狀的.Sharma[4]發(fā)現(xiàn)考慮大天體扁率的限制性三體系統(tǒng)存在周期性軌道.Singh[5]把諸如天體扁率、天體的輻射等多種擾動影響結(jié)合在一起，分析平衡點的非線性穩(wěn)定性.Douskos[6]對希爾型三體問題的平衡點及其穩(wěn)定性進行了系統(tǒng)研究.遺憾的是，以上文獻沒有涉及到天體扁率對太陽帆軌道的影響.因此，筆者主要研究天體扁率在希爾型三體系統(tǒng)中對太陽帆運動的影響，以及共線平衡點附近的懸浮軌道設(shè)計.

1 太陽帆力學(xué)模型

“IKAROS”太陽帆上的反射控制設(shè)備，簡稱RCD，是一種柔性多層液晶薄片通過改變其上的電壓來改變反射率的裝置[7].龔勝平等[8]提出了一種RCD，其中包含兩種模塊：鏡面反射模塊和吸收光線模塊，并把這種技術(shù)成功運用在研究地月系統(tǒng)平衡點太陽帆軌道運動和研究太陽帆在小行星附近的動力學(xué)[9-10].筆者將會利用這一技術(shù)并進行一些改進.改進型的RCD包含兩種模塊，一種是反射模塊，反射系數(shù)為ρr；另外一種是吸收模塊，吸收系數(shù)為ρa=1-ρr.和之前版本的RCD不同的是，筆者把散射和熱輻射一并考慮在內(nèi).假設(shè)RCD可以轉(zhuǎn)換成開啟和關(guān)閉兩種狀態(tài).太陽帆總的有效面積為Atot，設(shè)定RCD的面積為A0=Atot/8，吸收模塊的面積為Aa=pAtot,p∈[0,0.125]，其中p為比例系數(shù)，吸收模塊對應(yīng)的狀態(tài)為關(guān)閉；同樣反射模塊的面積為Ar=Atot-Aa，對應(yīng)開啟狀態(tài).于是，因吸收、反射和熱輻射作用而產(chǎn)生的太陽帆加速度可以表示為：

αns；

(1)

(2)

(3)

式中，P表示太陽光壓；太陽帆帆面法向量為n=(cosαcos (ωst),-cosαsin (ωst),sinα)T；ns=(cos (ωst),-sin (ωst),0)T表示太陽光線方向矢量；太陽帆質(zhì)量為m；α是帆面法向量與光線的夾角；ωs是光線在無量綱的相合坐標系下的角速度[11]，筆者視ωs為一常數(shù).

太陽帆力學(xué)模型的相關(guān)參數(shù)詳見表 1，詳見文獻[12]，太陽帆受到光壓加速度為：

a=aa+ar+ath. (4)

2 動力學(xué)系統(tǒng)

圖1 圓型限制性三體問題Fig.1 Circular restricted three-body problem

考慮大天體是扁球的情況下，太陽帆在該系統(tǒng)的運動微分方程為：

(5)

(6)

(7)

其中Ω表示勢函數(shù)，見文獻[6].

Ω=n2(x2+y2)/2+(1-μ)/r1+μ/r2+

(1-μ)A1/2/r13-3(1-μ)A1z2/2/r15.

(8)

ax=a0[f(p,α)+g(p,α)]cos (ωst)；

(9)

ay=-a0[f(p,α)+g(p,α)]sin (ωst)；

(10)

az=a0h(p,α)，

(11)

式中，f(·)、g(·)、h(·)是p和α的函數(shù).

f(p,α)=(1-p)(1.654 4cos2α+

0.041 7cosα-0.827 2)cosα;

(12)

g(p,α)=p(1-0.052 6cosα)cosα;

(13)

h(p,α)= (1.654 4cosα-1.654 4pcosα-

0.094 3p+0.041 7)sinαcosα.

(14)

運用文獻[14]中提到的方法，即把原點移動到小天體并進行坐標變換：

x=1-μ+μ1/3X;y=μ1/3Y;z=μ1/3Z.

(15)

令μ→0，則含扁率的希爾型限制性三體問題可以表示為：

(16)

(17)

(18)

式中:W為新系統(tǒng)的勢函數(shù).

(19)

(20)

3 共線平衡點

求解以下方程組可以得到系統(tǒng)的平衡點：

WX+ax=0；

(21)

WY+ay=0；

(22)

WZ+az=0.

(23)

經(jīng)典的希爾型三體問題只有兩個對稱的共線平衡點，筆者主要研究希爾型三體系統(tǒng)的共線平衡點和其附近鄰域的軌道運動情況.假設(shè)該系統(tǒng)的L1、L2兩共線平衡點位于X軸上，用Re=(Xe,0,0)T表示共線平衡點.很顯然共線平衡點滿足式(21)、(22)和(23)，整理可以得到：

3Xe+ 15A1Xe/2-Xe/Xe3+
a0[f(p,α)+g(p,α)]=0;

(24)

-a0[f(p,α)+g(p,α)]sin(ωst)=0;

(25)

a0h(p,α)=0.

(26)

α=0,ωst=2kπ(k∈Ζ)是式(25)和(26)成立的必要條件，把其代入式(24)，觀察式(24)，發(fā)現(xiàn)共線平衡點的位置由p、a0和A1共同決定.運用類似文獻[15]的數(shù)值計算方法，在3個參數(shù)p、a0、A1當中，固定一個參數(shù)不變，然后觀察其他兩個參數(shù)變化對共線平衡點位置的影響.圖2表示的是，當太陽帆特征加速度為一常值時，參數(shù)p和A1的改變對L1、L2位置的影響.經(jīng)典的希爾型限制性三體問題的共線平衡點位置為X=±3-1/3≈±0.693 361 ，然而含扁率的希爾系統(tǒng)的L1的橫坐標最小值為-0.693 459 ，L2的最大值為0.693 264.然后，按照同樣的方法，固定一個參數(shù)，設(shè)定扁率A1為一個常數(shù)，A1=0.004 23，則共線平衡點位置變化顯示在圖3.

圖2 含扁率的希爾系統(tǒng)共線平衡點位置變化圖，a0=0.001,p∈[0,0.125],A1∈[0,1] Fig.2 Positions of collinear equilibrium points in the Hill’s system with oblateness with variational p and A1 with fixed a0=0.001 ，p∈[0,0.125]，A1∈[0,1]

圖3 含扁率的希爾系統(tǒng)共線平衡點位置變化圖，A1=0.004 23,p∈[0,0.125],a0∈[0,1] Fig.3 Positions of collinear equilibrium points in the Hill’s system with oblateness with variational p and a0 with fixed A1=0.004 23, p∈[0,0.125]，a0∈[0,1]

4 系統(tǒng)線性化

為了更進一步研究太陽帆軌道，需要對系統(tǒng)微分方程進行線性化處理，即在平衡點處引入小的擾動，然后利用泰勒展開式，略去對系統(tǒng)影響很小的高階項，只保留線性項，建立變分方程.于是，在平衡點引入小擾動：

X=Xe+ξ,Y=η,Z=ζ.

(27)

然后把小擾動代入到系統(tǒng)方程，再進行線性化處理就可得到變分方程：

ξ+uξ；

(28)

(29)

(30)

其中，

WXX=3X2/R5-1/R3+3.0375；

(31)

WYY=3Y2/R5-1/R3；

(32)

WZZ=3Z2/R5-1/R3-1.0225；

(33)

u=[uξ,uη,uζ]T=[ax,ay,az]T.

(34)

此處的Wije(i,j=X,Y,Z)表示勢函數(shù)的二階偏導(dǎo)數(shù)在共線平衡點的取值，可以把變分方程寫成狀態(tài)空間的矩陣方程：

(35)

(36)

5 仿真計算

5.1 追蹤參考軌道

求解共線平衡點附近的太陽帆懸浮軌道，運用Simo[11]提出的方法，引入?yún)⒖架壍溃?/p>

ξref=ξ0cos (ωst)；

(37)

ηref=η0sin (ωst)；

(38)

ζref=ζ0.

(39)

參考軌道滿足系統(tǒng)變分方程，把參考軌道代入式(28)～式(30)，從而得到:

[f(p,α)+g(p,α)];

(40)

[f(p,α)+g(p,α)];

(41)

(42)

確定了參考軌道，就可以利用線性二次調(diào)節(jié)器(LQR)對太陽帆進行主動控制，從而控制太陽帆追蹤給定的參考軌道Xref=(ξref,ηref,ζref)T，主動控制可以使太陽帆軌道趨于漸近穩(wěn)定.考慮太陽帆軌道與參考軌道之間的誤差ΔX=X-Xref，然后應(yīng)用線性反饋控制，Δu=u-uref=-K(X-Xref)，使誤差ΔX滿足性能指標函數(shù)J取最小值：

ΔXTQΔX+ΔuTRΔu)dt.

(43)

其中,Q和R為系統(tǒng)的加權(quán)矩陣，是對稱正定或半正定的，可以自由選擇.通過求解代數(shù)黎卡提方程(Algebraic Riccati Equation)，見文獻[16]：

ATP+PA-PBR-1BTP+Q=0,

(44)

可以獲得收益矩陣K=R-1BTP，這樣就把非線性系統(tǒng)(35)轉(zhuǎn)化為穩(wěn)定的線性系統(tǒng):

(45)

判斷太陽帆軌道是否穩(wěn)定取決于矩陣A-BK特征值的實部是否是小于或等于0.表 2給出系統(tǒng)的相關(guān)參數(shù)和仿真的初始條件，這些參數(shù)的選取可根據(jù)具體的任務(wù)要求而確定，具有任意性.

表2 系統(tǒng)參數(shù)和仿真初始條件

5.2 仿真計算

通過仿真計算，可以求出理想太陽帆模型下共線平衡點L1的橫坐標位置-0.693 361 274 35.圖 4表示的是在主動控制下太陽帆周期懸浮軌道，可以看出太陽帆從擾動點出發(fā)，通過主動控制，逐漸接近參考軌道，最終幾乎與參考軌道重合，圖 5給出懸浮軌道與參考軌道的相對位置與速度誤差，最大位置誤差為2.947 68×10-5，最大的速度誤差為7.464 65×10-5，當軌道趨于漸近穩(wěn)定后，位置誤差保持在10-8左右，速度誤差保持在10-15左右.仿真結(jié)果表明通過追蹤參考軌道而設(shè)計出來的懸浮軌道可以達到漸近穩(wěn)定的狀態(tài).

圖4 太陽帆懸浮軌道Fig.4 Solar sail displaced orbit

太陽光壓加速度沿三軸的分量隨時間的變化情況見圖6.可以看出,太陽光壓加速度在三軸的分量以x軸變化幅度最大，y軸的變化次之，z軸的變化幅度最小，而且沿x,y兩軸的加速度分量變化具有周期性：x軸的加速度最大值為3.404 9×10-4，最小值為-3.404 8×10-4；y軸的加速度最大值為3.415 7×10-4，最小值為-6.401 1×10-4；z軸的光壓加速度大致在區(qū)間[0.005 920,0.005 958]小幅度變化，在太陽帆運行大概4.7個單位時間后，沿該軸的加速度趨于平穩(wěn)，穩(wěn)定在0.005 931附近.圖 7表示的是RCD控制參數(shù)p和太陽帆姿態(tài)角α隨時間的變化圖.由圖可以看出，兩個變化的參數(shù)都是在開始的一段時間大幅度變化，之后進入平穩(wěn)期，參數(shù)p最終穩(wěn)定在0.202附近，而姿態(tài)角α最終穩(wěn)定在44.99°左右，這也表明兩個參數(shù)最終穩(wěn)定值取決于參考軌道的對應(yīng)參數(shù)值.

圖5 太陽帆懸浮軌道和參考軌道間的相對位置誤差和相對速度誤差Fig.5 Relative position error and velocity error between solar sail displaced orbit and refer orbit

圖6 太陽帆光壓加速度Fig.6 Solar radiation pressure acceleration

圖7 RCD控制參數(shù)p和太陽帆姿態(tài)角α隨時間變化圖Fig.7 Time histories of control parameter of RCD and the pitch angle of the solar sail

6 結(jié)論

在考慮大天體含扁率的情況下，研究了太陽帆在希爾型限制性三體問題的懸浮軌道設(shè)計方案，并利用近似簡化，建立太陽帆在希爾型限制性三體問題系統(tǒng)中的動力學(xué)模型.研究發(fā)現(xiàn)，改變反射控制設(shè)備中吸收光線模塊和熱輻射模塊的面積大??；或者改變太陽帆特征加速度；或者改變大天體扁率都會引起系統(tǒng)的共線平衡點的位置發(fā)生變化.對系統(tǒng)進行了線性化處理，運用LQR對不穩(wěn)定的系統(tǒng)進行主動控制.通過調(diào)節(jié)太陽帆姿態(tài)角和改變吸收光線模塊的面積可以得到太陽帆懸浮軌道并且使其達到漸近穩(wěn)定.

參考文獻：

[1] TSUDA Y, MORI O, FUNASE R, et al. Flight status of IKAROS deep space solar sail demonstrator [J]. Acta astronaut, 2011, 69(9/10): 833-840.

[2] BAOYIN H, MCINNES C. Solar sail equilibria in the elliptical restricted three-body problem [J]. J Guid Control Dynam, 2006, 29(3): 538-543.

[3] SOSNITSKII S. On the hill stable motions in the three-body problem [J]. Adv Space Res, 2015, 56(5): 859-864.

[4] SHARMA R K. Periodic-orbits of the 3rd kind in the restricted 3-body problem with oblateness [J]. Astrophys & space scicence, 1990, 166(2): 211-218.

[5] SINGH J. Combined effects of oblateness and radiation on the nonlinear stability of L4 in the restricted three-body problem [J]. Astron J, 2009, 137(2): 3286-3292.

[6] DOUSKOS C N, MARKELLOS V V. Out-of-plane equilibrium points in the restricted three-body problem with oblateness [J]. Astronomy and astrophysics, 2006, 446(1): 357-360.

[7] TSUDA Y, MORI O, FUNASE R, et al. Achievement of IKAROS - Japanese deep space solar sail demonstration mission [J]. Acta astronaut, 2013, 82(2): 183-188.

[8] MU J S, GONG S P, LI J F. Reflectivity-controlled solar sail formation flying for magnetosphere mission [J]. Aerosp Sci Technol, 2013, 30(1): 339-348.

[9] GONG S P, LI J F, SIMO J. Orbital motions of a solar sail around the L2earth-moon libration point [J]. Journal of guidance, control, and dynamics, 2014, 37(4): 1349-1356.

[10] GONG S P, LI J F. Equilibria near asteroids for solar sails with reflection control devices [J]. Astrophys & space science, 2015, 355(2): 213-223.

[11] SIMO J, MCINNES C. Solar sail orbits at the Earth-Moon libration points [J]. Commun nonlinear scicence, 2009, 14(12): 4191-4196.

[12] MCINNES C R. Solar sailing: technology, dynamics and mission applications[M]. London: Springer, 1999.

[13] SHARMA R K, RAO P V S. A case of commensurability induced by oblateness [J]. Celestial mechanics and dynamical astronomy, 1978, 18(2): 185-194.

[14] MARKELLOS V V, ROY A E, PERDIOS E A, et al. A hill problem with oblate primaries and effect of oblateness on hill stability of orbits [J]. Astrophys & space science, 2001, 278(3): 295-304.

[15] TSIROGIANNIS G A, DOUSKOS C N, PERDIOS E A. Computation of the Liapunov orbits in the photogravitational RTBP with oblateness [J]. Astrophys & space science, 2006, 305(4): 389-398.

[16] WIE B. Space vehicle dynamics and control[J]. Aircraft engineering and oerospace technology, 1998,70(5): 2077-2078.