行星形变理论（2）

这是行星形变理论的第二部分。第一部分在这里。

啊……

总感觉这篇文章很难写清楚，主要是那个\(k\)值的问题。为了得到真正的形变量，必须先假设行星自身势是质点势来简化计算，等到把各种与形变\(C\)有关的量都算得差不多了，才能考虑自身势变化，并证明真正的形变量是原本的\(k\)倍。结果原来的那些符号现在到底指的是乘\(k\)之前的量还是乘\(k\)之后的量呢？是不是有点混乱……

为了之后讨论简便，从自转和潮汐力理论计算的扁率\(f_\omega\)，\(f_m\)和形变矩阵\(C\)都定义为乘\(k\)之前的量；而形变的变换矩阵\(T\)和势能变化\(\Delta\phi\)，引力场变化\(\Delta\vec g\)等都会定义成乘\(k\)之后真实的量。为了避免混乱，开始第二部分之前首先重新捋清一遍符号定义。

前情提要

考虑半径为\(R\)的球对称行星，其内部距中心为\(r\)处的密度为\(\rho(r)\)，行星的总质量为
\[M=4\pi\int_0^R r^2\rho(r)dr\]

如果该行星自转角速度为\(\vec\omega\)，并受到多个其他天体的潮汐力，其中质量为\(m_i\)的天体与其相对位置为\(\vec r_i\)\((i=1,2,\cdots)\)，则可以定义形变矩阵
\[C=C(\vec\omega,f_\omega)+k_t\sum_i C(\vec r_i,f_{m_i})\]

其中
\[f_\omega=\frac{\omega^2R^3}{2GM}\\
f_{m_i}=-\frac{3Gm_iR^3}{2GMr_i^3}\\
C(\vec r,f)=\left(\frac{I}{3}-\frac{\vec r\vec r}{r^2}\right)f
\]

\(k_t\)是行星对潮汐的形变响应因子，为一经验数值，描述行星的刚性。\(f_\omega\)与\(f_{m_i}\)分别为行星因自转与来自\(m_i\)的潮汐力而产生的理论扁率，一般来说都是远小于\(1\)的小量。

行星在自转角速度及其他天体潮汐力的影响下会产生形变。行星上原本位置为\(\vec r\)的点经过形变后会到达\(T\vec r\)，其中\(T=I+kC\)。无量纲量\(k\)是一个与行星内部物质分布有关的量，定义为
\[k=\frac{1}{1-AA_r}\\
A=\frac{4\pi}{MR^2}\int_0^R r^4\rho(r)dr\]

其中\(A_r\)为一个描述行星内部形变的不均匀性的经验系数。由于\(C\)是一个零迹对称矩阵，而且是小量，所以\(T\)的行列式几乎是\(1\)，也即这形变过程是不改变行星体积和密度的。

实际中，如果实测的行星扁率\(f\)可以认为完全归因于其自转，则可由此计算\(k\)值
\[k=f/f_\omega=\frac{2GMf}{\omega^2R^3}\]

形变后，行星对其外部一点\(\vec r\)处的引力势为\(\phi(\vec r)=\phi_0(\vec r)+\Delta\phi(\vec r)\)，其中
\[\phi_0(\vec r)=-GM/r\\
\Delta\phi(\vec r)=-\frac{GM(k-1)R^2}{r^5}\vec r\cdot C\vec r\]

分别为行星形变前的质点势以及因形变而产生的附加势。

第四章 形变行星的转动惯量

之后将要讨论行星在因形变产生的额外力矩下如何运动，所以首先得计算转动惯量\(J\)。相比起引力势，转动惯量的计算是比较简单的。直接写出转动惯量张量：
\[J=\int_{|\vec r’|\lt R}(|T\vec r’|^2I-(T\vec r’)(T\vec r’))\rho(|\vec r’|)dV’\]

注意到\(|T\vec r’|^2I\)可以写成\((I:(T\vec r’)(T\vec r’))I\)的形式，故\(J\)可以写成
\[J=(I:J_1)I-J_1\]

其中
\[J_1=\int_{|\vec r’|\lt R}(T\vec r’)(T\vec r’)\rho(|\vec r’|)dV’\]

继续计算之前先化简其中的被积张量：
\[\begin{aligned}
&(T\vec r’)(T\vec r’)\\
=&(\vec r’+kC\vec r’)(\vec r’+kC\vec r’))\\
=&\vec r’\vec r’+k(C\vec r’)\vec r’+k\vec r'(C\vec r’)+k^2(C\vec r’)(C\vec r’)
\end{aligned}
\]

上式中最后一项为\(C\)的高阶小量，可略去。中间含\(C\)的两项可写成张量\(C\)与\(\vec r’\vec r’\)点积的形式，并注意到\(C\)为对称矩阵，故有
\[(T\vec r’)(T\vec r’)=\vec r’\vec r’+kC\cdot\vec r’\vec r’+k\vec r’\vec r’\cdot C\]

故\(J_1\)又可以写成
\[J_1=J_2+kC\cdot J_2+kJ_2\cdot C\]

其中
\[J_2=\int_{|\vec r’|\lt R}\vec r’\vec r’\rho(|\vec r’|)dV’\]

在球坐标系中将\(\vec r’\vec r’\)写出
\[\vec r’\vec r’=r^2\begin{pmatrix}
\cos^2\phi\sin^2\theta & \cos\phi\sin\phi\sin^2\theta & \cos\phi\cos\theta\sin\theta\\
\cos\phi\sin\phi\sin^2\theta & \sin^2\phi\sin^2\theta & \cos\theta\sin\phi\sin\theta\\
\cos\phi\cos\theta\sin\theta & \cos\theta\sin\phi\sin\theta & \cos^2\theta
\end{pmatrix}
\]

可以继续计算\(J_2\)
\[\begin{aligned}
J_2=&\int_0^R r^2\rho(r)dr\int_0^\pi \sin\theta d\theta\int_0^{2\pi}\vec r’\vec r’d\phi\\
=&\int_0^R r^4\rho(r)dr\int_0^\pi \pi\sin\theta d\theta\begin{pmatrix}
\sin^2\theta & 0 & 0\\
0 & \sin^2\theta & 0\\
0 & 0 & 2\cos^2\theta
\end{pmatrix}\\
=&\frac{4\pi}{3}\int_0^R r^4\rho(r)dr\begin{pmatrix}
1 & 0 & 0\\
0 & 1 & 0\\
0 & 0 & 1
\end{pmatrix}\\
=&\frac{AMR^2I}{3}&
\end{aligned}\]

可以看到这里又出现了无量纲量\(A\)。注意这里的\(A\)是计算行星转动惯量时出现的，和行星内部形变无关，故无需用\(A_r\)做修正。将\(J_2\)代入\(J_1\)及转动惯量\(J\)的表达式并化简，有
\[J_1=\frac{AMR^2}{3}(I+2kC)\\
J=\frac{2AMR^2}{3}(I-kC)
\]

这就是形变行星的转动惯量张量。

一般来说，计算行星的转动惯量需要测量其内部物质分布。同样半径和质量的行星，物质分布越聚集于核心，转动惯量就越低。但本章的结论表明只需求得\(A\)值就可以简单地计算出转动惯量。假设行星形变均匀，\(A_r\approx 1\)，\(A\)值可以通过测量行星的实际扁率\(f\)和自转角速度\(\vec\omega\)测出，这要比测量行星内部物质分布容易得多。

第五章 形变产生的行星间额外的力与力矩

之前已经计算出形变行星的引力势\(\phi\)，从引力势可以直接求得引力场：
\[\vec g=\vec g_0+\Delta\vec g\]

其中\(\vec g_0=-\nabla\phi_0\)为行星的质点势产生的引力场，\(\Delta\vec g=-\nabla\Delta\phi\)为形变产生的附加引力场。计算结果为
\[\vec g_0(\vec r)=-\frac{GM\vec r}{r^3}\\
\Delta\vec g(\vec r)=-\frac{GM(k-1)R^2}{r^7}(5(\vec r\cdot C\vec r)\vec r-2r^2C\vec r)
\]

这个引力场会作用于所有其他天体。若有质量为\(m_i\)的天体与形变行星的相对位置为\(\vec r_i\)，则天体\(m_i\)受到引力\(\vec F_i=m_i\vec g(\vec r_i)\)，形变行星自身受到反作用力\(-\vec F_i\)。由于\(\vec F_i\)因形变产生的\(\Delta\vec g\)而不再与两天体连线\(\vec r_i\)平行，故这一对引力会引起一个非零的力矩\(\vec r_i\times\vec F_i\)。由力矩平衡，形变天体本身必然受到相反的力矩。其角动量\(\vec L\)将发生变化
\[\begin{aligned}
\frac{d\vec L}{dt}=&-\sum_i\vec r_i\times\vec F_i\\
=&-\sum_i m_i\vec r_i\times\Delta\vec g(\vec r_i)
\end{aligned}\]

第六章 形变延时与潮汐摩擦

上一章讨论了行星角动量\(\vec L\)会逐渐发生变化，这一变化也许可以解释行星的进动与章动，但并不能解释潮汐锁定现象。原因是显然的：整个讨论只考虑了万有引力，它本身是保守力，在它的作用下能量一定是守恒的。而潮汐锁定过程中必然出现能量的损耗，也就是潮汐摩擦。这一损耗从何而来？

答案其实非常简单，就是行星的顽抗力：行星并不能随着天体间的运动随时呈现完美的平衡形变，而是有坚持其过去的形变的“粘滞性”。我们从现在开始用\(C(t)\)表示\(t\)时刻行星的形变矩阵，这样本文之前根据\(C\)计算出来的关于\(\phi\)，\(J\)，\(\vec g\)，\(d\vec L/dt\)的各个公式仍均成立。而用\(C_0\)表示从\(t\)时刻的行星自转角速度和其他天体的潮汐力计算出来的理论平衡形变矩阵。

\(C_0\)是随着天体的运动而不断变化的，在行星没有粘滞的情况下\(C\)会等于\(C_0\)，但考虑行星粘滞，可以给出如下\(C\)随时间变化的经验公式
\[\frac{dC}{dt}=\frac{C_0-C}{\tau}\]

其中\(\tau\)是行星形变滞后于天体运动的时标，这也是一个经验量。对常见天体，比如地球和月球，大概在几个小时左右。这一公式当然是在行星的自转坐标系中成立的，想象一下粘滞性\(\tau\rightarrow\infty\)，行星的自转就会拖着形变\(C\)转动，所以在自转参考系中\(dC/dt=0\)。

显然在\(t\)时刻，行星处于平衡形变\(C_0\)时具有最低势能，而行星的实际形变却为\(C\)，并且处于向\(C_0\)变化之中，这个过程当然是损耗能量的。能量的损耗速率也可依此计算。首先写出行星总势能的表达式
\[E=\int_{|\vec r’|\lt R}\left(\frac{\phi(T\vec r’)}{2}+\phi_{ext}(T\vec r’)\right)\rho(|\vec r’|)dV’\]

算出来后再计算\(dE/dt\)即可。

……

………………

→_→

别看我，看我也没用。我不会算这个\(E\)。

不开玩笑，我目前是真的算不出来。不过我简单说一下难度都在哪。首先这个\(\phi(T\vec r’)\)是行星内部的势，而我们之前只算了行星外部的势，所以首先这个是要重算的。另外这个\(\phi_{ext}\)，根据第三章的某个公式我们可以将行星表面上的\(\phi_{ext}\)用\(\Delta\phi\)表示，可它在整个行星内部是如何变化的？当然由于它是小量，可以证明它满足\(\phi_{ext}(T\vec r’)=(r’/R)^2\phi_{ext}(RT\vec r’/r’)\)。总之都比较复杂。算出来还没完，要算\(E\)还得把它再积一遍，我暂时是积不动了……

所以我直接猜结果为
\[\frac{dE}{dt}=-\frac{8GM^2(k-1)(C_0-C)^2}{75R\tau}\]

其中\((C_0-C)^2=(C_0-C):(C_0-C)\)。

我猜的不一定对，但量级肯定是正确的，而形式我有\(80\%\)的把握是正确的。所以它顶多与正确答案差一个常数项系数，这系数可能还和\(k\)有关。以后如果有一天我严格算出来了，我再更新这一章……

————————更新（2020.9.7）

我想到怎么优化这个延迟形变模型了。

前面的模型有几个问题。第一是它的数值积分有困难，因为它是一个非保守方程（指数衰减），所以数值积分时，时间步长不能为负值（不然会指数发散），这样就不能使用ForestRuth等更好的数值积分方法了；第二是很难根据这个模型算出潮汐耗散功率的公式，之前就一直没算出来。

但考虑到这个耗散模型本身就是经验的，那不妨再对它进行一次近似，在一个时刻直接计算延迟之后的潮汐形变矩阵\(C’\)，而不是计算\(C_0\)再去数值积分\(C\)。

假设质量为\(M\)，半径为\(R\)的行星的自转角速度为\(\vec\omega\)，另有质量为\(m\)的天体相对\(M\)以速度\(\vec v=\vec v_m-\vec v_M\)运动，则可以计算它们的相对角速度为
\[\Delta\vec\omega=\frac{\vec r\times (\vec\omega\times\vec r-\vec v)}{r^2}\]
其中\(\vec r=\vec r_m-\vec r_M\)为它们的相对位移。

以\(M\)为中心，\(\vec r\)为\(x\)轴正向，\(\Delta\vec\omega\)为\(z\)轴正向建系，并认为\(\Delta\vec\omega\)和\(r\)在\(\tau\)的时标内近似不变，则在\(M\)的自转系中\(m\)的相对运动可以近似为稳恒的圆周运动，半径为\(r\)，角频率为\(-\Delta\omega\hat z\)。此时再根据\(C_0\)的公式和\(C\)满足的延迟演化方程，就能直接解析积分出\(C\)来。

积分过程省略，总之这样解析积分出来的延迟矩阵表达式为
\[C'(\vec r,f;\Delta\vec\omega,\tau)=\left(\left(\frac{1}{3}-z_0^2\right)I-\frac{\vec r’\vec r’}{r^2}+z_0^2\frac{\Delta\vec\omega\Delta\vec\omega}{\Delta\omega^2}\right)f\]
其中
\[
\vec r’=(x_0\hat x+y_0\hat y)r=x_0\vec r+y_0\frac{\Delta\vec\omega\times\vec r}{\Delta\omega}\\
x_0=\sqrt{\frac{1+\delta}{2\delta^2}},\quad
y_0=\sqrt{\frac{2\alpha^2}{(1+\delta)\delta^2}},\quad
z_0=\sqrt{\frac{2\alpha^2}{(1+\delta)\delta}}\\
\delta=\sqrt{1+4\alpha^2},\quad
\alpha=\Delta\omega\tau
\]

到此为止，延迟形变矩阵已经算完了。可以使用这个\(C’\)代替原先的\(C\)，结果就自然计入了潮汐耗散。最后从这个延迟矩阵算一下潮汐耗散功率
\[W=\vec\omega\cdot(-\vec r\times m\Delta\vec g)+\vec v\cdot m\Delta\vec g\]

计算过程留作习题（懒得写了）。总之最后结果是可以化简掉\(\vec\omega\)和\(\vec v\)，只用\(\Delta\vec\omega\)表示的：
\[W=-\frac{3Gm^2k_t(k-1)R^5\Delta\omega^2\tau}{(1+4\Delta\omega^2\tau^2)r^6}\]

————————更新完毕

小结 另一些数值计算

写到这里本文基本上结束了。最后再用一些数值计算举例吧。

首先根据本文的结论理论计算一下地轴的进动周期。地球的黄赤交角，也就是地球自转轴与黄道坐标系\(z\)轴的夹角为\(\epsilon=23^\circ26’\)。在本节计算中选择的参考系以地球中心为原点，黄道平面为\(xy\)平面，并使地球北极位于\(+x,+z\)平面中。这样地球自转角速度可以表示为\(\vec\omega=(\omega\sin\epsilon,0,\omega\cos\epsilon)^T\)。

仍然假设地球的形变几乎全部由自转引起，可以计算地球的形变矩阵与角动量
\[C=\left(\frac{I}{3}-\frac{\vec\omega\vec\omega}{\omega^2}\right)f_\omega\\
\begin{aligned}
\vec L=&J\vec\omega\\
=&\frac{2AMR^2}{3}\left(I-kC\right)\vec\omega\\
=&\frac{2AMR^2}{3}\left(1-k\left(\frac{1}{3}-1\right)f_\omega\right)\vec\omega\\
=&\frac{2AMR^2}{3}\left(1+\frac{2f}{3}\right)\vec\omega\\
=&J_\omega\vec\omega
\end{aligned}\]

其中\(f=kf_\omega\)是地球的实测扁率。代入实际数值，可得地球相对其自转轴的转动惯量为
\[J_\omega=0.331539MR^2=8.0368\times 10^{37}\mathrm{kg\cdot m^2}\]

地球自转角速度的变化主要有两个来源，月球和太阳。其中月球的轨道平面和黄道的夹角很小，而且比较稳定（相对于它和地球赤道的夹角），目前在\(5.145^\circ\)左右。故为简便起见，将月球运动近似看作是黄道平面上的圆周运动。同样，太阳相对地球的运动也是黄道平面上的圆周运动。

现在计算这样的圆周运动会如何影响地球的角动量。假设质量为\(m\)的天体以角速度\(\Omega\)在\(xy\)平面上做半径为\(r\)的圆周运动，周期\(T=2\pi/\Omega\)。其位置随时间的变化为\(\vec r(t)=(r\cos(\Omega t),r\sin(\Omega t),0)^T\)。则在一个周期内，它对地球的平均力矩为
\[\begin{aligned}
\vec\tau_m=&-\frac{1}{T}\int_0^T m\vec r(t)\times\Delta\vec g(\vec r(t))dt\\
=&-\frac{2GMm(k-1)R^2}{r^5T}\int_0^T \vec r(t)\times C\vec r(t)dt\\
=&-\frac{2GMm(k-1)R^2f_\omega}{r^5T}\int_0^T \vec r(t)\times\left(\frac{\vec r(t)}{3}-\frac{\vec\omega\cdot\vec r(t)}{\omega^2}\vec\omega\right)dt\\
=&\frac{m(k-1)R^5}{r^5T}\int_0^T (\vec r(t)\cdot\vec\omega)(\vec r(t)\times\vec\omega)dt\\
\end{aligned}\]

带入\(\vec r(t)\)与\(\vec\omega\)，写出被积向量的分量形式
\[(\vec r(t)\cdot\vec\omega)(\vec r(t)\times\vec\omega)=\omega^2r^2\sin\epsilon\cos\Omega t\begin{pmatrix}
\cos\epsilon\sin\Omega t\\
-\cos\epsilon\cos\Omega t\\
-\sin\epsilon\sin\Omega t
\end{pmatrix}\]

继续计算平均力矩
\[\begin{aligned}
\vec\tau_m=&\frac{m(k-1)R^5}{r^5T}\int_0^T (\vec r(t)\cdot\vec\omega)(\vec r(t)\times\vec\omega)dt\\
=&-\frac{m(k-1)\omega^2R^5}{2r^3}\begin{pmatrix}0\\1\\0\end{pmatrix}\cos\epsilon\sin\epsilon\\
=&-\frac{m(k-1)\omega^2R^5\hat y}{2r^3}\cos\epsilon\sin\epsilon
\end{aligned}
\]

注意到\(\vec\tau_m\)可以写成\(\Omega_p\hat z\times\vec L\)的形式，其中
\[\Omega_p=-\frac{m(k-1)\omega R^5}{2 r^3J_\omega}\cos\epsilon\]

而力矩是角动量随时间的变化，即\(\vec\tau_m=d\vec L/dt\)。对照可得，\(\Omega_p\hat z\)正是角动量\(\vec L\)的进动角速度。

代入月球重力参数\(Gm_m=4.9048695\times 10^{12}\mathrm{m^3/s^2}\)，月球轨道平均半径\(r_m=384400\mathrm {km}\)，可得月球引起的进动角速度
\[\Omega_{pm}=-5.32177\times 10^{-12}\mathrm{rad/s}\]

再代入太阳重力参数\(Gm_s=1.32712440018\times 10^{20}\mathrm{m^3/s^2}\)和地球轨道平均半径\(r_e=149597870.7\mathrm{km}\)，可得太阳引起的进动角速度
\[\Omega_{ps}=-2.44295\times 10^{-12}\mathrm{rad/s}\]

于是可得地球自转轴进动的周期
\[T_p=\frac{2\pi}{|\Omega_{pm}+\Omega_{ps}|}=25642\mathrm{yr}\]

在网上可以查到地轴进动的实际周期是\(25772\mathrm{yr}\)，可见我们的结果还是比较准确的，误差为\(0.50\%\)。主要误差来源应该是对月球轨道做的近似。

————————更新（2020.9.2）

第一部分引入了\(A_r\)解释了从行星密度径向分布计算\(k\)值的误差，但是并没有验证引入该参数的合理性。这里使用地轴进动的数值计算结果做一说明。

建立太阳-地球-月球三体运动的数值模型，初值使用实际天体的轨道参数（J2000），进行长度为100年的数值积分。数值积分考虑牛顿引力、广义相对论后牛顿近似、以及本文所描述的行星自转和潮汐形变作用及这些作用的附加力和力矩。当取\(A_r=1\)，即不考虑\(A_r\)的影响时，模拟所得地轴进动周期为\(25202\mathrm{yr}\)，误差为\(2.21\%\)。当取合理的\(A_r\)修正了扁率和物质分布的分歧时，模拟所得地轴进动周期为\(25787\mathrm{yr}\)，误差仅为\(0.06\%\)，大部分误差被消除，剩余的进动速率可能由其他行星引起。据此可以确认\(A_r\)具有实际的意义，因为一处无意义的修正不太可能同时精确修复两处误差。

————————更新完毕

最后提一下地球的潮汐损耗功率吧。由于第六章中的公式是纯粹经验的，计算结果只是估计，所以就不在这里拿公式具体计算了。对于比较合适的\(\tau\)的取值，比如\(10\)分钟，分别考虑太阳和月亮的潮汐，可以估计出地球的潮汐损耗功率在\(-3.8\times 10^{12}\mathrm W\)左右，这个估计和网上找到的数值是同一个数量级的。

终章

啊，终于把这篇东西糊完了。

本文讨论了行星在自转和潮汐力作用下如何形变，以及这形变造成的各种天体力学效果。最后说说我为什么要算这些玩意。

众所周知，有一个著名的天体力学模拟游戏，Universe Sandbox！

我写这篇东西的目标和它差不多，但和它的侧重点不同。Universe Sandbox是为了实时模拟行星运动供人游戏，它为了性能会牺牲精度。而我是为了能尽可能精确地积分一个天体系统的星历表，这样我就能造出一个虚拟的恒星系，并且赋予其真实而精准的历法！

——为了创造一个虚拟而真实的世界。

之后还有第三部分，讨论本文的模型在实现中遇到的问题。