『Python数値計算ノート』ではアフィリエイトプログラムを利用して商品を紹介しています。

万有引力の法則と重力加速度

万有引力

$2$ 個の物体があって、それぞれの質量を $M,\ m$、物体間の距離を $r$ とするとき、それぞれの物体には
 \[F=G\frac{Mm}{r^2}\tag{1}\]
の大きさの引力がはたらきます (万有引力の法則)。
 万有引力の法則の説明図
 
$G$ は 万有引力定数 とよばれる値です。2018 年の CODATA (科学技術委員会) が提供する推奨値は
 \[G=6.67430\times 10^{-11}[\mathrm{m}^{3}\mathrm{kg}^{-1}\mathrm{s}^{-2}]\tag{2}\]
となっており、SciPy 1.10.1 はこの値を採用しています:

# In[1]

import numpy as np
import scipy.constants as const
from scipy.constants import physical_constants as p_const

# 配列要素は小数点以下5桁まで表示
np.set_printoptions(precision=5)

# 万有引力定数の値
print(const.G)
print(const.gravitational_constant)
# 6.6743e-11
# 6.6743e-11

物理定数を収めた辞書 scipy.constants.physical_constants を使うと、値の他に単位や測定の不確かさを取得できます。

# In[2]

# 万有引力定数の値,単位,不確かさ
print(p_const["Newtonian constant of gravitation"])
# (6.6743e-11, 'm^3 kg^-1 s^-2', 1.5e-15)

現代においても万有引力定数の正確な測定は技術的に困難なので、CODATA 推奨値も変遷を続けています(ただし、最新推奨値が真値に近いという保証はありません)。SciPy のバージョンが更新されると、scipy.constants.G の値も変更される可能性があります。
 
(1) をベクトル表式で書き直してみましょう。質量 $M$ の物体 A が $\boldsymbol{r}_0$ の位置に、質量 $m$ の物体 B が $\boldsymbol{r}$ の位置にあるとき、物体 B にはたらく引力の絶対値は
 \[F=G\frac{Mm}{|\boldsymbol{r}_0-\boldsymbol{r}|^2}\tag{3}\]
と表せます (下図)。
 Python 万有引力の法則のベクトル表式
 
力の方向を単位ベクトルで表すと
 \[\frac{\boldsymbol{r}_0-\boldsymbol{r}}{|\boldsymbol{r}_0-\boldsymbol{r}|}\tag{4}\]
となるので、(3) と (4) を掛ければ、ベクトルで表された万有引力の法則が得られます:
 \[\boldsymbol{F}=GMm\frac{\boldsymbol{r}_0-\boldsymbol{r}}{|\boldsymbol{r}_0-\boldsymbol{r}|^3}\tag{5}\]

地球と太陽の間にはたらく万有引力

ベクトル表式の万有引力の法則を Python で実装してみます。

# In[3]

# 万有引力関数
# m0,r0:物体Aの質量(スカラー),位置(ベクトル)
# m,r:物体Bの質量(スカラー),位置(ベクトル)
def gravitation_vector(m0, r0, m, r):
    vec = np.array(r0) - np.array(r)
    return const.G * m0 * m * vec / np.linalg.norm(vec)**3

地球と太陽の質量をそれぞれ $M_e=5.972\times 10^{24}$ kg、$M_s=1.989\times 10^{30}$ kg とし、地球は近日点 (太陽に最も近い場所) にあるとして、その位置をベクトルで
 \[(1.471\times 10^{11}, 0, 0)\]
のように設定します。
Python 地球が近日点の位置にあるときに太陽から受ける引力ベクトル
このとき、地球にはたらく引力の大きさは以下のように計算できます。

# In[4]

# 地球の質量
me = 5.972e24

# 太陽の質量
ms = 1.989e30

# 太陽の位置を原点とする
rs = [0, 0, 0]

# 地球の近日点
r = [1.471e11, 0, 0]

# 地球が近日点にあるときにはたらく引力
fe = gravitation_vector(ms, rs, me, r)

print(fe)
# [-3.66371e+22  0.00000e+00  0.00000e+00]

重力加速度

天体の質量を $M$、天体の影響を受ける物体の質量を $m$ とするとき、運動方程式の右辺に万有引力を与えると、天体の重力によって生じる物体の加速度 (重力加速度) を求めることができます。すなわち、引力のはたらく方向 (物体の中心から天体の中心へ向かう方向) を $x$ 軸としたとき、物体の運動方程式は
 \[m\frac{d^2 x}{dt^2}=G\frac{Mm}{r^2}\tag{6}\]
となります。この式の両辺から $m$ を落とすと加速度の式
 \[\frac{d^2 x}{dt^2}=G\frac{M}{r^2}\tag{7}\]
が得られます。この式の $M$ に地球質量 $M_e=5.972\times 10^24$ kg、$r$ に地球の赤道半径 $R=6.378\times 10^6$ m を与えると、地球の赤道上における重力加速度を求めることができます:
 \[g=G\frac{M_e}{R^2}\tag{8}\]

重力加速度の計算

Python で天体の地表面付近の重力加速度を求める関数を定義して、地球の赤道上における重力加速度を計算してみましょう。

# In[5]

# 天体の重力加速度を求める関数
def g_acceleration(m, r):
    return const.G * m / r**2

# 地球の赤道半径
re = 6.378e6

# 地球の赤道上における重力加速度
g = g_acceleration(me, re)

print("{:.5f} [m/s^2]".format(g))
# 9.79811 [m/s^2]

地表から $h$ m だけ離れたところにある物体にはたらく重力を正確に計算するためには、(9) の $R$ の代わりに $R+h$ を使う必要があります。しかし、$h$ が $R$ に比べて十分小さいときは (9) で近似しても十分な精度が得られます。

遠心力による補正

地球は自転しているので、地表にあるすべての物体には自転軸から離れる方向に遠心力がはたらきます。
Python 遠心力の効果を含めた重力加速度
地球の自転角速度を $\omega$ で表すと、赤道上では遠心力 $f_e=mR\omega^2$ が重力に抗する方向にはたらくので、遠心力効果によって補正された重力加速度は
 \[g_e=g-R\omega^2\tag{9}\]
となります。緯度 $\varphi$ の地点における遠心力による加速度は
 \[f=(R\cos \varphi)\omega^2\tag{10}\]
で与えられます。緯度が大きいほど遠心力の効果は小さくなり、極地で $0$ となります。遠心力効果によって補正された重力加速度は余弦定理を使って、
 \[g’=\sqrt{g^2+f^2-2gf\cos\varphi}\tag{11}\]
によって計算できます。
 
遠心力の効果を含めた重力加速度を求める関数を Python で実装してみます。

# In[6]

# 天体表面の重力加速度を計算する関数

# m:天体の質量(mass)
# r:天体の半径(radius)
# rp:天体の自転速度(rotation period)
# phi:重力加速度を計算する緯度
# degrees:phiを度数単位で指定するか否か

# 戻り値:
# (遠心力を考慮しない重力加速度,
#  遠心力による加速度,
#  遠心力効果によって補正された重力加速度)

def surface_acceleration(m, r, rp, phi, degrees=True):

    if degrees == True:
        phi = np.deg2rad(phi)

    # 自転角速度
    w = 2 * np.pi / rp
    
    # 万有引力による重力加速度
    g = const.G * m / r**2
    
    # 遠心力による加速度
    f = r * w**2 * np.cos(phi)
    
    # 遠心力効果によって補正された重力加速度
    gf = np.sqrt(g**2 + f**2 - 2 * g * f * np.cos(phi))
                 
    return np.array([g, f, gf])

surface_acceleration() は

 ・遠心力効果を含めない重力加速度
 ・遠心力による加速度
 ・遠心力効果によって補正された重力加速度
 
をまとめて返すので、遠心力がどの程度の影響を及ぼしているのかを知ることができます。地球の赤道上での重力加速度を計算してみましょう。rp には地球の恒星に対する自転周期 $T=86164$ 秒を渡します。この値は太陽に対する自転周期 (すなわち $1$ 日の長さ) $24$ 時間よりも短くなっています (ここでは両者の違いについて深入りはしませんが、興味のある人は天文学の本や地学の教科書などを調べてみてください)。自転角速度は $\omega=2\pi/T$ で計算できます。

# In[7]

# 地球の対恒星自転周期
rp = 86164

# 赤道上の重力加速度
g_eq = surface_acceleration(me, re, rp, 0)

print(g_eq)
# [9.79811 0.03392 9.76419]

次に北緯 $60^{\circ}$ の地点における重力加速度を計算してみます。

# In[8]

# 北緯60°における重力加速度
g_60 = surface_acceleration(me, re, rp, 60)

print(g_60)
# [9.79811 0.01696 9.78964]

赤道よりも遠心力の効果が小さくなっていることがわかりますね。実際には、遠心力以外にも、地球が完全な球ではなく楕円体であること、表面地形や密度の違いなどによっても重力加速度は影響を受けます。このように、重力加速度は場所によって異なる値が測定されるので、標準重力加速度 を $g_0=9.80665$ [m/s2] に定めて、物体の標準重量 $mg_0$ を定義しています。物理学の理論計算においても、一般に $g=9.8$ m/s2 もしくは $g=9.81$ m/s2 が使用されます。

コメント

  1. HNaito より:

    In[1] プログラムを実行すると、SciPy(1.10.1) の const.G は 6.6743e-11 になっていました。下記は誤植と思われますので、ご確認ください。
    In[3] ブログラムで、norm(vec) → np.linalg.norm(vec)
    (7)、(8) 式で、d^2x/dt → /d^2x/dt^2

    • あとりえこばと より:

      ありがとうございます。
      最新の推奨値に更新されていたんですね。
      記事を SciPy 1.10.1 に合わせて修正しておきました。