含輻射和扁率的圓型限制性三體問(wèn)題的軌道穩(wěn)定性研究?

石紹伍 馬大柱

(湖北民族大學(xué)信息工程學(xué)院恩施445000)

1 引言

圓型限制性三體問(wèn)題描述的是在三體中把第3體的質(zhì)量看成無(wú)限小,不考慮第3體對(duì)主星和伴星的引力作用,伴星在主星的作用下近似作圓周運(yùn)動(dòng).該系統(tǒng)不可積,但存在一個(gè)雅克比積分常數(shù)和5個(gè)平動(dòng)點(diǎn)解,具有廣泛的天文應(yīng)用背景.如在研究月球運(yùn)動(dòng)時(shí),不考慮太陽(yáng)軌道偏心率、太陽(yáng)視差和月球軌道傾角,這就是一個(gè)平面圓型限制性三體問(wèn)題.然而,大多數(shù)行星由于其赤道半徑和極半徑并不相等,應(yīng)作為扁球處理.另一方面,太陽(yáng)作為輻射源產(chǎn)生的輻射作用不能忽略.因此,在經(jīng)典模型基礎(chǔ)上加入主星的輻射作用以及伴星的扁狀攝動(dòng)引力勢(shì)對(duì)第3體的影響更接近實(shí)際太陽(yáng)系模型.由于在旋轉(zhuǎn)坐標(biāo)系下,第3體除受到引力作用之外,還受到旋轉(zhuǎn)坐標(biāo)系產(chǎn)生的非慣性力的影響[1].因此,含輻射和扁球的圓型限制性三體問(wèn)題的哈密頓量實(shí)際上可分解為動(dòng)能對(duì)應(yīng)的Euler流、由坐標(biāo)旋轉(zhuǎn)導(dǎo)致的正交旋轉(zhuǎn)和勢(shì)能對(duì)應(yīng)的Euler流.為減少數(shù)值誤差對(duì)長(zhǎng)期積分的影響,數(shù)值積分工具宜采用保辛結(jié)構(gòu)算法和保流形結(jié)構(gòu)算法.

辛算法可用于研究哈密頓系統(tǒng)長(zhǎng)期演化問(wèn)題.辛算法有顯式和隱式之分,顯式方法常用于處理可分離式哈密頓系統(tǒng),但在右函數(shù)形式復(fù)雜時(shí)其計(jì)算效率較低.隱式方法常需與梯形法、龍格庫(kù)塔法、低階顯辛法結(jié)合使用,計(jì)算過(guò)程較為繁瑣,一般不宜采用.辛算法還可以分為力梯度辛算法和非力梯度辛算法.Ruth[2]、Wisdom等[3?4]、Forest和Ruth[5]以及Yoshida[6]等作者構(gòu)造了多種非力梯度辛算法.力梯度辛算法主要的代表作者是Chin[7].國(guó)內(nèi)研究人員也設(shè)計(jì)了一種含3階導(dǎo)數(shù)項(xiàng)的4階力梯度辛算法.該文指出力梯度辛算法在計(jì)算限制性三體問(wèn)題時(shí)比非力梯度算法效果要好,且4階力梯度辛算法明顯優(yōu)于非力梯度辛算法和其他常用的力梯度辛算法[8].顯然,上述結(jié)論可以推廣到含輻射以及扁球的圓型限制性三體問(wèn)題.

基于最小二乘法原理的流形改正方法利用標(biāo)度因子將每一步數(shù)值積分結(jié)果強(qiáng)迫拉回到初始積分超曲面上,有效減少了數(shù)值積分時(shí)局部截?cái)嗾`差的積累和人工耗散.此類方法無(wú)需對(duì)哈密頓量進(jìn)行分解,操作簡(jiǎn)單效率高,已在相對(duì)論天體力學(xué)、太陽(yáng)系動(dòng)力學(xué)、量子力學(xué)和宇宙學(xué)中[9]得到廣泛應(yīng)用.近年來(lái),研究人員相繼發(fā)展了多種流形改正方法[10–18],如單因子法[10]、雙因子法[11]、速度因子方法[12]等等.然而將大多數(shù)流形改正方法直接應(yīng)用到非開(kāi)普勒問(wèn)題卻存在一些困難,原因在于大多數(shù)方法均要求嚴(yán)格或近似地滿足開(kāi)普勒能量、拉普拉斯積分、角動(dòng)量關(guān)系.但是,只要某系統(tǒng)存在任意一個(gè)守恒量,仍然能夠使用速度因子方法.文獻(xiàn)[19–27]表明速度因子方法在攝動(dòng)開(kāi)普勒問(wèn)題、H′enon-Heiles系統(tǒng)、致密雙星問(wèn)題、標(biāo)量場(chǎng)模型等問(wèn)題中均可以改進(jìn)低階算法在精度和穩(wěn)定性方面的不足,同時(shí)在CPU占用率方面也遠(yuǎn)遠(yuǎn)低于高階算法.簡(jiǎn)言之,速度因子方法具有操作簡(jiǎn)單、耗時(shí)少、應(yīng)用范圍廣的特點(diǎn).對(duì)于含輻射和扁球的圓型限制性三體問(wèn)題,借助積分不變關(guān)系[28],也可以改正每一體軌道根數(shù)的精度.

辛算法和流形改正方法均可以長(zhǎng)期保持系統(tǒng)的能量積分,目前無(wú)法從理論上區(qū)分(如代數(shù)精度或泰勒展開(kāi)式)哪種方式占優(yōu).若在相同代數(shù)精度下,基于4階龍格庫(kù)塔法的速度因子算法在數(shù)值精度、穩(wěn)定性和計(jì)算效率方面明顯優(yōu)于4階力梯度辛算法.因此,本文利用速度因子方法討論主星為輻射源,伴星為扁球的平面圓型限制性三體問(wèn)題的動(dòng)力學(xué)問(wèn)題.

2 物理模型

與經(jīng)典平面圓型限制性三體問(wèn)題一樣,假定第3體的質(zhì)量無(wú)限小,第3體對(duì)兩主天體的運(yùn)動(dòng)沒(méi)有任何影響.約定兩主天體的總質(zhì)量為1,伴星的質(zhì)量為m2=μ,則主星的質(zhì)量為m1=1?μ.旋轉(zhuǎn)坐標(biāo)系中坐標(biāo)原點(diǎn)選為兩主天體的質(zhì)心,x軸為兩主天體的連線.如引言所述,第3體受到兩主天體的引力作用、主星的輻射和伴星的扁狀攝動(dòng)影響,其運(yùn)動(dòng)方程[1]為:

此處q ∈[0,1]、A2∈[0,0.2]、n、?、r1和r2分別為主星的輻射因子、伴星的扁率、平運(yùn)動(dòng)角頻率、等效勢(shì)、第3體到主星的距離和第3體到伴星的距離.q= 1且A2= 0時(shí)即為經(jīng)典平面圓形限制型三體問(wèn)題.此時(shí)系統(tǒng)存在一個(gè)雅克比積分常數(shù):

為討論方便,本文僅限于2維情況(令z= 0).由于?(x,y)= ?(x,±y),勢(shì)能曲線關(guān)于x軸對(duì)稱.平動(dòng)點(diǎn)位于零速度曲線上(),即:

系統(tǒng)存在3個(gè)共線平動(dòng)點(diǎn)L1、L2和L3,兩個(gè)三角平動(dòng)點(diǎn)L4和L5.平動(dòng)點(diǎn)的穩(wěn)定性與經(jīng)典情況類似,此處不再討論.第3體的運(yùn)動(dòng)也可以用哈密頓函數(shù)表示:

其中px和py分別表示第3體在x和y方向動(dòng)量的分量.力函數(shù)的形式為:

哈密頓式可以分解為動(dòng)能T和勢(shì)能V兩部分:

因動(dòng)能中含有旋轉(zhuǎn)坐標(biāo)系為非慣性系而對(duì)第3體附加的影響ypx+xpy,導(dǎo)致動(dòng)能不是動(dòng)量的嚴(yán)格2次型.但若將位置算子變?yōu)槲恢门c動(dòng)量混合型算子,則可以使用力梯度辛算法求解[29].

3 數(shù)值方法

3.1 流形改正原理

對(duì)于含一個(gè)運(yùn)動(dòng)積分?(x)=C的n′維微分動(dòng)力系統(tǒng),其中x是坐標(biāo)和動(dòng)量的函數(shù)(x=x(p,q)),C、p和q分別代表常數(shù)、動(dòng)量和坐標(biāo).理論上應(yīng)滿足下述關(guān)系,

［本刊訊］2012年2月28日，上海市護(hù)理學(xué)會(huì)組織有關(guān)護(hù)理專家召開(kāi)了《上海市護(hù)理事業(yè)發(fā)展規(guī)劃(2011－2015年)》研討會(huì)。會(huì)上，來(lái)自上海各大醫(yī)療機(jī)構(gòu)的護(hù)理專家系統(tǒng)學(xué)習(xí)了《中國(guó)護(hù)理事業(yè)發(fā)展規(guī)劃綱要(2011－2015年)》，就上海市未來(lái)5年護(hù)理事業(yè)發(fā)展與實(shí)施方案進(jìn)行了熱烈討論，并向相關(guān)衛(wèi)生行政職能部門提出了積極有效的建議。

但數(shù)值計(jì)算時(shí)很難滿足?(x)= 0.按照某數(shù)值算法得到t時(shí)刻數(shù)值解為η=η(p′,q′),標(biāo)準(zhǔn)值和數(shù)值解的差記為??,則有

定義矩陣E,

根據(jù)最小二乘法原理得到改正矢量?η為

上標(biāo)T代表轉(zhuǎn)置,W是任意給定的權(quán)重矩陣.?η也可以寫成另外一種形式

將改正矢量?η作用在低階算法得到的數(shù)值解η后,可得新的改正解ξ為

將改正解ξ作為下一步積分的初值,如此循環(huán).理論上,改正后值ξ的精度為改正前η的平方倍,即:

該方法是Nacozy在1971年提出的[30].

3.2 速度因子方法

流形改正原理提供了一種處理數(shù)值解誤差的近似方法,但由于坐標(biāo)和速度的量綱不同,(15)式的物理意義并不明確.但如果只考慮對(duì)動(dòng)量(速度)的改正,就可以解決量綱問(wèn)題.并且Wu等[13]已嚴(yán)格證明只改正速度與同時(shí)改正速度和坐標(biāo)是等效的.速度因子方法[12]正是在此基礎(chǔ)之上提出來(lái)的,該方法定義一個(gè)拉格朗日乘子γ將數(shù)值積分時(shí)偏離的部分拉回到初始流形.γ的表達(dá)式為:

注意γ可取正根或負(fù)根,大樣本數(shù)值模擬表明負(fù)根容易出錯(cuò),導(dǎo)致不能長(zhǎng)期積分,故推薦采用正根.該方法將作為基本數(shù)值工具用于研究本文模型的動(dòng)力學(xué)問(wèn)題.

4 軌道穩(wěn)定性研究

速度因子方法在每一步積分后將偏離值以最近的距離拉回到積分初值附近,這樣可以嚴(yán)格保證積分常數(shù)始終在初值附近波動(dòng),從而保證數(shù)值積分精度.模型中的速度因子為:伴星的質(zhì)量為μ= 0.001 (參照木星與太陽(yáng)的質(zhì)量比).一般而言,主星的輻射和引力作用對(duì)行星或其他小天體的運(yùn)動(dòng)占主導(dǎo)[1],伴星的扁球攝動(dòng)影響較小.以太陽(yáng)-木星-小行星系統(tǒng)為例,可取輻射因子q=0.8,扁率A2=0.005.本文主要討論輻射因子和扁率的存在對(duì)軌道穩(wěn)定性的影響.

4.1 單條軌道穩(wěn)定性研究

Poincar′e截面只適用于研究不超過(guò)4維的保守系統(tǒng),可以直觀反映系統(tǒng)的動(dòng)力學(xué)性質(zhì).嚴(yán)格周期運(yùn)動(dòng)在截面上只出現(xiàn)一個(gè)點(diǎn); 擬周期運(yùn)動(dòng)會(huì)得到一些封閉曲線; 而混沌運(yùn)動(dòng)則出現(xiàn)一些呈隨機(jī)分布的點(diǎn).首先取一條軌道(軌道1):Cj= 3.06,x0= 0.29.若不考慮太陽(yáng)輻射和扁狀攝動(dòng),該軌道為擬周期軌道[29].計(jì)算結(jié)果見(jiàn)圖1 (a),截面上出現(xiàn)了由5條封閉曲線構(gòu)成的島鏈結(jié)構(gòu).如果只加入扁狀攝動(dòng),令q= 1,A2在0–0.2每間隔0.001變化,我們發(fā)現(xiàn)只要滿足條件0< A2<0.1,軌道總是處于混沌狀態(tài).圖1 (b)是其中選取的一個(gè)樣本,此時(shí)q= 1,A2= 0.009.圖中出現(xiàn)了許多雜亂無(wú)章的點(diǎn),原本有序的軌道變成混沌軌道.當(dāng)A20.1時(shí),扁狀攝動(dòng)太大使得軌道逃逸.如果只考慮輻射因子的影響,保持A2= 0不變,令q在區(qū)間[0,1]內(nèi)每間隔0.01緩慢變化.模擬發(fā)現(xiàn)當(dāng)q0.48時(shí),所有軌道均有序.圖1 (c)所示軌道為q= 0.9,該軌道是擬周期運(yùn)動(dòng)軌道.當(dāng)q <0.48時(shí),輻射太強(qiáng)使得軌道逃逸.同時(shí)具有輻射和扁狀攝動(dòng)時(shí),所有軌道運(yùn)動(dòng)均有序(要求q0.48,A2<0.1).圖1 (d)給出的q= 0.8,A2= 0.005即為這種情況.再考慮另外一條軌道(軌道2):Cj=3.06,x0=0.48.在無(wú)輻射因子和扁率的情況下該軌道為混沌軌道,截面上留下1000個(gè)雜亂無(wú)章的點(diǎn),如圖2 (a)所示.在無(wú)輻射(q= 1)但有扁率的情況下,當(dāng)0< A2<0.1時(shí),軌道混沌.圖2 (b)給出的A2= 0.009就是其中一例,截面上分布著大量雜亂無(wú)規(guī)則的點(diǎn),軌道處于混沌運(yùn)動(dòng).當(dāng)A20.1時(shí),軌道逃逸.有輻射無(wú)扁率(A2=0)時(shí),當(dāng)q0.68時(shí),軌道有序.當(dāng)q <0.68時(shí),軌道逃逸.與圖1的擬周期軌道相比,圖2中較小的輻射作用就可以使得混沌軌道逃逸.圖2 (c)給出的是q=0.9,A2=0的結(jié)果.該結(jié)果與同時(shí)含有輻射和扁率的結(jié)果(前提是q0.68,0

4.1.2 Lyapunov指數(shù)

為驗(yàn)證上述結(jié)果,我們計(jì)算另外一個(gè)混沌指標(biāo)Lyapunov指數(shù).Lyapunov指數(shù)通過(guò)衡量?jī)舌徑壍离S時(shí)間平均分離比來(lái)區(qū)分周期(擬周期)運(yùn)動(dòng)和混沌運(yùn)動(dòng),對(duì)保守系統(tǒng)和耗散系統(tǒng)均適用,而且對(duì)系統(tǒng)的維數(shù)不作限制.計(jì)算方法主要有變分法[20–22](同時(shí)積分變分方程和運(yùn)動(dòng)方程)和兩粒子法[23?24](積分運(yùn)動(dòng)方程兩次),后者效率稍高,本文采用兩粒子法.一個(gè)n′維動(dòng)力系統(tǒng)包含n′個(gè)Lyapunov指數(shù),為計(jì)算方便,一般只計(jì)算最大Lyapunov指數(shù)

ξ′(0)和ξ′(t)分別為兩鄰近軌道初始和時(shí)刻t時(shí)的距離.變分法中,兩鄰近軌道的初始變分距離為單位1.而兩粒子法情況下初始時(shí)刻的距離必須適當(dāng)選擇[23],一般在雙精度情況下取為10?8.對(duì)有界軌道,若系統(tǒng)混沌,則λ >0; 若系統(tǒng)有序,則λ= 0.用數(shù)值方法計(jì)算Lyapunov指數(shù)時(shí),理論上時(shí)間應(yīng)為無(wú)窮,本文最終的積分時(shí)間為105.仍然考慮兩條不同的軌道,選取步長(zhǎng)τ= 0.01,軌道初值和相關(guān)參數(shù)與Poincar′e截面中保持一致.每一步積分后,均實(shí)施一次重整化.兩條軌道的Lyapunov指數(shù)如圖3所示,圖中我們將Lyapunov指數(shù)和積分時(shí)間都取了對(duì)數(shù).圖3 (a)給出軌道1的運(yùn)動(dòng)情況.容易看出q=1,A2=0.009與其他3種情況完全不同(該軌道是混沌軌道).無(wú)輻射和扁狀攝動(dòng)情況下,該軌道為擬周期軌道; 無(wú)輻射但扁狀攝動(dòng)不太強(qiáng)時(shí),軌道1混沌; 無(wú)扁狀攝動(dòng)但輻射不太強(qiáng)時(shí),軌道1有序; 既有輻射又有扁狀攝動(dòng)時(shí),軌道1有序.該結(jié)論與圖1的結(jié)果一致.圖3 (b)給出的是原本為混沌的軌道2加入輻射因子和扁狀攝動(dòng)后的結(jié)果.無(wú)輻射但扁狀攝動(dòng)不太強(qiáng)時(shí),軌道2混沌; 無(wú)扁狀攝動(dòng)但輻射不太強(qiáng)時(shí),軌道2有序; 既有輻射又有扁狀攝動(dòng)時(shí),軌道2有序.該結(jié)論與圖2的結(jié)果一致.

圖1 軌道1的Poincar′e截面圖.(a) q = 1, A2 = 0; (b) q = 1, A2 = 0.009; (c) q = 0.9, A2 = 0; (d) q = 0.8,A2 = 0.005Fig.1 Poincar′e sections for the orbit 1.(a) q = 1, A2 = 0; (b) q = 1, A2 = 0.009; (c) q = 0.9, A2 = 0; (d)q = 0.8, A2 = 0.005

計(jì)算Lyapunov指數(shù)時(shí)要求時(shí)間盡量取得足夠大,這樣才能得到一個(gè)收斂于穩(wěn)定數(shù)值的值,如果在給定的時(shí)間段內(nèi)不采用重整化、不對(duì)時(shí)間作平均化處理,可以得到另外一個(gè)隨時(shí)間變化的指標(biāo)[23]—快速Lyapunov指數(shù)(FLI).快速Lyapunov指數(shù)也是一種方便快捷的區(qū)分有序和混沌的指標(biāo).有界軌道中,切向量長(zhǎng)度隨時(shí)間演化的速度是完全不同的.若兩鄰近軌道呈指數(shù)式的偏離,則該軌道是混沌的; 若呈代數(shù)式變大,則是有序軌道.與Lyapunov指數(shù)相比,快速Lyapunov指數(shù)只需要在較短的積分時(shí)間內(nèi)就可以區(qū)分有序和混沌軌道.該模型用快速Lyapunov指數(shù)計(jì)算的結(jié)果與前面給出的Lyapunov指數(shù)結(jié)果非常相似,此處不再重復(fù)列出.

圖2 軌道2的Poincar截面圖.(a) q = 1, A2 = 0; (b) q = 1, A2 = 0.009; (c) q = 0.9, A2 = 0; (d) q = 0.8,A2 = 0.005Fig.2 Poincar sections for the orbit 2.(a) q = 1, A2 = 0; (b) q = 1, A2 = 0.009; (c) q = 0.9, A2 = 0; (d)q = 0.8, A2 = 0.005

圖3 兩條軌道的Lyapunov指數(shù)(λ)圖.(a)軌道1,(b)軌道2Fig.3 Lyapunov exponents λ for the two orbits.(a)The orbit 1; (b)The orbit 2

4.2 多條軌道穩(wěn)定性研究

以上只考慮了兩條特殊軌道的長(zhǎng)期演化情況,此處進(jìn)行多樣本數(shù)值模擬.考慮到實(shí)際情況下伴星的扁率都很小,大范圍討論沒(méi)有必要,令伴星的扁狀攝動(dòng)A2= 0.009.關(guān)于輻射因子,因其與第3體的面質(zhì)比有關(guān),取值可以在較大范圍內(nèi)變化,是一個(gè)可以調(diào)節(jié)的物理量,所以我們著重討論輻射因子對(duì)軌道穩(wěn)定性的影響.首先考慮一組初值條件:Cj= 3.06,x0= 0.29,y0在0–0.51每間隔0.01變化,總共產(chǎn)生了52條軌道.時(shí)間步長(zhǎng)和相關(guān)參數(shù)與前面計(jì)算過(guò)程保持一致.圖4分別給出了52條軌道的Lyapunov指數(shù)和快速Lyapunov指標(biāo)隨y0的演化.需要強(qiáng)調(diào)的是,圖中已經(jīng)剔除了逃逸的軌道.從圖4(a)–(b)中明顯可以看出q= 1,A2= 0.009的模擬結(jié)果與無(wú)輻射無(wú)扁率的結(jié)果多數(shù)情況下保持一致,但也不完全一致.扁狀攝動(dòng)的存在使得部分軌道從有序變?yōu)榛煦?輻射單獨(dú)存在與輻射和扁狀攝動(dòng)同時(shí)存在的結(jié)果基本一致.在經(jīng)典平面圓形限制性三體問(wèn)題中的混沌軌道如果只考慮輻射或者輻射與扁率同時(shí)考慮時(shí)則變成了擬周期軌道,這說(shuō)明輻射和弱扁率項(xiàng)的存在可以增加小天體有序運(yùn)動(dòng)的機(jī)會(huì).

圖4 大樣本軌道的Lyapunov指數(shù)λ (a)和快速Lyapunov指數(shù)FLI (b)Fig.4 Lyapunov exponents λ (a)and fast Lyapunov exponents FLI (b)for the large samples of orbits

除此之外,我們?cè)儆懻摿硗庖唤M初始條件:Cj=3.06,y0=0,x0在0–1變化,結(jié)果見(jiàn)表1所示.(1)當(dāng)x00.04時(shí),第3體離主星太近,無(wú)論有無(wú)輻射和扁率作用,所有軌道均逃逸.(2)當(dāng)x0= 0.05時(shí),原軌道為混沌軌道,輻射的作用可以使軌道逃逸(0q <0.1),也可以使軌道處于混沌(0.1q1),但只要輻射因子項(xiàng)存在,這些軌道都不穩(wěn)定.然而扁率的存在對(duì)軌道的穩(wěn)定性影響卻不大.(3)當(dāng)0.05x00.49時(shí),輻射因子對(duì)軌道穩(wěn)定性的改變比較大.無(wú)論初始軌道是有序還是混沌軌道,較小的輻射因子值都可以使得軌道在短時(shí)間內(nèi)逃逸或瓦解.隨著輻射因子的逐漸增加,短期內(nèi)沒(méi)有逃逸的軌道在較長(zhǎng)的積分時(shí)間內(nèi)都會(huì)處于混沌狀態(tài).當(dāng)輻射因子取值達(dá)到一定程度時(shí),原本處于混沌狀態(tài)的軌道變得有序.需要指出的是,當(dāng)x0= 0.2和x0= 0.4時(shí),輻射因子項(xiàng)可以讓原本有序的軌道逃逸或混沌,也可以繼續(xù)保持有序狀態(tài).其變化規(guī)律與上述其他軌道相似,輻射因子值越大,軌道有序運(yùn)動(dòng)的機(jī)會(huì)越多.扁率對(duì)上述軌道的影響不大,但也有少量軌道除外,比如第3體非?？拷魈祗w但不至于使得軌道逃逸時(shí)(x00.1),扁率可以使得混沌軌道變?yōu)橛行?太陽(yáng)系內(nèi)很難找到與以上區(qū)域?qū)?yīng)的樣本(實(shí)際天體或人造天體),所以上述結(jié)論只存在理論上的參考價(jià)值.(4)當(dāng)x00.5時(shí),輻射因子對(duì)軌道穩(wěn)定性的影響與(1)、(2)和(3)不一樣.輻射因子的作用使軌道要么逃逸要么有序.第3體與主星的距離越遠(yuǎn),逃逸區(qū)間的取值越大,有序運(yùn)動(dòng)要求輻射因子的取值也越大.扁率基本上沒(méi)什么影響,僅在x0= 0.6附近的區(qū)域,可以讓部分有序軌道變?yōu)榛煦畿壍?(5)當(dāng)x00.86時(shí),與(1)情況一樣,輻射和扁率存在與否都不會(huì)產(chǎn)生影響,因?yàn)樵搮^(qū)域不存在穩(wěn)定軌道.從以上多種情況來(lái)看,扁率對(duì)軌道的穩(wěn)定性影響是比較小的,這點(diǎn)與輻射因子完全不同.因?yàn)檩椛湟蜃涌梢允钩跏架壍捞右?也可以使之保持長(zhǎng)期穩(wěn)定有序.

表1 輻射因子q和扁率A2對(duì)軌道穩(wěn)定的影響Table 1 Effects of q and A2 on the stability of the orbits

對(duì)應(yīng)實(shí)際的太陽(yáng)-木星-小行星系統(tǒng),文獻(xiàn)[1]中取輻射因子q= 0.8,扁率A2= 0.005.令Cj= 3.06,結(jié)果如圖5所示.在無(wú)輻射和扁率時(shí),小行星的運(yùn)動(dòng)多數(shù)表現(xiàn)出混沌現(xiàn)象(x>0.5).只考慮扁率時(shí),混沌區(qū)域的范圍有所增加,積分的時(shí)間越長(zhǎng),混沌的區(qū)域還會(huì)擴(kuò)大,逃逸軌道的數(shù)量也越來(lái)越多.只有輻射作用時(shí),不存在混沌區(qū)域,所有軌道都是周期軌道或擬周期軌道.但如果輻射因子q <0.8,逃逸軌道的數(shù)量會(huì)逐漸增多.扁率和輻射同時(shí)作用與輻射項(xiàng)單獨(dú)存在的結(jié)果一致,雖然理論上扁率也有很多可能的取值,但實(shí)際情況下扁率都很小.

總之,無(wú)論是單條軌道演化還是多條軌道的數(shù)值實(shí)驗(yàn)結(jié)果均表明: 在軌道不至于逃逸的情況下,輻射和扁率的共同作用使得相空間中出現(xiàn)有序運(yùn)動(dòng)的幾率增加了,小天體的運(yùn)動(dòng)可以長(zhǎng)期處于穩(wěn)定狀態(tài).但需要強(qiáng)調(diào)的是,引力作用和主星的輻射作用才是主要因素,伴星的扁球攝動(dòng)相對(duì)主星的輻射和引力作用來(lái)說(shuō)是小量,基本上可以忽略不計(jì).總的來(lái)說(shuō),形狀攝動(dòng)和輻射的存在的確可以改變軌道的運(yùn)動(dòng)狀態(tài),增加部分軌道有序運(yùn)動(dòng)的幾率.這點(diǎn)在計(jì)算太陽(yáng)系行星歷表時(shí)若有所體現(xiàn)的話,則太陽(yáng)系各天體穩(wěn)定的時(shí)標(biāo)可能會(huì)更長(zhǎng).

圖5 太陽(yáng)-木星-小行星系統(tǒng)的Poincar截面圖.(a) q = 1, A2 = 0; (b) q = 1, A2 = 0.005; (c) q = 0.8, A2 = 0;(d) q = 0.8, A2 = 0.005Fig.5 Poincar sections for the Sun-Jupiter-Asteroid system.(a) q = 1, A2 = 0; (b) q = 1, A2 = 0.005;(c) q = 0.8, A2 = 0; (d) q = 0.8, A2 = 0.005

5 結(jié)論

本文在經(jīng)典平面圓型限制性三體問(wèn)題中加入了恒星輻射和伴星的扁球攝動(dòng)因素,使該模型更具有實(shí)際物理意義.本文以保流形的速度因子方法為基本工具,研究了該模型的動(dòng)力學(xué)問(wèn)題.數(shù)值模擬結(jié)果表明單純的形狀攝動(dòng)或單純的輻射作用都會(huì)改變軌道的運(yùn)動(dòng)狀態(tài),使得有序運(yùn)動(dòng)變得無(wú)序甚至逃逸,也可以使得混沌運(yùn)動(dòng)變得有序.本文的數(shù)值研究結(jié)果表明: 當(dāng)小天體的運(yùn)動(dòng)沒(méi)有受到強(qiáng)烈的輻射作用或扁狀攝動(dòng)時(shí),特別是輻射因子的取值不能低于一定的限度,一般情況下小天體的運(yùn)動(dòng)不會(huì)逃逸到無(wú)窮遠(yuǎn)處.此時(shí)輻射和扁率的共同作用可以增加系統(tǒng)的可積區(qū)域,增加小天體有序運(yùn)動(dòng)的機(jī)會(huì).對(duì)于實(shí)際的太陽(yáng)-木星-小天體系統(tǒng),數(shù)值模擬證實(shí)了上述理論結(jié)果.因此在計(jì)算太陽(yáng)系行星歷表時(shí),建議考慮太陽(yáng)輻射的影響.