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

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

0 引言

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

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

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

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

αns；

(1)

(2)

(3)

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

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

a=aa+ar+ath. (4)

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

圖1 圓型限制性三體問(wèn)題Fig.1 Circular restricted three-body problem

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

(5)

(6)

(7)

其中Ω表示勢(shì)函數(shù)，見(jiàn)文獻(xiàn)[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)

運(yùn)用文獻(xiàn)[14]中提到的方法，即把原點(diǎn)移動(dòng)到小天體并進(jìn)行坐標(biāo)變換：

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

(15)

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

(16)

(17)

(18)

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

(19)

(20)

3 共線平衡點(diǎn)

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

WX+ax=0；

(21)

WY+ay=0；

(22)

WZ+az=0.

(23)

經(jīng)典的希爾型三體問(wèn)題只有兩個(gè)對(duì)稱的共線平衡點(diǎn)，筆者主要研究希爾型三體系統(tǒng)的共線平衡點(diǎn)和其附近鄰域的軌道運(yùn)動(dòng)情況.假設(shè)該系統(tǒng)的L1、L2兩共線平衡點(diǎn)位于X軸上，用Re=(Xe,0,0)T表示共線平衡點(diǎn).很顯然共線平衡點(diǎn)滿足式(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)共線平衡點(diǎn)的位置由p、a0和A1共同決定.運(yùn)用類似文獻(xiàn)[15]的數(shù)值計(jì)算方法，在3個(gè)參數(shù)p、a0、A1當(dāng)中，固定一個(gè)參數(shù)不變，然后觀察其他兩個(gè)參數(shù)變化對(duì)共線平衡點(diǎn)位置的影響.圖2表示的是，當(dāng)太陽(yáng)帆特征加速度為一常值時(shí)，參數(shù)p和A1的改變對(duì)L1、L2位置的影響.經(jīng)典的希爾型限制性三體問(wèn)題的共線平衡點(diǎn)位置為X=±3-1/3≈±0.693 361 ，然而含扁率的希爾系統(tǒng)的L1的橫坐標(biāo)最小值為-0.693 459 ，L2的最大值為0.693 264.然后，按照同樣的方法，固定一個(gè)參數(shù)，設(shè)定扁率A1為一個(gè)常數(shù)，A1=0.004 23，則共線平衡點(diǎn)位置變化顯示在圖3.

圖2 含扁率的希爾系統(tǒng)共線平衡點(diǎn)位置變化圖，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)共線平衡點(diǎn)位置變化圖，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)線性化

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

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

(27)

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

ξ+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ì)函數(shù)的二階偏導(dǎo)數(shù)在共線平衡點(diǎn)的取值，可以把變分方程寫(xiě)成狀態(tài)空間的矩陣方程：

(35)

(36)

5 仿真計(jì)算

5.1 追蹤參考軌道

求解共線平衡點(diǎn)附近的太陽(yáng)帆懸浮軌道，運(yùn)用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)對(duì)太陽(yáng)帆進(jìn)行主動(dòng)控制，從而控制太陽(yáng)帆追蹤給定的參考軌道Xref=(ξref,ηref,ζref)T，主動(dòng)控制可以使太陽(yáng)帆軌道趨于漸近穩(wěn)定.考慮太陽(yáng)帆軌道與參考軌道之間的誤差ΔX=X-Xref，然后應(yīng)用線性反饋控制，Δu=u-uref=-K(X-Xref)，使誤差ΔX滿足性能指標(biāo)函數(shù)J取最小值：

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

(43)

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

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

(44)

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

(45)

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

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

5.2 仿真計(jì)算

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

圖4 太陽(yáng)帆懸浮軌道Fig.4 Solar sail displaced orbit

太陽(yáng)光壓加速度沿三軸的分量隨時(shí)間的變化情況見(jiàn)圖6.可以看出,太陽(yáng)光壓加速度在三軸的分量以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]小幅度變化，在太陽(yáng)帆運(yùn)行大概4.7個(gè)單位時(shí)間后，沿該軸的加速度趨于平穩(wěn)，穩(wěn)定在0.005 931附近.圖 7表示的是RCD控制參數(shù)p和太陽(yáng)帆姿態(tài)角α隨時(shí)間的變化圖.由圖可以看出，兩個(gè)變化的參數(shù)都是在開(kāi)始的一段時(shí)間大幅度變化，之后進(jìn)入平穩(wěn)期，參數(shù)p最終穩(wěn)定在0.202附近，而姿態(tài)角α最終穩(wěn)定在44.99°左右，這也表明兩個(gè)參數(shù)最終穩(wěn)定值取決于參考軌道的對(duì)應(yīng)參數(shù)值.

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

圖6 太陽(yáng)帆光壓加速度Fig.6 Solar radiation pressure acceleration

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

6 結(jié)論

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

參考文獻(xià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.