司震 錢霙婧 楊曉東 張偉

(北京工業(yè)大學(xué)材料與制造學(xué)部,機(jī)械結(jié)構(gòu)非線性振動(dòng)與強(qiáng)度北京市重點(diǎn)實(shí)驗(yàn)室,北京 100124)

引言

小行星形成于太陽(yáng)系早期,在幾十億年的演化進(jìn)程中,完整保留了太陽(yáng)系早期物質(zhì)等原始信息,不僅有助于人類探究太陽(yáng)系起源及生命起源等謎題,也便于人類采取有效的應(yīng)對(duì)措施來防御、減緩或消除碰撞威脅[1].如今,各國(guó)正逐步開展小行星的探測(cè)采樣任務(wù),JASA 的隼鳥號(hào)(Hayabusa)探測(cè)器于2005 年登陸小行星25143 Itokawa,成為人類第一個(gè)小行星采樣返回任務(wù)[2].自此以后,JASA 的隼鳥II 號(hào)(Hayabusa II) 探測(cè)器、NASA 的OSIRIS-Rex 探測(cè)器相繼對(duì)小行星162173 Ryugu、101955 Bennu 開展探測(cè)采樣任務(wù)[3-4].

目前,小行星的探測(cè)采樣任務(wù)主要采用懸停探測(cè)軌道,但由于探測(cè)器所處的動(dòng)力學(xué)環(huán)境的差異性,針對(duì)小行星探測(cè)任務(wù)類型或任務(wù)的不同階段,需要對(duì)各階段的動(dòng)力學(xué)進(jìn)行分析進(jìn)而設(shè)計(jì)相應(yīng)的懸停軌道[5-6].隼鳥號(hào)探測(cè)器在小行星Itokawa 附近執(zhí)行任務(wù)的大部分時(shí)間里都采用慣性系懸停方式,即探測(cè)器跟隨小行星進(jìn)行公轉(zhuǎn)運(yùn)動(dòng),在小行星軌道系中進(jìn)行位置保持.而在小行星表面附近進(jìn)行采樣時(shí),則采取本體系懸停操作,即相對(duì)于小行星表面某位置點(diǎn)進(jìn)行位置保持,以便于在采樣返回任務(wù)中進(jìn)行下降和上升機(jī)動(dòng)[7-9].

太陽(yáng)系中已知小行星自身的引力場(chǎng)普遍呈現(xiàn)十分復(fù)雜的不規(guī)則特性,其附近探測(cè)器的運(yùn)行軌道受不規(guī)則引力場(chǎng)影響[10-11].Scheeres 最早提出了考慮不規(guī)則引力場(chǎng)影響下本體系懸停軌道的受控運(yùn)動(dòng)[12].Sawai 考慮小行星引力場(chǎng)不規(guī)則特性,應(yīng)用球諧、橢球諧模型、多面體模型等研究了均勻自旋和非均勻自旋小行星附近探測(cè)器的本體系懸停軌道[13-14].Broschart 在不考慮第三體引力及外界干擾情況下,對(duì)探測(cè)器施加控制推力來抵消小行星引力及慣性力,研究了本體系懸停和慣性系懸停的穩(wěn)定區(qū)域[9,15].Zhang 等[16]參照限制性三體問題,考慮了小行星不規(guī)則引力對(duì)探測(cè)器軌道運(yùn)動(dòng)的影響,研究了翻滾小行星附近的慣性系懸?？刂茊栴}.Zeng 等[17]從非理想太陽(yáng)帆模型的角度,分析了小行星不規(guī)則引力分布對(duì)懸停軌道的影響.

由于小行星普遍呈現(xiàn)質(zhì)量小、引力場(chǎng)弱的特點(diǎn),研究表明對(duì)于在小行星附近懸停運(yùn)動(dòng)的航天器而言,太陽(yáng)輻射光壓以及太陽(yáng)引力攝動(dòng)不可忽略[18-21].在光壓攝動(dòng)方面,Morrow 等[22]利用太陽(yáng)光壓輻射設(shè)計(jì)并分析了小行星附近的太陽(yáng)帆懸停軌道.Williams等[23]在Morrow 等的基礎(chǔ)上,根據(jù)小行星的球形簡(jiǎn)化模型,研究了太陽(yáng)光壓輻射影響下慣性系和本體系懸停探測(cè)軌道設(shè)計(jì)問題.Giancotti 和Funase[24]對(duì)太陽(yáng)帆在特洛伊小行星附近的人工平衡點(diǎn)以及附近軌跡進(jìn)行了可行性研究.Feng 等[25]研究了非球面項(xiàng)C20和太陽(yáng)輻射壓力對(duì)小行星附近航天器軌道的影響.Xin 等[26-27]以橢球體小行星模型,考慮太陽(yáng)光壓影響,研究了小行星懸停軌道及其穩(wěn)定性.在太陽(yáng)引力攝動(dòng)方面,劉瑩瑩等[28]研究了非球形攝動(dòng)、太陽(yáng)光壓攝動(dòng)以及太陽(yáng)引力攝動(dòng)對(duì)繞飛小行星航天器軌道的影響.倪彥碩等[29]研究了考慮太陽(yáng)引力攝動(dòng)的小行星附近軌道動(dòng)力學(xué),結(jié)果表明太陽(yáng)引力對(duì)平衡點(diǎn)附近航天器運(yùn)動(dòng)的影響較大[29].

由于太陽(yáng)的周期運(yùn)動(dòng),光壓攝動(dòng)與引力攝動(dòng)在系統(tǒng)動(dòng)力學(xué)方程中均呈現(xiàn)為周期性激勵(lì)[30-33].根據(jù)文獻(xiàn)分析可知,光壓攝動(dòng)呈現(xiàn)外激勵(lì)的形式,即在系統(tǒng)動(dòng)力學(xué)中表現(xiàn)為非齊次項(xiàng).然而,引力攝動(dòng)則呈現(xiàn)參數(shù)激勵(lì)的形式,即在系統(tǒng)動(dòng)力學(xué)中表現(xiàn)為周期性剛度項(xiàng)系數(shù).通常,發(fā)生外激勵(lì)共振時(shí),外激頻率需要接近系統(tǒng)的固有頻率;而發(fā)生參激共振時(shí),參數(shù)激勵(lì)頻率需要接近系統(tǒng)的兩倍固有頻率或與兩倍固有頻率相關(guān)的加減組合關(guān)系.由非線性振動(dòng)理論可知,兩種形式的共振并不能同時(shí)發(fā)生[34-35].因此,本文著重研究小行星不規(guī)則引力場(chǎng)與太陽(yáng)引力攝動(dòng)下,系統(tǒng)參數(shù)激勵(lì)共振軌道設(shè)計(jì)問題,為小行星平衡點(diǎn)附近的懸停軌道設(shè)計(jì)提出新思路.

為了考慮不規(guī)則小行星引力場(chǎng)和太陽(yáng)引力,本文采用受攝粒桿模型,即將小行星視為由無(wú)質(zhì)量桿連接的質(zhì)量粒子組成[36-39],并加入引力攝動(dòng)天體部分.本文把太陽(yáng)引力視為受攝小行星系統(tǒng)的組成部分,而非單純的攝動(dòng)源,采用非線性振動(dòng)理論中的參激共振概念,推導(dǎo)了小行星平衡點(diǎn)附近穩(wěn)定的參數(shù)共振軌道.研究了受攝小行星系統(tǒng)的幅頻響應(yīng)曲線、力頻響應(yīng)曲線以及長(zhǎng)短周期運(yùn)動(dòng)之間的能量轉(zhuǎn)移現(xiàn)象等非線性特性.此外,由于粒桿模型將小行星視為由無(wú)質(zhì)量桿連接的質(zhì)量粒子組成,本文的研究成果可以拓展現(xiàn)有小行星附近周期軌道族.

1 受攝不規(guī)則小行星系統(tǒng)建模

本文采用粒桿模型模擬小行星的不規(guī)則引力場(chǎng),即根據(jù)小行星形狀特點(diǎn)及質(zhì)量分布特征,將小行星系統(tǒng)等效看作由3 個(gè)質(zhì)量粒子(m1,m2和m3) 組成,并假定m3m1=m2.3 個(gè)質(zhì)量粒子由兩個(gè)恒定距離的無(wú)質(zhì)量桿連接,在空間上呈等腰三角形構(gòu)型.小行星系統(tǒng)的質(zhì)心圍繞攝動(dòng)天體M4做圓形軌道運(yùn)動(dòng),且M4到小行星的距離遠(yuǎn)大于小行星內(nèi)部每個(gè)質(zhì)量粒子之間的距離.假定所有天體都在同一個(gè)平面內(nèi)運(yùn)動(dòng),且探測(cè)器質(zhì)量忽略不計(jì).

以小行星系統(tǒng)的質(zhì)心為原點(diǎn)O,建立小行星本體坐標(biāo)系(O-xyz),受攝粒桿模型示意圖如圖1 所示.Oz軸沿著小行星的自轉(zhuǎn)角速度? 方向,Ox軸平行于從粒子m1指向粒子m2的連線,Oy軸由右手坐標(biāo)系確定.將系統(tǒng)無(wú)量綱化,即設(shè)單元長(zhǎng)度為粒子m1與m2之間的距離(表示為L(zhǎng)),粒子m3與m1的距離在Oy軸上的投影為h,單位質(zhì)量為M=m1+m2+m3.引入μ=m1/M作為系統(tǒng)的質(zhì)量參數(shù),則各粒子m1,m2,m3和M4的質(zhì)量可以表示為

圖1 受攝粒桿模型示意圖Fig.1 Perturbed particle-linkage model

設(shè)定長(zhǎng)度比σ=h/L為距離參數(shù).當(dāng)粒子m3在小行星的質(zhì)心上方時(shí),σ 數(shù)值為正,否則為負(fù).由于小行星的質(zhì)心以恒定距離Lso繞M4做圓形軌道運(yùn)動(dòng),因此攝動(dòng)天體M4的位置可以寫為(Lsocos θ,Lsosin θ).其中,θ 是Ox軸與攝動(dòng)天體M4和小行星質(zhì)心的連線之間的角度(逆時(shí)針測(cè)量).探測(cè)器(x,y)在小行星本體坐標(biāo)系O-xy中的動(dòng)力學(xué)方程可以表示為

其中,k是無(wú)量綱參數(shù),代表小行星表面引力和離心力之比

各質(zhì)量粒子到探測(cè)器的距離ri可以展開為

定義攝動(dòng)天體的特征參數(shù)

其中,以無(wú)量綱參數(shù)k=1 為邊界,k>1 時(shí),引力大于離心力,小行星處于聚集狀態(tài);k< 1 時(shí),引力小于離心力,小行星處于離散狀態(tài).

觀察式(6)可知,攝動(dòng)天體M4對(duì)探測(cè)器的引力由兩部分組成:由于平均引力而產(chǎn)生的直接引力部分(以線性剛度項(xiàng)表示為kβx/2 和kβy/2),以及因攝動(dòng)天體位置的周期性變化激勵(lì)而產(chǎn)生的參數(shù)激勵(lì)部分.相應(yīng)地,類似著名的Mathieu 方程[30],式(6) 可被分為兩部分:第一部分,非擾動(dòng)系統(tǒng)即小行星系統(tǒng),表示為

其中,有效勢(shì)V表示為

以及第二部分,參數(shù)激勵(lì)部分,即式(6) 中包含具有非自治正弦函數(shù)和余弦函數(shù)的項(xiàng).值得注意的是,無(wú)擾動(dòng)系統(tǒng)(7)與經(jīng)典粒桿模型系統(tǒng)并不相同.無(wú)擾動(dòng)系統(tǒng)的有效勢(shì)函數(shù)包含攝動(dòng)天體的特征參數(shù)β,而β最終會(huì)影響系統(tǒng)的固有頻率.

在本研究中,攝動(dòng)小行星系統(tǒng)的平衡點(diǎn)坐標(biāo)(xo,yo)與小行星的質(zhì)量參數(shù)μ,距離參數(shù)σ,自旋參數(shù)k和攝動(dòng)參數(shù)β 均相關(guān).在給定參數(shù)[μ,σ,k,β]的情況下,在xy平面上生成Vx=?V/?x和Vy=?V/?y的零值曲線,從中可以得出交點(diǎn)的初始值.再通過Newton-Raphson 迭代法得到平衡點(diǎn)的精確位置[40].

當(dāng)β=0 時(shí),系統(tǒng)退化為經(jīng)典的粒桿模型.根據(jù)Yang 等[36]的研究,可得小行星433 Eros 和243 Ida 的參數(shù),按照上述迭代方法求解平衡點(diǎn),如圖2(a)和圖2(b)所示.小行星433 Eros 和243 Ida 都有4 個(gè)外部平衡點(diǎn).當(dāng)考慮攝動(dòng)天體的引力即β0 時(shí),小行星433 Eros 和243 Ida 的平衡點(diǎn)的位置將略有改變.

為了分析式(6)中的探測(cè)器動(dòng)力學(xué)行為,有必要將坐標(biāo)系的原點(diǎn)從原有的小行星本體坐標(biāo)系移動(dòng)到無(wú)擾動(dòng)系統(tǒng)的平衡點(diǎn)處.相應(yīng)的,單位長(zhǎng)度調(diào)整為平衡點(diǎn)與其最接近的質(zhì)量粒子(m1,m2或m3)之間的距離,表示為γo.這里,選擇任意一平衡點(diǎn)Lo作為原點(diǎn),在新的坐標(biāo)系中,ξ-和η-的方向與小行星本體坐標(biāo)系的方向平行,參見圖1.需要注意的是,選擇作為新坐標(biāo)系原點(diǎn)的平衡點(diǎn)應(yīng)該是穩(wěn)定的.對(duì)于以Lo?ξη(o=1,2,3,···)表示的以穩(wěn)定平衡點(diǎn)Lo為中心的坐標(biāo)系和小行星本體坐標(biāo)系O–xy之間的轉(zhuǎn)換關(guān)系可以表示為

圖2 平衡點(diǎn)隨參數(shù)變化的分布示意圖Fig.2 Distributions of equilibrium points with variations of parameters

根據(jù)式(9),將各粒子到探測(cè)器的距離ri(i=1,2,3)進(jìn)行Legendre 展開到四階,引入非線性項(xiàng),得到非線性系統(tǒng).進(jìn)而可將式(6)表示為

其中,有效勢(shì)函數(shù)V的第二階導(dǎo)數(shù)由下標(biāo)ξ 和η 表示,上下標(biāo)中的o 表示求導(dǎo)過程取在平衡點(diǎn)處.Ei,Fi,Di(i=1,2,3)定義為新坐標(biāo)系下平衡點(diǎn)與各質(zhì)量粒子的位置關(guān)系,其表達(dá)式列于附錄A 中.由于圖1 中系統(tǒng)模型關(guān)于y軸對(duì)稱,對(duì)于平衡點(diǎn)L2點(diǎn)和L4點(diǎn),有D1=D2,E1=?E2,E3=0,F1=F2.M1,M2,···,M7,N1,N2,···,N7是有效勢(shì)函數(shù)V經(jīng)Legendre 展開獲得的多項(xiàng)式系數(shù),其表達(dá)式見附錄A.ε 為小量符號(hào),用于標(biāo)記非線性項(xiàng)的量級(jí),在計(jì)算中取值為1.

2 參數(shù)激勵(lì)共振條件分析

根據(jù)Nayfeh 和Mook[34]的觀點(diǎn),參激共振發(fā)生有兩個(gè)重要的前提條件.首先,參激共振周期軌道設(shè)計(jì)問題應(yīng)在穩(wěn)定平衡點(diǎn)附近進(jìn)行分析.平衡點(diǎn)的穩(wěn)定性與其周圍軌道的穩(wěn)定性密切相關(guān),平衡點(diǎn)的穩(wěn)定性決定了其周圍不變流形的穩(wěn)定性.其次,參數(shù)激勵(lì)頻率與系統(tǒng)固有頻率應(yīng)存在相應(yīng)二倍關(guān)系.結(jié)合這兩個(gè)前提條件,可以找到發(fā)生參數(shù)共振的參數(shù)區(qū)域.下面依次分析兩條件.

條件一:平衡點(diǎn)穩(wěn)定性條件.對(duì)于ξ–η 分量,等式(10)中系統(tǒng)線性部分的特征方程為

其中

如果特征方程(12)具有純虛部或具有負(fù)實(shí)部的復(fù)數(shù)根,則平衡點(diǎn)將是穩(wěn)定的.因此,需滿足以下條件

條件二:參數(shù)激勵(lì)頻率與固有頻率關(guān)系條件.根據(jù)式(12),得到平衡點(diǎn)Lo附近的非擾動(dòng)系統(tǒng)的固有頻率為

此外,M4繞小行星系統(tǒng)質(zhì)心的圓周運(yùn)動(dòng)是二體問題,且滿足開普勒第三定律.在慣性坐標(biāo)系中,M4的圓周運(yùn)動(dòng)頻率?s可表示為

由于θ 是圖1 中x軸與從質(zhì)心到M4的線之間的夾角,因此,在平衡點(diǎn)坐標(biāo)系中,θ 相對(duì)于x軸旋轉(zhuǎn)的頻率為

這樣,可以將式(10)中的θ 改寫為θ=ω0t,即

其中,ω=2ω0.

通常,二維系統(tǒng)存在兩個(gè)固有頻率ω1和ω2(假設(shè)ω2> ω1),分別對(duì)應(yīng)于長(zhǎng)周期和短周期的軌道運(yùn)動(dòng).因此,在理論上可能存在第一類和第二類參激的兩種主共振類型,即滿足頻率關(guān)系

以及和型參激共振和差型參激共振的兩種組合共振,即滿足頻率關(guān)系

由于平衡點(diǎn)的位置、穩(wěn)定性以及固有頻率均與參數(shù)[μ,σ,k,β]有關(guān),而參數(shù)激勵(lì)頻率與k,β 參數(shù)相關(guān),因此,通過條件一和條件二,可獲得滿足參激共振條件的參數(shù)空間.

為了方便觀察攝動(dòng)參數(shù)β 對(duì)系統(tǒng)的影響,在仿真中固定參數(shù)σ 和k,研究系統(tǒng)參數(shù)μ與β 的取值范圍.圖3 為σ=0.8,k=0.9 時(shí),滿足條件一和條件二的μ與β 取值范圍.在圖3 中,灰色區(qū)域表示滿足平衡點(diǎn)穩(wěn)定條件的取值區(qū)域,紅線表示滿足第二類參激主共振條件的取值區(qū)域,藍(lán)線則表示滿足和型參激共振條件的取值區(qū)域.

圖3 滿足參激共振條件的參數(shù)取值區(qū)域σ=0.8,k=0.9Fig.3 Parametric resonance region when σ=0.8,k=0.9

圖4 不同系統(tǒng)參數(shù)下滿足參激共振條件的參數(shù)取值Fig.4 Parametric resonance regions with different k and σ

圖4 為更多組合取值變化條件下的分析結(jié)果,圖中紅線和藍(lán)線分別為滿足第二類參激主共振條件以及和型參激共振條件的共振線.為了觀察小行星自旋特性對(duì)系統(tǒng)共振區(qū)域的影響,固定參數(shù)σ,改變參數(shù)k,得到滿足不同參激共振情況的參數(shù)取值范圍,如圖4(a)所示.由圖4(a)可知,當(dāng)k=0.6 時(shí),系統(tǒng)只能發(fā)生第二類參激主共振,隨著參數(shù)k增大,系統(tǒng)同時(shí)出現(xiàn)第二類參激主共振以及和型參激共振.

同樣,為了觀察小行星形狀對(duì)系統(tǒng)共振區(qū)域的影響,固定參數(shù)k,改變參數(shù)σ,得到滿足不同參激共振情況的參數(shù)取值范圍,如圖4(b).在圖4(b)中,隨著參數(shù)σ 減小,滿足共振條件所對(duì)應(yīng)的系統(tǒng)μ值增大.因此,對(duì)于不規(guī)則形狀的小行星,形狀越接近細(xì)長(zhǎng)型,且處于聚集狀態(tài)時(shí),即形狀參數(shù)σ 越小,自旋參數(shù)k越大,系統(tǒng)更容易找到滿足參激共振的周期軌道族.需要說明的是,在數(shù)值模擬的過程中,沒有發(fā)現(xiàn)滿足第一類參激主共振和差型參激共振關(guān)系的參數(shù)取值區(qū)域.

結(jié)合圖3 和圖4 的分析可知,對(duì)于攝動(dòng)小行星系統(tǒng),普遍存在第二類參激主共振,在一定條件下,存在和型參激共振,但是普遍不存在第一類參激主共振和差型參激共振.

此外,本文分析了系統(tǒng)存在內(nèi)共振的可能性,即系統(tǒng)固有頻率之間存在關(guān)系

同樣,以σ=0.8 和k=0.9 時(shí)為例,結(jié)果如圖3 所示.圖中,綠線和黑線分別表示滿足1:2 和1:3 內(nèi)共振條件的參數(shù)范圍.

根據(jù)分析可知,系統(tǒng)存在同時(shí)發(fā)生參激共振以及內(nèi)共振的可能性.在此區(qū)域內(nèi),參數(shù)激勵(lì)共振使得系統(tǒng)在微小引力攝動(dòng)激勵(lì)條件下產(chǎn)生較大幅值的響應(yīng),即參激共振周期軌道,而內(nèi)共振則可能使得系統(tǒng)整體的能量在長(zhǎng)短周期的軌道運(yùn)動(dòng)之間轉(zhuǎn)換,由此產(chǎn)生更為復(fù)雜的非線性動(dòng)力學(xué)行為.

3 參數(shù)激勵(lì)共振軌道求解

基于上述分析可知,受攝小行星系統(tǒng)普遍存在內(nèi)共振以及第二類參激主共振,本節(jié)以第二類參激主共振和1:3 內(nèi)共振為例,對(duì)參數(shù)激勵(lì)共振軌道進(jìn)行求解.

為了確定非線性和參數(shù)激勵(lì)對(duì)系統(tǒng)振幅和相位的組合影響,采用多尺度方法,假定式(18)具有如下形式的解

其中,上述解均與快時(shí)間尺度T0=t以及慢時(shí)間尺度T2=ε2t有關(guān).約定算子

把式(22)代入式(18),比較ε 的同次冪,得到

不難發(fā)現(xiàn),式(24)為齊次方程.因此,式(24)的解可以表示為

其中,A1和A2分別為長(zhǎng)短周期軌道振幅,cc表示其前面各項(xiàng)的共軛

將式(27)代入式(24),可以獲得式(25)中ξ2和η2的解.考慮參激頻率和固有頻率間的關(guān)系以及1:3 的內(nèi)共振情況,引入調(diào)諧參數(shù)τ 和κ,將頻率之間的關(guān)系表示為

由此可得

為了確定式(26)的可解條件,尋求ξ3和η3的特解,以此來消除長(zhǎng)期項(xiàng),可設(shè)

然后,將ξ1,η1和ξ2,η2代入式(26)右端,式(31)代入到式(26)左端,并使兩側(cè)exp(iω1T0)和exp(iω2T0)的系數(shù)相等,得到

且W11,W22,S11,S22的表達(dá)式列于附錄B 中,上標(biāo)“—”表示相應(yīng)參數(shù)的共軛.

由此,確定式(27) 的可解性條件的問題被簡(jiǎn)化為確定線性代數(shù)方程(32)的系數(shù)行列式為零的條件,考慮式(28)中的關(guān)系,可得

把式(33)和式(34)代入式(35),得到

式中,G20與Gr1,Gr3(r=1,2)都是關(guān)于系統(tǒng)非線性系數(shù)M1,M2,...,M7,N1,N2,...,N7的復(fù)數(shù),分離其實(shí)虛部為R20與Rr1,Rr2,Rr3,Ir1,Ir2,Ir3,(r=1,2)表達(dá)式列于附錄B 中,上標(biāo)“—”表示相應(yīng)參數(shù)的共軛.

將式(36) 中的A1和A2表示為實(shí)數(shù)ar,βr的形式

結(jié)合歐拉公式,分離實(shí)部和虛部,并引入新的相位角

把關(guān)系式(39)代入式(36),得到一次近似解式(27)中振幅和相位應(yīng)滿足的自治微分方程

其中,方程(40)和(41)明確了系統(tǒng)非線性的影響.非線性項(xiàng)不僅可以直接影響振幅,還可以通過改變相位間接影響振幅.由此,可將系統(tǒng)的一次近似解表示為形式

并根據(jù)由式(28)、式(30) 和式(39) 所給出的相位關(guān)系

得到參激共振運(yùn)動(dòng)的一次近似解的通解形式

其中,ω 是參數(shù)激勵(lì)的頻率.需要注意的是,與式(27)不同,式(44)并不是系統(tǒng)的線性解,其包含瞬態(tài)和穩(wěn)態(tài)兩種運(yùn)動(dòng).振幅ar和相位φr(r=1,2)的取值必須嚴(yán)格遵循方程(40) 和(41),而式(27) 中系統(tǒng)線性解的振幅可以選擇為任意值.

此外,由式(44) 可知,長(zhǎng)周期軌道運(yùn)動(dòng)頻率為ω/6,而短周期軌道運(yùn)動(dòng)頻率為ω/2.由此可見,由于系統(tǒng)非線性的影響,使得短周期軌道的運(yùn)動(dòng)頻率與參數(shù)激勵(lì)頻率之間呈現(xiàn)精確的二倍關(guān)系,與此同時(shí),長(zhǎng)短周期運(yùn)動(dòng)之間呈現(xiàn)精確的1:3 內(nèi)共振.雖然共振響應(yīng)軌道同時(shí)包含長(zhǎng)短周期分量,但是由于頻率之間的諧振關(guān)系依然呈現(xiàn)周期性.

4 參數(shù)激勵(lì)共振軌道特性分析

4.1 系統(tǒng)穩(wěn)態(tài)解及其穩(wěn)定性判斷

由于通解式(44) 中包含有周期性的穩(wěn)態(tài)運(yùn)動(dòng),顯然,考慮到非線性項(xiàng)的影響后,根據(jù)式(40) 和式(41) 可知,系統(tǒng)可能存在定常振幅==0 和定常相位==0,即在攝動(dòng)天體的影響下,參數(shù)激勵(lì)共振響應(yīng)的振幅不再是無(wú)邊界,而是存在實(shí)現(xiàn)穩(wěn)態(tài)周期軌道運(yùn)動(dòng)的可能性.該軌道可為小行星探測(cè)任務(wù)的軌道設(shè)計(jì)提供參考.

根據(jù)式(40) 和式(41),穩(wěn)態(tài)運(yùn)動(dòng)的定常解振幅ar0和相位φr0應(yīng)滿足的代數(shù)方程

其中,a10,a20,φ10,φ20表示a1,a2,φ1,φ2對(duì)應(yīng)的穩(wěn)態(tài)解,由于模型的對(duì)稱性,I20,I21,I22,I23取值為0.

按照a10是否為零,上述方程組的解可分為平凡解和非平凡解兩種情況.當(dāng)a10取平凡解時(shí),系統(tǒng)退化為只與ω2相關(guān)的短周期運(yùn)動(dòng).本文考慮一般情況,即a10取非平凡解情況,a100,a200.由式(45)和式(47)可知

即系統(tǒng)存在的4 種相位組合情形.

給定系統(tǒng)參數(shù),采用Newton-Raphson 法針對(duì)不同相位組合情況可求得相應(yīng)振幅穩(wěn)態(tài)解a10,a20,從而得到振幅穩(wěn)態(tài)解與調(diào)諧參數(shù)τ 之間的關(guān)系圖,稱之為幅頻響應(yīng)曲線.由式(29)可知,調(diào)諧參數(shù)τ 可以用于表征參數(shù)激勵(lì)頻率的大小,由此,通過幅頻響應(yīng)曲線可以分析得知不同參數(shù)激勵(lì)頻率所對(duì)應(yīng)的系統(tǒng)振幅穩(wěn)態(tài)解.

此外,為了確保參激共振軌道的穩(wěn)定性,需要對(duì)其長(zhǎng)短周期運(yùn)動(dòng)進(jìn)行分析.根據(jù)式(45) 至式(48) 可知,長(zhǎng)周期運(yùn)動(dòng)穩(wěn)態(tài)解的穩(wěn)定性由式(45)和式(46)的奇點(diǎn)性質(zhì)判定,短周期運(yùn)動(dòng)穩(wěn)態(tài)解的穩(wěn)定性由式(47)和式(48)的奇點(diǎn)性質(zhì)判定.考慮到式(49),參激共振軌道的穩(wěn)定性可以通過式(46)和式(48)的奇點(diǎn)性質(zhì)來判斷.對(duì)式(46)和式(48)的奇點(diǎn)性質(zhì)進(jìn)行分析,得到用于判斷長(zhǎng)短周期振幅a10,a20穩(wěn)定性的Jacobi 矩陣為

矩陣A的特征值可以表示為

其中,p是A的跡,q是A的行列式

當(dāng)p2?4q> 0 時(shí),λ1,2為互異實(shí)特征值.若q< 0,則λ1λ2< 0,奇點(diǎn)為鞍點(diǎn).若q> 0,p< 0,則λ1λ2>0,奇點(diǎn)為漸近穩(wěn)定結(jié)點(diǎn).當(dāng)p2?4q< 0 時(shí),λ1,2為共軛復(fù)特征值.若p< 0,奇點(diǎn)為漸近穩(wěn)定焦點(diǎn).若p> 0,奇點(diǎn)為不穩(wěn)定焦點(diǎn).綜上所述,得到奇點(diǎn)漸近穩(wěn)定的判別條件,即同時(shí)滿足參激共振軌道長(zhǎng)短周期振幅穩(wěn)定性的條件如下

下面通過數(shù)值模擬,對(duì)系統(tǒng)穩(wěn)態(tài)解的性質(zhì)進(jìn)一步分析,μ=0.001 13,σ=0.8,k=0.9,β=0.014 為例,選取穩(wěn)定的平衡點(diǎn)L4,對(duì)4 種相位組合情況分別討論,繪制幅頻響應(yīng)曲線,如圖5 和圖6 所示.

圖5 幅頻響應(yīng)曲線(情況1、情況3)Fig.5 Frequency-amplitude response curves(Case 1,Case 3)

圖6 幅頻響應(yīng)曲線(情況2、情況4)Fig.6 Frequency-amplitude response curves(Case 2,Case 4)

圖6 幅頻響應(yīng)曲線(情況2、情況4)(續(xù))Fig.6 Frequency-amplitude response curves(Case 2,Case 4)(continued)

圖5 為情況1(φ10=0,φ20=0)和情況3(φ10=π,φ20=0)的幅頻響應(yīng)曲線.在圖5 中,兩種情況的曲線重合,黑線為長(zhǎng)周期振幅穩(wěn)態(tài)解a10,藍(lán)線為短周期振幅a20.不難發(fā)現(xiàn),穩(wěn)態(tài)解a10,a20均呈現(xiàn)隨調(diào)諧參數(shù)τ 的增大而單調(diào)遞增的趨勢(shì),即一個(gè)參數(shù)激勵(lì)頻率只對(duì)應(yīng)于一組穩(wěn)態(tài)解.當(dāng)調(diào)諧參數(shù)τ=0.000 7 時(shí),系統(tǒng)存在一組穩(wěn)態(tài)解a10=0.014 2,a20=0.000 3,如圖5中紅點(diǎn)所示.由式(52)和式(54)可得,a10=0.014 2,a20=0.000 3 時(shí),p=?0.010 6 < 0,q=?0.004 8 < 0,λ1=0.004 3,λ2=?0.074 7,該組穩(wěn)態(tài)解是不穩(wěn)定的.通過式(54) 驗(yàn)證穩(wěn)態(tài)解的穩(wěn)定性,發(fā)現(xiàn)在幅頻響應(yīng)曲線中,任何一組穩(wěn)態(tài)解所對(duì)應(yīng)的q值均小于0,穩(wěn)態(tài)解具有不穩(wěn)定性.因此,在設(shè)計(jì)參激共振周期軌道時(shí),情況1 和情況3 將不予以考慮.

圖6 為相位組合關(guān)系為情況2(φ10=0,φ20=π)和情況4(φ10=π,φ20=π)的幅頻響應(yīng)曲線.當(dāng)調(diào)諧參數(shù)τ=0.000 7,系統(tǒng)存在一組穩(wěn)態(tài)解a10=0.013 9,a20=0.000 4,如圖6 中紅點(diǎn)所示.基于式(52) 和式(54)可得,p=?0.114 1<0,q=0.003 1>0,p2?4q>0,λ1=?0.056 7,λ2=?0.070 0,該組穩(wěn)態(tài)解是穩(wěn)定的.為了進(jìn)一步描述調(diào)諧參數(shù)τ 對(duì)穩(wěn)態(tài)解的穩(wěn)定性的影響,將幅頻響應(yīng)曲線分為3 個(gè)區(qū)域,標(biāo)記為I,II,III,對(duì)長(zhǎng)短周期振幅變化軌跡分別進(jìn)行標(biāo)記.

當(dāng)τ 從0 增大到0.02 時(shí),長(zhǎng)周期振幅a10經(jīng)過A1至A4點(diǎn),軌跡方向由黑色箭頭表示;短周期振幅a20經(jīng)過B1至B4點(diǎn),軌跡方向由藍(lán)色箭頭表示,如圖6(a)所示.在A1、B1點(diǎn)(τ=0)到A2,B2點(diǎn)(τ=0.010 4)的區(qū)域I 內(nèi),系統(tǒng)只存在一組非平凡解,且通過式(46)和式(48)可分別求出.基于式(54)可得,p<0,q>0,即區(qū)域I 中的各組穩(wěn)態(tài)解均是穩(wěn)定的,如圖6 區(qū)域I中的實(shí)線所示.在A2,B2點(diǎn)到A3,B3點(diǎn)(τ=0.012 3)的區(qū)域II 中,長(zhǎng)周期振幅a10和短周期振幅a20分別存在兩組非平凡解.一組非平凡解所對(duì)應(yīng)的q值小于0,為不穩(wěn)定穩(wěn)態(tài)解,如圖6 區(qū)域II 中的虛線所示.一組非平凡解的p< 0,q> 0,則該組穩(wěn)態(tài)解是穩(wěn)定的,如圖6 區(qū)域II 中的實(shí)線所示.但奇點(diǎn)種類有所不同,當(dāng)τ ∈(0.010 4,0.011 8)時(shí),p2?4q> 0,λ1,λ2為負(fù)實(shí)數(shù),呈現(xiàn)“穩(wěn)定結(jié)點(diǎn)”特性,當(dāng)τ ∈(0.011 8,0.012 3)時(shí),p2?4q< 0,λ1,λ2為實(shí)部為負(fù)的復(fù)數(shù),呈現(xiàn)“穩(wěn)定焦點(diǎn)”特性.特別注意的是,在A3,B3點(diǎn),兩個(gè)非平凡解不再以原軌跡變化,出現(xiàn)“跳躍現(xiàn)象”,a10從A3點(diǎn)跳躍到A4點(diǎn),a20從B3點(diǎn)跳躍到B4點(diǎn).在區(qū)域III 內(nèi),即τ ∈(0.012 3,+∞)時(shí),p< 0,q< 0,呈現(xiàn)“鞍點(diǎn)”特征,該組穩(wěn)態(tài)解是不穩(wěn)定的,如圖6 區(qū)域III 中的虛線所示.

當(dāng)τ 從0.02 減小到0 時(shí),長(zhǎng)周期振幅a10經(jīng)過C1至C4點(diǎn),軌跡方向由黑色箭頭表示; 短周期振幅a20經(jīng)過D1至D4點(diǎn),軌跡方向由藍(lán)色箭頭表示,如圖6(b)所示.在τ 經(jīng)過C2,D2點(diǎn)到達(dá)C3,D3點(diǎn)之前,系統(tǒng)存在兩組非平凡解,之后可以觀察到跳躍到C4,D4點(diǎn)的現(xiàn)象.與上述τ 增大過程類似,從C4,D4點(diǎn)之后,唯一可實(shí)現(xiàn)的穩(wěn)定的非平凡解由式(46) 和式(48) 給出.除此之外,系統(tǒng)不存在穩(wěn)定非平凡解.因此,調(diào)諧參數(shù)的變化方向不同,使得系統(tǒng)穩(wěn)定非平凡解所對(duì)應(yīng)的調(diào)諧參數(shù)區(qū)域也有所不同.

結(jié)合上述長(zhǎng)短周期振幅變化的特點(diǎn),發(fā)現(xiàn)了一個(gè)有趣的現(xiàn)象.在穩(wěn)定區(qū)域I 和區(qū)域II 中,隨著調(diào)諧參數(shù)τ 的增大,長(zhǎng)周期振幅a10先增大后減小,而短周期振幅a20不斷增大,直至超過長(zhǎng)周期振幅a10.根據(jù)式(29) 得知,參數(shù)激勵(lì)ω 與系統(tǒng)固有頻率ω2之間存在共振關(guān)系,通過參數(shù)激勵(lì),能量不斷輸入系統(tǒng),激發(fā)系統(tǒng)短周期運(yùn)動(dòng)(振幅a20相關(guān)運(yùn)動(dòng)).在區(qū)域I的初始階段,內(nèi)共振調(diào)諧參數(shù)κ < 0.005,使得在區(qū)域I 內(nèi)共振的影響顯著,作用于短周期運(yùn)動(dòng)的能量通過內(nèi)共振不斷被轉(zhuǎn)移到長(zhǎng)周期運(yùn)動(dòng),進(jìn)而出現(xiàn)長(zhǎng)周期振幅a10為主的穩(wěn)態(tài)運(yùn)動(dòng).隨著調(diào)諧參數(shù)τ 不斷增大,內(nèi)共振調(diào)諧參數(shù)κ 也隨之增大,內(nèi)共振能量轉(zhuǎn)移效應(yīng)減弱,從而在區(qū)域II 內(nèi)短周期運(yùn)動(dòng)的振幅a20不斷增大直至超過長(zhǎng)周期振幅a10.

為了進(jìn)一步描述參數(shù)激勵(lì)能量(引力攝動(dòng)β)對(duì)系統(tǒng)的影響程度,以μ=0.001 13,σ=0.8,k=0.9為例,固定調(diào)諧參數(shù)τ=0.000 7,繪制參數(shù)激勵(lì)振幅1.5kβ 和參激響應(yīng)振幅ar0的力頻關(guān)系曲線,如圖7所示.

圖7 力頻響應(yīng)曲線Fig.7 Amplitude of the excitation–amplitude of response curves

在情況1 和情況3 中,β ∈[0,0.016 0],穩(wěn)態(tài)解隨參激振幅的增大而減小.p< 0,q< 0,穩(wěn)態(tài)解呈現(xiàn)“鞍點(diǎn)”特性,該組穩(wěn)態(tài)解是不穩(wěn)定的,如圖7(a)中虛線所示.力頻關(guān)系分析結(jié)果和幅頻響應(yīng)曲線結(jié)果一致,因此,在設(shè)計(jì)參激共振周期軌道時(shí),情況1 和情況3 將不予以考慮.

在情況2 和情況4 中,類似于幅頻響應(yīng)曲線,根據(jù)穩(wěn)態(tài)解的個(gè)數(shù)和穩(wěn)定性,可以將力頻響應(yīng)曲線分為3 個(gè)區(qū)域,如圖7(b)中虛線所示.其中,實(shí)線表示穩(wěn)定的穩(wěn)態(tài)解,虛線表示不穩(wěn)定的穩(wěn)態(tài)解.區(qū)域邊界(紅線)為β=0.003 8 (1.5kβ=0.005 1)和β=0.006 9(1.5kβ=0.009 3).當(dāng)β ∈[0,0.003 8) 時(shí),長(zhǎng)周期振幅a10隨參激振幅的增大而減小,短周期振幅a20則呈現(xiàn)不斷增大,p< 0,q> 0,該穩(wěn)態(tài)解是穩(wěn)定的.當(dāng)β ∈[0.003 8,0.006 9] 時(shí),出現(xiàn)兩組穩(wěn)態(tài)解,一組穩(wěn)態(tài)解p<0,q<0,是不穩(wěn)定的;一組穩(wěn)態(tài)解p<0,q>0,是穩(wěn)定的.當(dāng)β ∈(0.006 9,0.016 0]時(shí),短周期振幅a20隨參激振幅的增大而減小,長(zhǎng)周期振幅a10則呈現(xiàn)先增大后減小的特點(diǎn),p< 0,q> 0,對(duì)應(yīng)的穩(wěn)態(tài)解是不穩(wěn)定的.力頻關(guān)系分析結(jié)果和幅頻響應(yīng)曲線結(jié)論一致.

4.2 參數(shù)激勵(lì)共振響應(yīng)軌道

通過上述分析判斷,根據(jù)參激共振響應(yīng)(44) 以及所求穩(wěn)定的穩(wěn)態(tài)解,最終得到了滿足第二主參數(shù)共振和1:3 內(nèi)共振條件下穩(wěn)定的參激周期軌道.由于在穩(wěn)定區(qū)域I 和區(qū)域II 中,探測(cè)器的動(dòng)力學(xué)行為有所差別.因此,將用兩組數(shù)據(jù)分別對(duì)參激共振軌道進(jìn)行模擬求解.

以參數(shù)μ=0.001 13,σ=0.8,k=0.9,β=0.014為例,得到調(diào)諧參數(shù)τ=0.000 7,κ=0.004 1,參激周期軌道長(zhǎng)短周期的振幅a10=0.013 9,a20=0.000 4.

圖8(a)中紅線和藍(lán)線分別表示長(zhǎng)短周期運(yùn)動(dòng)振幅隨時(shí)間的變化情況,黑線表示長(zhǎng)短周期疊加運(yùn)動(dòng)振幅隨時(shí)間變化曲線.不難發(fā)現(xiàn),由于非線性效應(yīng),長(zhǎng)短周期軌道頻率之間呈現(xiàn)精確的1:3 內(nèi)共振比例關(guān)系.

圖8(b)中,藍(lán)線表示攝動(dòng)天體隨時(shí)間t在x方向的運(yùn)動(dòng)分量,黑線表示參數(shù)激勵(lì)響應(yīng)即探測(cè)器周期軌道隨時(shí)間t在x方向的運(yùn)動(dòng)分量.由于在該組參數(shù)條件下,內(nèi)共振引起的能量轉(zhuǎn)移效應(yīng)顯著,長(zhǎng)周期軌道振幅遠(yuǎn)遠(yuǎn)大于短周期軌道振幅.雖然共振響應(yīng)中同時(shí)包含有長(zhǎng)短周期的分量,但是參激周期軌道幾乎完全呈現(xiàn)為長(zhǎng)周期運(yùn)動(dòng),如圖8(c)所示.此外,由于非線性因素影響,攝動(dòng)天體頻率之間精確保持為1:6的關(guān)系.

以參數(shù)μ=0.001,σ=0.8,k=0.9,β=0.012 為例,得到調(diào)諧參數(shù)τ=0.000 6,κ=0.076 7,以及兩組穩(wěn)態(tài)解,(1)a10=0.023 5,a20=0.036 9,p< 0,q> 0,p2?4q< 0,穩(wěn)態(tài)解呈穩(wěn)定特性; (2)a10=0.013 2,a20=0.049 9,p< 0,q< 0,p2?4q> 0,穩(wěn)態(tài)解呈現(xiàn)不穩(wěn)定特性.取穩(wěn)定穩(wěn)態(tài)解a10=0.023 5,a20=0.036 9為參激周期軌道長(zhǎng)短周期的振幅,繪制共振軌道如圖9 所示.圖9(a)中紅線和藍(lán)線分別表示長(zhǎng)短周期運(yùn)動(dòng)振幅隨時(shí)間的變化情況,黑線表示長(zhǎng)短周期疊加運(yùn)動(dòng)振幅隨時(shí)間變化曲線.由圖可知,該組參數(shù)條件下長(zhǎng)短周期軌道振幅量級(jí)相當(dāng),且由于非線性效應(yīng),長(zhǎng)短周期軌道頻率之間呈現(xiàn)精確的1:3 內(nèi)共振比例關(guān)系.圖9(b) 中藍(lán)線和黑線分別表示攝動(dòng)天體和探測(cè)器周期軌道隨時(shí)間t在x方向的運(yùn)動(dòng)分量.顯然共振響應(yīng)呈現(xiàn)為長(zhǎng)短周期疊加運(yùn)動(dòng),且由于1:3 內(nèi)共振的影響,使得探測(cè)器周期軌道振幅呈現(xiàn)亞諧共振狀態(tài).此外由于長(zhǎng)短周期軌道頻率之間的諧振關(guān)系,最終呈現(xiàn)的響應(yīng)依然呈現(xiàn)周期性,如圖9(c)所示.

圖8 當(dāng)σ=0.8,k=0.9,μ=0.001 13,β=0.014 時(shí)參激共振軌道示意圖Fig.8 Stable parametric resonance orbit when σ=0.8,k=0.9,μ=0.001 13,β=0.014

圖9 當(dāng)σ=0.8,k=0.9,μ=0.001,β=0.012 時(shí)參激共振軌道示意圖Fig.9 Stable parametric resonance orbit when σ=0.8,k=0.9,μ=0.001,β=0.012

通過上述模擬分析,不難發(fā)現(xiàn),參激周期軌道與系統(tǒng)參數(shù)選取密切相關(guān).最終參激周期軌道的形狀以及周期關(guān)系將嚴(yán)格遵循式(44)的表達(dá)形式.

值得一提的是:μ,σ,k是與小行星自身特征相關(guān)的參數(shù),β 與攝動(dòng)天體的引力擾動(dòng)有關(guān).當(dāng)參數(shù)滿足參激共振條件一和條件二時(shí),系統(tǒng)可得到一個(gè)第二類參激主共振軌道或一個(gè)和型參激共振軌道.此外,參激共振的條件頗為苛刻,并非所有受攝小行星系統(tǒng)都能滿足該條件.本文所提出的方法為軌道設(shè)計(jì)中處理引力攝動(dòng)提供了新的思路.

5 結(jié)論

本文將太陽(yáng)引力視為受攝小行星系統(tǒng)的組成部分,采用非線性振動(dòng)理論中的參激共振概念,創(chuàng)新性地設(shè)計(jì)了小行星平衡點(diǎn)附近穩(wěn)定的參數(shù)共振軌道.

(1)通過引入質(zhì)量參數(shù)μ,形狀參數(shù)σ,自旋參數(shù)k以及攝動(dòng)參數(shù)β,對(duì)受攝小行星系統(tǒng)進(jìn)行建模.研究發(fā)現(xiàn),系統(tǒng)平衡點(diǎn)位置以及穩(wěn)定性、平衡點(diǎn)附近運(yùn)動(dòng)的固有頻率均與4 個(gè)參數(shù)相關(guān).

(2)通過固定參數(shù)σ 和k,定性分析了攝動(dòng)參數(shù)β對(duì)系統(tǒng)的影響.對(duì)于攝動(dòng)小行星系統(tǒng)而言,普遍存在第二類參激主共振,在一定條件下,存在和型參激共振,但是普遍不存在第一類參激主共振和差型參激共振,且系統(tǒng)普遍存在1:2 和1:3 內(nèi)共振.

(3)以第二類參激主共振和1:3 內(nèi)共振為例,采用多尺度方法得到參激共振響應(yīng)的一次近似解,求解了振幅和相位的穩(wěn)態(tài)解并分析其穩(wěn)定性,得到由參激共振引起的穩(wěn)定周期軌道.通過受攝小行星系統(tǒng)的幅頻響應(yīng)曲線和力頻響應(yīng)曲線分析了系統(tǒng)非線性特性.通過分析可知,由于系統(tǒng)非線性的影響,短周期軌道的運(yùn)動(dòng)頻率與參數(shù)激勵(lì)頻率之間呈現(xiàn)精確的二倍關(guān)系,與此同時(shí),長(zhǎng)短周期運(yùn)動(dòng)之間呈現(xiàn)精確的1:3 內(nèi)共振.雖然共振響應(yīng)軌道同時(shí)包含長(zhǎng)短周期分量,但是由于頻率之間的諧振關(guān)系,使得軌道依然呈現(xiàn)周期性.

此外,由于粒桿模型將小行星視為由無(wú)質(zhì)量桿連接的質(zhì)量粒子組成,因此本文提出的參激共振軌道設(shè)計(jì)方法同樣適用于雙小行星系統(tǒng)或三小行星系統(tǒng)中參數(shù)共振軌道動(dòng)力學(xué)問題.本文的研究成果可以拓展現(xiàn)有小行星附近周期軌道族,為軌道設(shè)計(jì)中處理引力攝動(dòng)提供了新的思路.

附錄A

附錄B