二體問題下中心力場滿足二次平方反比之斥力問題

好久沒更新部落格了,今日留個爪印。

最近想徹底解決掉遺漏許久的彗星塵埃粒子在力場中的問題,如今所寫出的IDL代碼已經可以完美地解決了\beta = \frac{F_{rad}}{F_{grav}}< 1的情況。如下圖所示,模擬C/2011 L4 (PanSTARRS)彗星在今年4月1日的部分塵埃尾巴情況。

Synchrone and syndyne grids for the comet
Synchrone and syndyne grids for the comet

解決這個問題的本質其實就是不斷地計算不同塵埃粒子的Kepler’s Equation。由於粒子大小不同,釋放時間也不同,其mean anomaly和偏心率也會不同。計算方法參考的是Probstein和Finson於1968年所發表的一片論文《A theory of dust comets. I. Model and equations》。不過論文僅針對了軌道是拋物線運行的彗星的情況。但只要按照同樣的思路,橢圓軌道、雙曲線軌道的彗星情況同樣可以迎刃而解。然而論文中只討論了\beta < 1,對於大於及等於1之情況卻未有照顧。不等式取等號時對應塵埃粒子將以勻速直線運動方式運動,而大於時則是受到中心斥力。於是我很自然地聯想到,在中心斥力力場作用下,塵埃粒子的運動方程是否仍然有類似於中心吸引力的類似形式。

起初爲了方便了事,想在網上直接尋找答案,然而找了許久,並無收穫,於是決心自己動筆計算以探究竟。以下用極座標進行討論,認爲塵埃質量與太陽質量相比是無窮小量。原點設在中心力場中心。我們不難知道可以推導出以下兩組方程:

\ddot{r} - r\dot{\theta}^2 = \frac{\mu}{r^2}——(1)

r^2\dot{\theta}=h——(2)

其中\mu = \left ( \beta - 1 \right )GM_\odot 。首先先推導軌道類型。令u = \frac{1}{r}代入(1)式,並結合(2)式,可以得到如下一個二次非齊次線性微分方程:

\frac{\mathrm{d}^2 u}{\mathrm {d}\theta^2}+u=-\frac{\mu}{h^2}

可得其通解爲u=\frac{\mu}{h^2}\left [ e\cos \left ( \theta - \omega \right ) - 1\right ],其中e和\omega爲待定常數。將其表達回r,得到:

r=\frac{h^2/\mu}{ e\cos \left ( \theta - \omega \right ) - 1}——(3)

由解析幾何知,這是一條以外焦點爲原點的雙曲線。可記得Rutherford Scattering,被散射粒子正是沿着雙曲線軌道運行的!

我們定義q = r_{min} = \left | a \right |(1+e),a爲對應雙曲線之半長徑,爲負值。剩下的工作是找出r和\theta關於時間t的含數。

將(2)式代入(1)可以得到:

\ddot{r} - \frac{h^2}{r^3}=\frac{\mu}{r^2}——(4)

由於\ddot{r} = \dot{r}\frac{\mathrm{d} \dot{r}}{\mathrm{d} r},故(4)式可以化爲:

\dot{r}\mathrm{d} \dot{r}=(\frac{h^2}{r^3}+\frac{\mu}{r^2})\mathrm{d}r

兩邊積分得到:

\dot{r}^2=-\frac{2\mu}{r}-\frac{h^2}{r^2}+K_1——(5)

其中K_1是一個待定積分常數,可以用初始條件代入求出。方法如下。

對於r=q時,有\dot{r}=0,代入(5)可以解得K_1 = \frac{\mu}{\left | a \right |},即

\dot{r}^2=-\frac{2\mu}{r}-\frac{h^2}{r^2}+\frac{\mu}{\left | a \right |}——(6)

注意到v^2 = \dot{r}^2 + \dot{\theta}^2 r^2 ,對(6)變換下形式,於是立即獲得了一個與活力公式十分相似的一條描述粒子速率的一條公式:

v^2 = \mu (\frac{1}{\left | a \right |}-\frac{2}{r})——(7)

可以發現,只要把(7)式括號內的兩項調整位置,便是「正宗」的活力公式!

我們的任務還未完成,現繼續求r和\theta關於t的表達式。將(6)式變形,可以得到:

\mathrm{d}t = \frac{r\mathrm{d}r}{\sqrt{\frac{\mu}{\left | a \right |}}r^2 - 2\mu r - \mu \left | a \right |(e^2 - 1)}

即:

\sqrt{\frac{\mu}{\left | a \right |}}\mathrm{d}t=\frac{r\mathrm{d}r}{\left | a \right |\sqrt{(r - \left | a \right |)^2 - a^2 e^2}}

我們令r - \left | a \right |=\left | a \right |e\cosh{F},則有:

\sqrt{\frac{\mu}{\left | a \right |^3}}\mathrm{d}t = (e\cosh{F}+1)\mathrm{d}F

兩邊同時積分,得到:

\sqrt{\frac{\mu}{\left | a \right |^3}}t = e\sinh{F}+F+K_2

廣義Kepler’s Equation雛形初現!!上式中K_2爲一積分常數,用同樣的方式,考慮r = q時,選擇t=0爲過近日點之時,則可得到K_2 = 0 。於是廣義Kepler’s Equation爲

\sqrt{\frac{\mu}{\left | a \right |^3}}t = e\sinh{F}+F——(8)

可見,在斥力下的方程與引力下的方程形式上非常接近。仿照真正意義上的Kepler’s Equation,我們且將F亦稱爲eccentric anomaly(偏近點角)。這是個超越方程,只要將F解出,根據

r = \left | a \right |(1 + e\cosh{F})——(9)

則可求得r。r既得,\theta亦可根據(3)式求出。不過此處我們想看看\theta與F之間的關係爲何。爲簡單見,令 f = \theta - \omega

將(9)式代入(3)式,我們可以得到\cos{f} = \frac{e + \cosh{F}}{1 + e \cosh{F}}。利用半角公式2 \sin^2 {\frac{f}{2}} = 1 - \cos{f},得:

\sin{\frac{f}{2}} = \sqrt{\frac{(e-1)\sinh^2 {\frac{F}{2}}}{1+e \cosh{F}}}

類似地,利用2 \cos^2 {\frac{f}{2}} = 1 + \cos{f},可得:

\cos{\frac{f}{2}} = \sqrt{\frac{(e+1)\cosh^2 {\frac{F}{2}}}{1+e \cosh{F}}}

兩式相除,最後得到:

\tan{\frac{f}{2}} = \sqrt{\frac{e-1}{e+1}}\tanh{\frac{F}{2}}——(10)

至此,問題已經完美解決了!我們可以得到的結論是,中心力滿足二次平方反比條件的中心斥力情況下,被散射物體的運動軌跡方程與對應吸引力場下的雙曲線軌道運動方程十分相似,我們只需改變幾個地方的符號即可獲得。

Advertisements
二體問題下中心力場滿足二次平方反比之斥力問題