制限3体
円制限3体をルンゲクッタ法で解く
近接連星系でラグランジュ点から放たれた点の軌跡を積分する
jupyter-notebook にて program したもの。
import numpy as np
from matplotlib import pyplot as plt
from matplotlib import patches
import japanize_matplotlib
import warnings
import sys
warnings.simplefilter('ignore')
# 3体問題-1
# 2025.03.25
#
# global 変数
crk0=0.292893218813452
crk1=2-crk0
r1=np.zeros(4) # ルンゲクッタ計算用配列
r2=np.zeros(4)
r3=np.zeros(4)
r4=np.zeros(4)
q1=np.zeros(4)
q2=np.zeros(4)
q3=np.zeros(4)
q4=np.zeros(4)
yw=np.zeros(4)
y0=np.zeros(4)
y1=np.zeros(4)
y2=np.zeros(4)
y3=np.zeros(4)
y4=np.zeros(8)
a=np.zeros(8)
prm=np.zeros(16)
x1=0
x2=0
mu=0
#-----------------------------------------------------
def newton(x,a): # 5次方程式をニュートン法で解く
x3 = x # 初期値
for i in range(20):
yy1 = a[5]*x3**5 + a[4]*x3**4 + a[3]*x3**3 + a[2]*x3**2 + a[1]*x3 + a[0]
yy2 = 5.*a[5]*x3**4 + 4.*a[4]*x3**3 + 3.*a[3]*x3**2 + 2.*a[2]*x3 + a[1]
if abs(yy2) < 1.e-4:
return 0
dx = yy1/yy2
if abs(dx) < 1.e-4: # 差分が一定値以下なら
return x3
x3 = x3 -dx
return x3
#----------------------------------------------------------
def ufunc(k,yy) -> np.float64: # 計算処理
x=yy[0]
y=yy[1]
dx=yy[2]
dy=yy[3]
if k == 0:
return dx
elif k == 1:
return dy
elif k == 2:
rr1=np.sqrt((x-x1)*(x-x1)+y*y)
rr2=np.sqrt((x-x2)*(x-x2)+y*y)
fx=-mu*(x-x1)/rr1**3-(1-mu)*(x-x2)/rr2**3
ddx=dy*2.+x+fx
return ddx
elif k == 3:
rr1=np.sqrt((x-x1)*(x-x1)+y*y)
rr2=np.sqrt((x-x2)*(x-x2)+y*y)
fy=-mu*y/rr1**3-(1.-mu)*y/rr2**3
ddy=-dx*2+y+fy
return ddy
#----------------------------------------------------------
def drkgn(y0,y4,dt,nn): # ルンゲクッタ・ジル法
for k in range(nn):
q4[k]=0
yw[k]=y4[k]=y0[k]
for i in range(2):
for k in range(nn):
# y1, q1
r1[k]=ufunc(k,yw)*dt
ww=r1[k]*0.5-q4[k]
y1[k]=yw[k]+ww
q1[k]=q4[k]+ww*3-r1[k]*0.5
for k in range(nn):
# y2, q2
r2[k]=ufunc(k,y1)*dt
ww=(r2[k]-q1[k])*crk0
y2[k]=y1[k]+ww
q2[k]=q1[k]+ww*3.-r2[k]*crk0
for k in range(nn):
# y3, q3
r3[k]=ufunc(k,y2)*dt
ww=(r3[k]-q2[k])*crk1
y3[k]=y2[k]+ww
q3[k]=q2[k]+ww*3.-r3[k]*crk1
for k in range(nn):
# y4, q4
r4[k]=ufunc(k,y3)*dt
ww=(r4[k]-q3[k]*2.)/6.
yw[k]=y3[k]+ww
y4[k+i*nn]=yw[k]
q4[k]=q3[k]+ww*3.-r4[k]*0.5
return
#----------------------------------------------------------
def parm(yyy,xxx): # パラメータ計算
x=yyy[0]
y=yyy[1]
dx=yyy[2]
dy=yyy[3]
rr1=np.sqrt((x-x1)*(x-x1)+y*y)
rr2=np.sqrt((x-x2)*(x-x2)+y*y)
u=(x*x+y*y)*.5+mu/rr1+(1.-mu)/rr2
c=u*2.-(dx*dx+dy*dy)
h=(x-x1)*dy-y*dx+rr1*rr1
hkp=np.sqrt(mu*rr1)
xxx[0]=rr1
xxx[1]=rr2
xxx[2]=u
xxx[3]=c
xxx[4]=h
xxx[5]=hkp
return
#----------------------------------------------------------------
aa=np.zeros(8) # ラグランジュ点計算用係数配列
dt=0.0025 # 時間刻み幅
kk=0
# --- 質量比入力
q=float(input('* 質量比( M2/M1) : '))
if q < 1.e-2 or q > 0.95:
sys.exit() # 入力エラーで終了
# --- 質量比入力
n=int(input('* 繰り返し回数 : '))
ix=np.zeros(n)
iy=np.zeros(n)
# ---- 重心=原点とする
mu=1./(1.+q) # M1の質量割合 = 原点からM2までの距離
x1=mu-1.
x2=mu # M2の質量割合 = 原点からM1までの距離
# ---- ラグランジュ点の計算 (M2からの距離/軌道長)
aa[0]=-q
aa[1]=2*q
aa[2]=-q
aa[3]=3+q
aa[4]=-3-2*q
aa[5]=1+q
x3=0.5 # L1 の初期値
xl0 = newton(x3,aa)
xl1 = x2-xl0 # L1のX座標
print('\n ---> 質量比:','{:.3f}'.format(q),' (M2/M1) ---- 繰返し回数 :',n)
print(' ---> M1 の座標:','{:.3f}'.format(x1),' ---- M2 の座標:','{:.3f}'.format(x2))
print(' ---> L1 の座標:','{:.3f}'.format(xl1),' particle の出発点 計算刻み幅:','{:.3e}'.format(dt))
#
# ラグランジュ点の別な計算
z1=x1+0.0001
z2=x2-0.0001
zy=z1-mu/((z1-x1)*(z1-x1))+(1.-mu)/((z1-x2)*(z1-x2))
while abs(z1-z2)>1.e-5:
zm=(z1+z2)/2.
zq=zy*(zm-mu/((zm-x1)*(zm-x1))+(1.-mu)/((zm-x2)*(zm-x2)))
if zq >= 0.:
z1=zm
else:
z2=zm
xl2=zm
print(' ---> L1の別解:','{:.3f}'.format(xl2))
# ---- 初期値設定 & パラメータ計算
y0[0]=xl1 # x
y0[1]=0 # y
y0[2]=-0.2 # dx=-0.1
y0[3]=0 # dy = 0
t0=0. # 時間初期値
# ----- 初期値
parm(y0,prm)
rr1=prm[0]
rr2=prm[1]
u=prm[2]
c0=prm[3]
h=prm[4]
hkp=prm[5]
print('* init. x:','{:.3f}'.format(y0[0]),' y:','{:.3f}'.format(y0[1]),' r1:','{:.3f}'.format(rr1),' r2:','{:.3f}'.format(rr2),' h:','{:.3f}'.format(h),' u:','{:.3e}'.format(u),' c0:','{:.3e}'.format(c0),' hkp:','{:.3e}'.format(hkp))
for j in range(n):
drkgn(y0,y4,dt,4) # ルンゲクッタ・ジル法
for k in range(4):
y0[k]=y4[k+4]
parm(y0,prm)
ix[j]=x=y0[0]
iy[j]=y=y0[1]
dx=y0[2]
dy=y0[3]
rr1=prm[0]
rr2=prm[1]
u=prm[2]
c=prm[3]
h=prm[4] # 座標刻み幅
hkp=prm[5]
# if j%10==0:
# print(' loop :',j,' x:','{:.2f}'.format(x),' y:','{:.2f}'.format(y),' dx :','{:.2f}'.format(dx),' dy:','{:.2f}'.format(dy),' r1:','{:.2f}'.format(rr1),' r2:','{:.2f}'.format(rr2),' c :','{:.2e}'.format(c),' h:','{:.2e}'.format(h),' hkp:','{:.2e}'.format(hkp))
if h>=hkp:
kk=kk+1
# if abs(c-c0)> 1.e-4:
print('\n* loop : ',j+1)
title="q="+str(q)
fig=plt.figure(figsize=(12,12))
fig, ax = plt.subplots()
ax.text(0.2,0.45,title,fontsize="xx-large")
ax.set_xlim(-1.0,1.0)
ax.set_ylim(-0.8,0.8)
c0 = patches.Circle( (xl1,0), 0.005, facecolor="green", edgecolor="red")
ax.add_patch(c0)
ax.text(xl1-0.05,0.05,"L1")
c1 = patches.Circle( (x1,0), 0.01, facecolor="pink", edgecolor="red")
ax.add_patch(c1)
c2 = patches.Circle( (x2,0), 0.05, facecolor="yellow", edgecolor="black")
ax.add_patch(c2)
plt.plot(ix,iy,color="blue",lw=0.7)
#for i in range(n-1):
# plt.plot([ix[i],ix[i+1]], [iy[i],iy[i+1]], color='blue',lw=1)
plt.show()
fig.savefig('three-body.png')
---> 質量比: 0.500 (M2/M1) ---- 繰返し回数 : 1000 ---> M1 の座標: -0.333 ---- M2 の座標: 0.667 ---> L1 の座標: 0.237 particle の出発点 計算刻み幅: 2.500e-03 ---> L1の別解: 0.237 * init. x: 0.237 y: 0.000 r1: 0.571 r2: 0.429 h: 0.326 u: 1.973e+00 c0: 3.906e+00 hkp: 6.168e-01 * loop : 1000
<Figure size 1200x1200 with 0 Axes>
0.
制限3体
\(\def\bm{\boldsymbol}\)\(\def\mb{\mathbf}\)\(\def\di{\displaystyle}\)\(\def\ve{\varepsilon_0}\)\(\newcommand{\pdr}[2]{\dfrac{\partial #1}{\partial #2}}\)\(\newcommand{\ppdr}[2]{\dfrac{\partial^2 #1}{\partial #2}}\)
1.制限3体問題
制限3体問題とは、2つの質点が円運動している重力場の中で、
第3の微小質点の運動を求めるものである。
2つの質点\(\,\mathrm{P}_1,\mathrm{P}_2\,\)の質量を\(\,M_1,M_2\,\)とし質点間の距離を\(\,a\,\)とすると、
ケプラーの第3法則から
\(\quad G(M_1+M_2)=\omega^2a^3\qquad\cdots\,\)(1)
が成り立つ。
1.1.回転座標系
\(\hspace{40mm}\)
\(\hspace{80mm}\)図1.回転座標系
この系の重心を座標原点にとり、回転していない座標系\(\,(O,XYZ)\,\)『慣性座標系』での第3の質点\(\,\mathrm{P}\,\)の座標を\(\,(X,Y,Z)\,\)として、回転座標系の変換を行う。
\(\,(X,Y)\,\)平面内で角速度\(\,\omega\,\)で回転している座標系を\(\,(O,xyz)\,\)とし、
その座標系での\(\,\mathrm{P}\,\)の座標を\(\,(x,y,z)\,\)とすると
\(\left. \begin{array}{l}
X=x\cos \omega t-y\sin \omega t \\[4pt]
Y=x\sin \omega t+y\cos \omega t
\end{array} \right\}\qquad\cdots\,\)(2)
である。ここで、(2)式を時間で微分する。
\(\quad\dot{X}=\dot{x}\cos\omega t-\omega x\sin\omega t-\dot{y}\sin\omega t-\omega y\cos\omega t\)
\(\hspace{9mm}=(\dot{x}-\omega y)\cos\omega t-(\dot{y}+\omega x)\sin\omega t\qquad\cdots\,\)(3)
\(\quad\ddot{X}=(\ddot{x}-\omega\dot{y})\cos\omega t-\omega(\dot{x}-\omega y)\sin\omega t\)
\(\hspace{25mm}-(\ddot{y}+\omega x)\sin\omega t-\omega(\dot{y}+\omega x)\cos\omega t\)
\(\hspace{9mm}=(\ddot{x}-2\omega\dot{y}-\omega^2x)\cos\omega t-(\ddot{y}+2\omega\dot{x}-\omega^2y)\sin\omega t\qquad\cdots\,\)(4)
\(\quad\dot{Y}=\dot{x}\sin\omega t+\omega x\cos\omega t+\dot{y}\cos\omega t-\omega y\sin\omega t\)
\(\hspace{9mm}=(\dot{x}-\omega y)\sin\omega t+(\dot{y}+\omega x)\cos\omega t\qquad\cdots\,\)(5)
\(\quad\ddot{Y}=(\ddot{x}-\omega\dot{y})\sin\omega t+\omega(\dot{x}-\omega y)\cos\omega t\)
\(\hspace{25mm}+(\ddot{y}-\omega\dot{x})\cos\omega t-\omega(\dot{y}+\omega x)\sin\omega t\)
\(\hspace{9mm}=(\ddot{x}-2\omega\dot{y}-\omega x)\sin\omega t+(\ddot{y}+2\omega\dot{x}-\omega^2y)\cos\omega t\qquad\cdots\,\)(6)
ここで、行列表示をすると(2)式は
\(\quad\left(\begin{matrix}X\\Y\end{matrix}\right)=\left(\begin{matrix}\cos\omega t&-\sin\omega t\\\sin\omega t&\cos\omega t\end{matrix}\right)\left(\begin{matrix}x\\y\end{matrix}\right)\qquad\cdots\,\)(7)
微分の式((3)式)と二次微分の式((4)式は、それぞれ次のように表す事が出来る
\(\quad\left(\begin{matrix}\dot{X}\\\dot{Y}\end{matrix}\right)=\left(\begin{matrix}\cos\omega t&-\sin\omega t\\\sin\omega t&\cos\omega t\end{matrix}\right)\left(\begin{matrix}\dot{x}-\omega y\\\dot{y}+\omega x\end{matrix}\right)\qquad\cdots\,\)(8)
\(\quad\left(\begin{matrix}\ddot{X}\\\ddot{Y}\end{matrix}\right)=\left(\begin{matrix}\cos\omega t&-\sin\omega t\\\sin\omega t&\cos\omega t\end{matrix}\right)\left(\begin{matrix}\ddot{x}-2\omega\dot{y}-\omega^2x\\\ddot{y}+2\omega\dot{x}-\omega^2y\end{matrix}\right)\qquad\dots\,\)(9)
1.2.計算モデル
\(\hspace{40mm}\)
慣性系の\(\,X\,\)軸上に、主星(質量\(=M_1\))と伴星(質量\(=M_2\))を置き、重心に原点を置く。
ここで、\(\dfrac{M_1}{M_1+M_2}=\mu_1\,,\quad \dfrac{M_2}{M_1+M_2}=\mu_2\,\,\)とする。
主星-伴星の距離を\(\,a\,\)とすると、\(M_1\,\)の座標は\(\,(-\mu_2a,0)\,\)、\(\,M_2\,\)の座標は\(\,(\mu_1a,0)\,\)となる。
点\(\,P\,\)の座標を\(\,(x,y)\,\)とすると\(\,\bm{r_1}\,,\,\bm{r_2}\,\)の長さは
\(\quad |\bm{r_1}|=\sqrt{\left\{-\mu_2a-x\right\}^2+y^2}\qquad\cdots\,\)(10)
\(\quad |\bm{r_2}|=\sqrt{\left\{\mu_1a-x\right\}^2+y^2}\qquad\cdots\,\)(11)
なお、\(P-G_c\,\)間の\(\,\bm{r}\,\)は\(\,\,|\bm{r}|=\sqrt{x^2+y^2}\,\,\)である。
1.3.運動方程式
点\(\,P\,\)の慣性系での座標を\(\,(X,Y)\,\)とし、質量を\(\,m\,\)とすると運動方程式は
\(\quad m\frac{d^2X}{dt^2}=m\ddot{X}=F_x\qquad\cdots\,\)(12)
\(\quad m\frac{d^2Y}{dt^2}=m\ddot{Y}=F_y\qquad\cdots\,\)(13)
となる。ここでベクトル量も(2)式に従い、座標変換すると
\(\quad \left. \begin{array}{c}
f_X=f_x\cos \omega t-f_y\sin nt \\[4pt]
f_Y=f_x\sin \omega t+f_y\cos \omega t\quad
\end{array} \right\} \qquad\cdots\,\)(15)
よって、回転座標系での運動方程式は
\(\quad \left. \begin{array}{l}
\ddot{x}-2\omega\dot{y}-\omega^2x=f_x\\[4pt]
\ddot{y}+2\omega\dot{x}-\omega^2y=f_y\quad\\[4pt]
\ddot{z}=f_z
\end{array} \right\} \qquad\cdots\,\)(17)
となる。ここで、点\(\,P\,\)に対する\(M_1\,,\,M_2\,\)からの万有引力は次のようになる。
\(\quad\dfrac{GM_1m\bm{r_1}}{r_1^3}+\dfrac{GM_2m\bm{r_2}}{r_2^3}\qquad\cdots\,\)(18)
1.4.無次元化
- 長さの単位として、両星間の距離\(\,\bm{a}\,\)を採用する。
- 質量の単位として、両星の質量の和に万有引力定数を掛けた\(\,\bm{G(M_1+M_2)}\,\)を採用する。
- 時間の単位として、以下の値を採用する
ケプラーの第3法則から\(\quad\dfrac{1}{\omega}=\sqrt{\dfrac{a^3}{G(M_1+M_2}}\,\,\)で次元は\(\,[T]\,\)である。
よって、\(\quad t=\sqrt{\dfrac{a^3}{G(M_1+M_2)}}\,\tau\quad\)とおけば\(\,\bm{\tau}\,\)が無次元の量になる。
(17)式と無次元化により運動方程式は以下のようになる。
\(\quad\left. \begin{array}{l}
\ddot{x}=2\omega\dot{y}+\omega^2x+f_x\\[4pt]
\ddot{y}=-2\omega\dot{x}+\omega^2y+f_y\quad\\[4pt]
\ddot{z}=f_z
\end{array} \right\}\qquad\cdots\,\)(20)
この微分は時間\(\,\bm{\tau}\,\)についてのものである。そして、\(\,f_x\,,\,f_y\,,\,f_z\,\)は、次のようになる。
\(\quad\left. \begin{array}{l}
f_x=-\dfrac{\mu_1(x-\mu_2)}{r_1^3}-\dfrac{\mu_2(x+\mu_1)}{r_2^3}\quad\\[4pt]
f_y=-\dfrac{\mu_1y}{r_1^3}-\dfrac{\mu_2 y}{r_2^3}\\[4pt]
f_z=-\dfrac{\mu_1z}{r_1^3}+\dfrac{\mu_2 z}{r_2^3}\\[4pt]
r_1^2=(x-\mu_2)^2+y^2+z^2\\[4pt]
r_2^2=(x+\mu_1)^2+y^2+z^2
\end{array} \right\}\qquad\cdots\,\)(22)
1.5.方程式の数値積分
方程式の積分で差分法によるものは、最初に差分表の作成が必要となる。それで、差分表を作らずに数値積分する方法として、ルンゲクッタ法がある。
ルンゲクッタ法はどのような常微分方程式にも使用できる。ただ、そのためには、方程式の系を1階の常微分方程式の系に直さなければならない。そのために、次のように連立させて解法する形とする。また変数は\(\,y_1\sim y_4\,\)のラベルとする。
\(\quad\left. \begin{array}{l}
y_1=x\,,\quad y_2=y\\[3pt]
y_3=\dot{x}\,,\quad y_4=\dot{y}\quad\\[3pt]
z=0\,,\quad \dot{z}=0
\end{array} \right\} \qquad\cdots\,\)(23)
新しい変数(\(\,y_1\sim y_4\,\))での方程式は次のようになる
\(\quad\left. \begin{array}{l}
\dfrac{dy_1}{dt}=y_3\\[6pt]
\dfrac{dy_2}{dt}=y_4\\[7pt]
\dfrac{dy_3}{dt}=2y_4+y_1+f_x\\[6pt]
\dfrac{dy_4}{dt}=-2y_3+y_2+f_y\\[6pt]
f_x=-\dfrac{\mu_1(y_1-\mu_2)}{r_1^3}-\dfrac{\mu_2(y_1+\mu_1)}{r_2^3}\quad\\[7pt]
f_y=-\dfrac{\mu_1y_2}{r_1^3}-\dfrac{\mu_2 y_2}{r_2^3}\\[7pt]
r_1^2=(y_1-\mu_2)^2+y_2^2\\[4pt]
r_2^2=(y_1+\mu_1)^2+y_2^2
\end{array} \right\} \)(25)
2.方程式の解法
初期値問題として数値的に解く。
\(\quad \dfrac{dy}{dx}=f(x,y)\qquad\cdots\,\)(30)
\(x=x_0\,\)の周りで Taylor 展開すると
\(\quad y(x_0+h)=y(x_0)+hy'(x_0)+\dfrac{h^2}{2!}y”(x_0)+\cdots\qquad\cdots\,\)(31)
ここで、\(y'(x_0)=f(x_0,y_0)\,\)なので
\(\quad y(x_0+h)=y(x_0)+hf(x_0,y_0)+O(h^2)\qquad\cdots\,\)(32)
となる。この解法には次に述べる Euler法などの手法がある
2.1 Euler 法
(32)式に従って計算するのが Euler法である。
Euler 法では、点\(\,(x_ i , y_ i)\,\)での微係数をそのまま\(\,x_{i+1}\,\)まで外挿して\(\,y_{i+1}\,\)を定める。
精度を上げるためさらに高次の Taylor展開を達成する。ここで次のような式で表すことを考える。
\(y(x_0+h)=y(x_0)+\mu_1hf(x_0,y_0)+\mu_2hf(x_0+\alpha h,y_0+\beta k)\qquad\cdots\,\)(33)
右辺第二項は、点\(\,(x_0,y_0)\,\)での傾き\(\,h・f(x_0,y_0)\,\)に\(\,\mu_1\,\)の重みを掛けたものであり、第三項は、\(\,x_0\,\)と\(\,x_0+h\,\)の間のある点\((x_0+\alpha h,y_0+\beta k)\,\)での傾き\(\,h・f(x_0+\alpha h,y_0+\beta k)\,\)に\(\,\mu_2\,\)の重みを掛けたものである。即ち、第二項と第三項とで、\(x_0\,\)と\(\,x_0+h\,\)の間の、中間の\(\,x\,\)における導関数を\(\,f(x,y)\,\)の関数値の適当な平均を用いて、\(y(x_0)\,\)から\(\,y(x_0+h)\,\)を求める事を考える。ここで、
\(\quad y(x_0+h)=y(x_0)+hf(x_0,y_0)+\dfrac{h^2}{2}\left(\pdr{f}{x}+f(x_0,y_0)\pdr{f}{y}\right)+\dfrac{h^3}{6}\left[\ppdr{f}{x^2}+\cdots\right.\qquad\cdots\,\)(34)
これを、\(xy\,\)2変数同時に展開して
\(\quad f(x_0+\alpha h,y_0+\beta k)=f_0+\alpha h\left(\pdr{f}{x}\right)_0+\beta k\left(\pdr{f}{y}\right)_0\qquad\cdots\,\)(35)
この(33)式と(35)式から、
\(\quad y(x_0+h)=y(x_0)+\mu_1 hf_0+\mu_2 h\left[f_0+\alpha\left(\pdr{f}{x}\right)_0+\beta\left(\pdr{f}{y}\right)_0\right]\\[-38pt]\)
\(\hspace{21mm}=y(x_0)+(\mu_1+\mu_2)hf_0+\alpha\mu_2h^2\left(\pdr{f}{x}\right)_0+\beta\mu_2hk\left(\pdr{f}{y}\right)_0\quad\cdots\,\)(36)
となる。この(36)式と(31)式を比べれば
\(\quad \mu_1+\mu_2=1\quad,\quad \alpha\mu_2=\beta\mu_2=\dfrac{1}{2}\qquad\cdots\,\)(37)
ここで、\(k=hf_0\,\)ならば、式(37)は、Taylor展開に\(\,O(h^2)\,\)の精度まで一致することが分かる。
特に、\(\alpha=\beta=1\quad,\quad\mu_1=\mu_2=\dfrac{1}{2}\,\)のときは
\(\quad y(x_0+h)+y(x_0)+\dfrac{1}{2}\left[hf(x_0,y_0)+hf(x_0+h,y_0+hf_0)\right]\qquad\cdots\,\)(38)
これは改良 Euler 公式と呼ばれる。
2.2 Runge-Kutta法
さらに高次の項まで展開すると以下のようになる。
\(\quad y(x_0+h)=y(x_0)+\dfrac{1}{6}\left\{k_1+2(k_2+k_3)+k_4\right\}\quad\cdots\,\)(40)
より以下の形になる。
\(\quad k_1=hf(x_0,y_0)\qquad\cdots\,\)(41)
\(\quad k_2=hf\left(x_0+\dfrac{h}{2},y_0+\dfrac{k_1}{2}\right)\qquad\cdots\,\)(42)
\(\quad k_3=hf\left(x_0+\dfrac{h}{2},y_0+\dfrac{k_2}{2}\right)\qquad\cdots\,\)(43)
\(\quad k_4=hf(x_0+h,y_0+k_3)\qquad\cdots\,\)(44)
(40)式は、Taylor 展開の四次の精度までを持つ。通常、この式を Runge-Kutta の公式と呼ぶ。
\(\hspace{20mm}\)
Runge-Kutta法では、まず点\(\,(x_i,y_i)\,\)での微係数から\(\,y\,\)の増分\(\,k_1\,\)を推定する。次に、点\(\,(x_i,y_i)\,\)といま得た点\(\,(x_{i+1},y_i+k1)\,\)の中間点での微係数から、これに平行な破線で示す直線によって増分\(\,k_2\,\)の推定を行う。これを更に繰り返して、第3の増分推定\(\,k_3\,\)を得る。
次いで、第3の推定点\(\,(x_i+\varDelta x,y_i+k_3)\,\)における微係数から第4の増分推定\(\,k_4\,\)を算出する。そして最後にこれら4つの増分推定を\(\,1:2:2:1\,\)の重みを付けて平均化して、これを最終的な増分推定として\(\,y_{i+1}\,\)を定めている。
2.3 Runge-Kutta-Gill法
(40)式を改良したのが Runge-Kutta-Gill法である。この方法は式(40)~(44)の係数を変化させたものである。
\(\quad y(x_0+h)=y(x_0)+\dfrac{1}{6}\left\{k_1+(2-\sqrt{2})k_2+(2+\sqrt{2})k_3+k_4\right\}\qquad\cdots\,\)(46)
\(\quad k_1=hf(x_0,y_0)\hspace{90mm}\cdots\,\)(47)
\(\quad k_2=hf\left(x_0+\dfrac{h}{2},y_0+\dfrac{k_1}{2}\right)\hspace{60mm}\cdots\,\)(48)
\(\quad k_3=hf\left[x_0+\dfrac{h}{2},y_0+(-1+\sqrt{2})+\left(1-\dfrac{1}{\sqrt{2}}\right)k_2\right]\qquad\cdots\,\)(49)
\(\quad k_4=hf(x_0+h,y_0+k_3)\hspace{50mm}\cdots\,\)(50)
2.4 Runge-Kutta-Gill法の実際の処理
手順は以下の通り
- 出発点の\(\,(x,y)\,\)を\(\,(x_0,y_0)\,\)とする\(\qquad\)初期値として\(\,q_0=0\,\)とする
- 定数として\(\,crk0=1-1/\sqrt{2}\quad,\quad crk1=1+1/\sqrt{2}\,\)とする
- \(y(x_0)=y_0\,\)をもとに\(\,k_1=h\cdot f(x_0,y_0)\,\)を計算\(\quad k_1=h\cdot f(x_0,y_0)\)
- \(k_1\,,\,y_0\,,q_0\,\)から\(\,y_1\,\)を計算する: 1回目\(\,[q_0=0]\,\)2回目\(\,[q_0=q_4]\)(前回の)\(\,y_1=y_0+(k_1-2q_0)/2\)
- \(k_1\,,\,q_0\,\)から\(\,q_1\,\)を計算する\(\quad q_1=q_0+3(k_1-2q_0)/2-k_1/2\)
- \(x_0+h/2\,\)と\(\,y1\,\)から\(\,k_2\,\)を計算する\(\quad k_2=h\cdot f(x_0+h/2,y_1)\)
- \(y_1\,,\,k_2\,,\,q_1\,\)から\(\,y_2\,\)を計算する\(\quad y_2=y_1+\,crk0\,(k_2-q_1)\)
- \(q_1\,,\,k_2\,\)から\(\,q_2\,\)を計算する\(\quad q_2=q_1+3\,crk0\,(k_2-q_1)-\,crk0\,k_2\)
- \(x_0+h/2\,\)と\(\,y2\,\)から\(\,k_3\,\)を計算する\(\quad k_3=h\cdot f(x_0+h/2,y_2)\)
- \(y_2\,,\,k_3\,,\,q_2\,\)から\(\,y_3\,\)を計算する\(\quad y_3=y_2+\,crk1\,(k_3-q_2)\)
- \(q_2\,,\,k_3\,q_3\,(\)前回の\(),\,k2\,\)から\(\,q_3\,\)を計算する\(\quad q_3=q_2+3\,crk1\,(k_3-q_3)-crk1\,k_2\)
- \(x_0+h\,\)と\(\,y3\,\)から\(\,k_4\,\)を計算する\(\quad k_4=h\cdot f(x_0+h,y_3)\)
- \(y_3\,,\,k_4\,,\,q_3\,\)から\(\,y_4\,\)を計算する\(\quad y_4=y_3+(k_4-2q_3)/6\)
- \(q_3\,,\,k_4\,\)から\(\,q_4\,\)を計算する\(\quad q_4=q_3+(k_4-2q_3)/2-k_4/2\quad\)次回の\(\,q_0\,\)になる
- (3)から(14)の処理を繰り返して、最終出力となる
3.ラグランジュ点
3体目を放出する点として、第1ラグランジュ点を採用する。
3.1 ラグランジュ点とは
ラグランジュ点とは, 質量差のある二つの天体が共通重心の周りをそれぞれ円軌道を描いて回っているとき, この二天体に比べて質量が無視できるほど小さな宇宙機をある速度を与えてこの軌道面内に置くと, 最初の二天体との相対位置を変えずに回り続けることができる位置のことであり, 五つの点が存在する。
3.2 第1ラグランジュ点 L\(_1\)
L\(_1\,\)は質量 M\(_1\,\)と M\(_2\,\)の互いに公転する2天体を結ぶ直線上で 2体の間に存在する。太陽系の場合、一般に地球よりも太陽に近い軌道を回る物体は地球よりも短い公転周期を持つ。
これは地球から重力で引かれる効果を無視した場合の話である。物体がちょうど地球と太陽の間にあると、地球の重力の効果によって、物体が太陽に引かれる力は弱められる。物体が地球に近ければ近いほどこの効果は大きい。この効果によって L\(_1\,\)では物体の公転周期が地球の公転周期と等しい。
3.3 L\(_1\,\)の計算
2天体間の距離を\(\,R_\,\)とする。天体の質量比 \(M_2/M_1\,\)を\(\,q\,\)とする。各々の天体から重心への距離をそれぞれ\(\,R_1\,,\,R_2\,\)とすると
\(\quad R_1=R\dfrac{M_2}{M_1+M_2}=R\dfrac{q}{1+q}\qquad\cdots\,\)(60)
\(\quad R_2=R\dfrac{M_1}{M_1+M_2}=R\dfrac{1}{1+q}\qquad\cdots\,\)(61)
\(\hspace{40mm}\)
\(M_1\,\)と\(\,M_2\,\)は互いの万有引力で引き合いながら共通重心\(\,G\,\)の周りを周回する。\(M_1\,\)の軌道が円軌道とすると、ケプラーの第3法則から
\(\quad M_1R_1\omega^2=M_1\dfrac{M_2}{M_1+M_2}R\omega^2=\dfrac{GM_1M_2}{R^2}\qquad\cdots\,\)(62)
L\(_1\,\)では両星からの万有引力と重心の周りを回転する遠心力が釣り合っている。L\(_1\,\)と重心\(\,G\,\)の距離を\(\,\ell_1\,\)L\(_1\,\)と\(\,M_2\,\)の距離を\(\,\ell_x\,\)とする。L\(_1\,\)点に質量\(\,m\,\)の物体があるとすると、それの運動方程式は
\(\quad\dfrac{GmM_1}{(R-\ell_1)^2}-\dfrac{GmM_2}{\ell_1^2}=m\omega^2(R_2-\ell_1)\)
\(\hspace{40mm}=mG\dfrac{M_1+M_2}{R^3}\left(R-\dfrac{M_2}{M_1+M_2}R-R\ell_x\right)\qquad\cdots\,\)(63)
ここで、式(63)の両辺を\(\,\dfrac{GmM_1}{R^2}\,\)で割り、\(M_2/M_1=q\,\)とおくと
\(\dfrac{R^2}{(R-\ell_1)^2}-\dfrac{R^2q}{\ell_1^2}=\dfrac{1+q}{R}\left(R-\dfrac{q}{1+q}R-\ell_1\right)\quad\cdots\,\)(65)
今回のプログラムでは、ここから5次方程式を立てて、Newton法で解くやり方と、両星の中間から引力と遠心力のバランスでの繰り返し計算の2通りで解いた。
4.運動軌跡
この連立方程式をルンゲクッタで解いたものが、下の図である。条件は伴星/主星の質量比を \(q=0.5\)とし、ラグランジュ点(L\(_1\,\))から物体を放した。物体が主星の周りでリサージュを描く様子が分かる。
\(\hspace{20mm}\)
座標系は主星伴星間の距離が1となっている。主星の半径が小さい(白色矮星)場合、このような形で降着円盤が形成されるが、アルゴル型などの場合、主星の半径が 0.1(比率) を超えるので、物体は回り込まずに中心からみて右上に衝突する形になる。
![]() |
![]() |
ケプラーの法則
\(\def\bm{\boldsymbol}\)\(\def\di{\displaystyle}\)\(\def\ve{\varepsilon_0}\)\(\newcommand{\pdr}[2]{\dfrac{\partial #1}{\partial #2}}\)\(\newcommand{\ppdr}[2]{\dfrac{\partial^2 #1}{\partial #2}}\)
1.準備
1.1.極座標における速度、加速度
二次元での極座標の扱いは\(\qquad x=r\cos\theta\qquad y=r\sin\theta\qquad\cdots\,\)(1)
座標の回転では、点\(\,(x,y)\,\)を原点中心に\(\,\theta\,\)だけ回転させて、\((X,Y)\,\)となったとすると
\(\quad\left(\begin{matrix}x\\y\end{matrix}\right)=\left(\begin{matrix}\cos\theta&-\sin\theta\\\sin\theta&\,\cos\theta\end{matrix}\right)\left(\begin{matrix}X\\Y\end{matrix}\right)\qquad\cdots\,\)(2)
のように表す事が出来る。
速度は、\(x\,,\,y\,\)を時間微分することで求められる。
極座標では\(\quad\dot{x}=\dot{r}\cos\theta-r\theta\sin\theta\qquad \dot{y}=\dot{r}\sin\theta+r\dot{\theta}\cos\theta\quad\)なので
\(\quad\left(\begin{matrix}\dot{x}\\ \dot{y}\end{matrix}\right)=\left(\begin{matrix}\cos\theta&-\sin\theta \\ \sin\theta&\,\cos\theta\end{matrix}\right)\left(\begin{matrix}\dot{r}\\ r\dot{\theta}\end{matrix}\right)\quad\cdots\,\)(3)
のように表す事が出来る。変換行列は回転と同じものである。
加速度は速度の(3)式をさらに時間微分する。\(\ddot{x}\,\)と\(\,\ddot{y}\,\)は
\(\quad\ddot{x}=\ddot{r}\cos\theta-\dot{r}\dot{\theta}\sin\theta-\dot{r}\dot{\theta}\sin\theta-r\ddot{\theta}\sin\theta-r\dot{\theta}^2\cos\theta\)
\(\hspace{7mm}=(\ddot{r}-r\dot{\theta}^2)\cos\theta-(r\ddot{\theta}+2\dot{r}\dot{\theta})\sin\theta\quad\cdots\,\)(4)
\(\quad\ddot{y}=\ddot{r}\sin\theta+\dot{r}\dot{\theta}\cos\theta+\dot{r}\dot{\theta}+r\ddot{\theta}\cos\theta-r\dot{\theta}\dot{\theta}\sin\theta\)
\(\hspace{7mm}=(\ddot{r}-r\dot{\theta}^2)\sin\theta+(r\ddot{\theta}+2\dot{r}\dot{\theta})\cos\theta\quad\cdots\,\)(5)
これを行列で表すと
\(\quad\left(\begin{matrix}\ddot{x}\\ \ddot{y}\end{matrix}\right)=\left(\begin{matrix}\cos\theta&-\sin\theta \\ \sin\theta&\,\cos\theta\end{matrix}\right)\left(\begin{matrix}\ddot{r}-r\dot{\theta}^2\\ r\ddot{\theta}+2\dot{r}\dot{\theta}\end{matrix}\right)\quad\cdots\,\)(6)
1.2.運動方程式
\(x\,,\,y\,\)座標での運動方程式は、\(F_x=m\ddot{x}\,\)になる。
極座標では(6)式から、\(r\,\)方向の加速度が\(\,\ddot{r}-r\dot{\theta}^2\,\,\)となることから、
\(\quad F_r=m(\ddot{r}-r\dot{\theta}^2)\qquad\cdots\,\)(7)
となる。同じく\(\,\theta\,\)方向では加速度が\(\,r\ddot{\theta}+2\dot{r}\dot{\theta}\,\,\)であるので、
\(\quad F_{\theta}=m(r\ddot{\theta}+2\dot{r}\dot{\theta})\qquad\cdots\,\)(8)
となる。
2.ケプラーの第2法則(面積速度一定の法則)
\(\hspace{40mm}\)
ケプラーの法則は(先に、第1法則および第2法則が発見されて1609年に発表され、後に第3法則が発見されて1619年に発表された)は、ティコ・ブラーエの観測記録から、太陽に対する火星の運動を推定し、3つの法則として定式化した。
第2法則は惑星と太陽とを結ぶ線分が単位時間に掃く面積(面積速度度)は、一定であるというものである。
2.1 導出
運動方程式で\(\,\theta\,\)方向は(8)式から、\(F_{\theta}=m(r\ddot{\theta}+2\dot{r}\dot{\theta})\,\)となるが、
\(\quad\dfrac{d}{dt}\left(r^2\dot{\theta}\right)=r^2\ddot{\theta}+2r\dot{r}\dot{\theta}=r(r\ddot{\theta}+2\dot{r}\dot{\theta})\quad\) となるので
\(\quad F_{\theta}=m\dfrac{1}{r}\dfrac{d}{dt}\left(r^2\dot{\theta}\right)\quad\cdots\,\)(9)
と表す事が出来る。
2.2 角運動量保存
角運動量は運動量\(\,\bm{p}\,\)と位置ベクトル\(\,\bm{x}\,\)を使って以下のように表わすことが出来る
\(\quad\bm{p}\times\bm{x}=\bm{L}\qquad\cdots\,\)(10)
ここで単位質量当たりの角運動量\(\,\bm{k}\,\)を考える。
\(\quad\bm{k}=\dot{\bm{x}}\times\bm{x}\hspace{10mm}\cdots\hspace{30mm}\)これを時間微分すると
\(\quad\dfrac{d\bm{k}}{dt}=\ddot{\bm{x}}\times\bm{x}+\dot{\bm{x}}\times\dot{\bm{x}}\quad\cdots\,\)(11)
(11)式の右辺の第二項は、同じベクトルの外積なので\(\,0\,\)となる。第一項も中心力での加速度は位置ベクトルと平行なので\(\,0\,\)になる。
よって、単位角運動量\(\,\bm{k}\,\)の時間微分は\(\,0\,\)になるので、定数ベクトルとなる。この単位角運動量に質量\(\,m\,\)を乗じた\(\,\bm{L}\,\)も定数となるので、角運動量は保存される。
2.3 面積速度
極座標における面積速度は、単位時間に\(\bm{r}\,\)と\(\,r\theta\,\)によって作られる三角形の面積となる。(9)式からも、中心力では角度成分への力はないので、\(F_{\theta}=0\,\)となる。
時間\(\,t\,\)で積分すると、\(\quad r^2\dot{\theta}=h\qquad h\,\):積分定数\(\quad\cdots\,\)(12)
が導かれる。(\(\,h\,\)は時間項を含まないので定数となる)
3.ケプラーの第1法則(楕円軌道の法則)
惑星は、太陽を焦点のひとつとする楕円軌道上を動く。太陽の位置を原点に取り、太陽と惑星の距離を\(\,r\) 、 真近点角を\(\,\theta\,\)をパラメータとする極座標では、惑星の軌道は次の式で与えられる。
\(\quad r=\dfrac{p}{1+\ve\,\cos\theta}\qquad\cdots\,\)(13)
\(\hspace{47mm}\)
ここで、 \(p\,\)は半通径(semi-latus rectum)、\(\ve\,\)は楕円の離心率である。ただし\(\,0\,\ge\,\ve\,>\,1\,\)であり、\(\ve\,=0\,\)のとき、太陽中心の円軌道を表す。
3.1 導出
3.1.1 step-1
惑星の質量を\(\,m\,\)、太陽の質量を\(\,M\,\)、万有引力定数を\(\,G\,\)とすると、惑星の運動方程式は
\(\quad F_r=m\ddot{\bm{x}}=-GMm\dfrac{\bm{x}}{|\bm{x}|^3}\qquad\cdots\,\)(15)
これを極座標で表すと、(9)式から
\(\quad m(\ddot{r}-r\dot{\theta}^2)=-\dfrac{GMm}{r^2}\qquad\cdots\,\)(16)
この式に (12)式を代入すると
\(\quad m\ddot{r}-\dfrac{mh^2}{r^3}=-\dfrac{GMm}{r^2}\hspace{15mm}\)両辺を\(\,m\,\)で割ると
\(\quad\ddot{r}-\dfrac{h^2}{r^3}=-\dfrac{GM}{r^2}\qquad\cdots\,\)(17)
3.1.2 step-2
(14)式と合成関数の微分より
\(\quad\dot{r}=\dfrac{dr}{d\theta}=\dfrac{d\theta}{dt}\dfrac{dr}{d\theta}=\dot{\theta}\dfrac{dr}{d\theta}=\dfrac{h}{r^2}\dfrac{dr}{d\theta}\qquad\cdots\,\)(18)
ここで\(\,p=\dfrac{1}{r}\,\)という変数を導入する。これを\(\,\theta\,\)で微分すると
\(\quad\dfrac{dp}{d\theta}=-\dfrac{1}{d\theta}\left(\dfrac{1}{r}\right)=\dfrac{dr}{d\theta}\dfrac{1}{r^2}\qquad\cdots\,\)(19)
この(19)式を(18)式に代入すると、以下の式が得られる
\(\quad\dot{r}=-h\dfrac{dp}{d\theta}\qquad\cdots\,\)(20)
さらに時間微分すると
\(\quad\ddot{r}=-h\dfrac{d}{dt}\left(\dfrac{dp}{d\theta}\right)=-h\dfrac{d\theta}{dt}\dfrac{d}{d\theta}\left(\dfrac{dp}{d\theta}\right)=-h\dfrac{d\theta}{dt}\dfrac{d^2p}{d\theta^2}\)
\(\hspace{9mm}=-h\,\dot{\theta}\,\dfrac{d^2p}{d\theta^2}=-h\dfrac{h}{r^2}\dfrac{d^2p}{d\theta^2}\qquad\because\quad r^2\dot{\theta}=h\)
\(\hspace{9mm}=-\dfrac{h^2}{r^2}\dfrac{d^2p}{d\theta^2}\qquad\cdots\,\)(21)
この(21)式を(17)式に代入すると
\(\quad-\dfrac{h^2}{r^2}\dfrac{d^2p}{d\theta^2}-\dfrac{h^2}{r^3}=-\dfrac{GM}{r^2}\qquad\cdots\,\)(22)
整理すると
\(\quad\dfrac{d^2p}{d\theta^2}=-p+K\hspace{15mm}\)ここで\(\quad \dfrac{GM}{h^2}=K\quad\)と置いた。
このような振動方程式の一般解は
\(\quad p=A\cos(\theta+B)+K\qquad\Longrightarrow\quad r=\dfrac{1}{A\cos(\theta+B)+K}\qquad\cdots\,\)(23)
ここで長軸を\(\,y\,\)軸に取り、\(\,B=0\,\)、\(A/K=\ve\,\,,\,\,1/K=\ell\,\)と置くと
\(\quad r=\dfrac{\ell}{1+\ve\cos\theta}\qquad\cdots\,\)(24)\(\hspace{20mm}\)となり、第一法則が導出された。
3.ケプラーの第3法則
\(\hspace{45mm}\)
近日点と遠日点での距離と速度をそれぞれ\(\,R_1\,\,,\,v_1\,\)と\(\,R_2\,\,,\,\,v_2\,\)とする。
近日点と遠日点の面積速度\(\,V_s\,\)は同じなので
\(\quad V_s=\dfrac{1}{2}v_1R_1=\dfrac{1}{2}v_2R_2\qquad\cdots\,\)(26)
またエネルギー保存則より、近日点と遠日点のエネルギーは等しいので
\(\quad E_s=\dfrac{1}{2}mv_1^2-\dfrac{GMm}{R_1}=\dfrac{1}{2}mv_2^2-\dfrac{GMm}{R_2}\qquad\cdots\,\)(27)
(26)式と(27)式から
\(\quad v_1R_1=v_2R_2\qquad\longrightarrow\qquad\dfrac{v_2}{v_1}={R_1}{R_2}\qquad\cdots\)(28)
\(\quad\dfrac{1}{2}v_1^2\left(1-\dfrac{v_2^2}{v_1^2}\right)=\dfrac{GM}{R_1}\left(1-\dfrac{R_1}{R_2}\right)\qquad\)(28)式を代入すると
\(\quad\dfrac{1}{2}v_1^2\left(1-\dfrac{R_1^2}{R_2^2}\right)=\dfrac{GM}{R_1}\left(1-\dfrac{R_1}{R_2}\right)\qquad\cdots\,\)(29)
(29)式を\(\,\left(1-\frac{R_1}{R_2}\right)\,\)で割ると
\(\quad\dfrac{1}{2}v_1^2\left(1+\dfrac{R_1}{R_2}\right)=\dfrac{GM}{R_1}\qquad\)これにより
\(\quad v_1=\sqrt{\dfrac{2R_2GM}{R_1(R_1+R_2)}}\qquad\cdots\,\)(30)
面積速度は
\(\quad V_s=\dfrac{1}{2}v_1R_1=\frac{1}{2}R_1\sqrt{\dfrac{2R_2GM}{R_1(R_1+R_2)}}=\sqrt{\dfrac{R_1R_2GM}{2(R_1+R_2)}}\qquad\cdots\,\)(31)
\(\quad\)ここで楕円の性質より \(\qquad R_1+R_2=2a\quad , \quad R_1R_2=b^2\)
\(\quad\) ここで楕円の面積より \(\qquad S=\pi ab \qquad\)これらを代入すると
\(\quad V_s=\dfrac{b}{2}\sqrt{\dfrac{GM}{a}}\qquad\cdots\,\)(33)
\(\quad S=V_sT=\dfrac{b}{2}\sqrt{\dfrac{GM}{a}}T=\pi ab\)
\(\quad \dfrac{GM}{a}=\left(\dfrac{2\pi a}{T}\right)\quad\longrightarrow\quad\dfrac{T^2}{a^3}=\dfrac{4\pi^2}{GM}\qquad\cdots\,\)(34)\(\qquad\)が導かれる。
![]() |
![]() |
音速
\(\def\bm{\boldsymbol}\)\(\def\mb{\mathbf}\)\(\def\di{\displaystyle}\)\(\def\ve{\varepsilon_0}\)\(\newcommand{\pdr}[2]{\dfrac{\partial #1}{\partial #2}}\)\(\newcommand{\ppdr}[2]{\dfrac{\partial^2 #1}{\partial #2}}\)
音速を流体力学的に導出する。中央大学 棚橋先生の「流体力学」を参考させてもらいました。
1.音波
Wikipediaでは、「気体、液体、固体を問わず、弾性体を伝播するあらゆる弾性波の総称をさす。」となっている。流体として考えると、密度の変化が音速で伝わっていくものということになる。密度が上がると圧力が上がるが、そのために生じる圧力勾配力は流体の圧力を下げる向きにかかるため、密度の値をもとの値に戻そうとする復元力としての役割を果たす。この復元力のために流体の振動が発生するが、これが空間的に伝わっていくのが音波の伝搬である。
1-1: モデル化
\(x\,\)軸方向に伝搬する平面波の音波について考える。全体の密度は\(\,\rho_0\,\)一定とする。ここで、密度が微小に変化した場合を考える。
\(\quad \rho(t,x)=\rho_0+\delta\rho(t,x)\)
\(\qquad\,\,|\delta\rho(t,x)|\ll\rho_0\qquad\cdots\,\)(1)
微小の密度(圧力)変化に対応して、速度場も微小に変化する。音波がない場合は流体は静止しているので、静止速度場はゼロである。
\(\qquad\bm{v}(t,x)=v_0+\delta v(t,x)=\delta v(t,x)\qquad\cdots\,\)(2)
微小圧力変化\(\,\delta\rho(t,x)\,\)とそれに対応する流体の平衡位置からのズレが\(\,\delta v(t,x)\,\)である。
状態方程式を\(\,P=P(\rho)\,\)とする。
1-2: 運動方程式
流体力学の連続の式(質量保存則)とオイラー方程式(粘性を含まない運動量保存則)を適用する。
・連続の式
\(\qquad\pdr{\rho}{t}+\nabla\cdot(\rho\bm{v})=0\qquad\cdots\,\)(3)
\(\qquad\pdr{\rho}{t}+\nabla\cdot(\rho\bm{v})=\pdr{}{t}(\rho_0+\delta\rho)+\nabla\cdot[(\rho_0+\delta\rho)\delta\bm{v}]\)
\(\qquad\quad=\pdr{\delta\rho}{t}+\rho_0\nabla\cdot\delta\bm{v}+\nabla\cdot(\delta\rho\delta\bm{v})\simeq\pdr{\delta\rho}{t}+\rho_0\pdr{\delta v}{x}=0\qquad\cdots\,\)(4)
\(\qquad\quad\delta\rho\delta\bm{v}\,\)項は2次の微少量のため無視した。\(\quad(\quad\because\quad\delta\rho\ll\rho_0\quad)\)
・オイラー方程式
\(\qquad\pdr{\bm{v}}{t}+(\bm{v}\cdot\nabla)\bm{v}=-\dfrac{1}{\rho}\nabla P\qquad\cdots\,\)(5)
まず 左辺は\(\quad\pdr{\delta\bm{v}}{t}+(\delta\bm{v}\cdot\nabla)\delta\bm{v}=\pdr{\delta v}{t}+\pdr{\delta\bm{v}}{x}\delta\bm{v}\sim\pdr{\delta v}{t}\qquad\cdots\,\)(6)
ここで\(\,\partial\bm{v}/\partial t \delta\bm{v}\,\)は2次の微少量なので無視した。
右辺は\(\quad-\dfrac{1}{\rho}\nabla P(\rho)=\dfrac{1}{\rho}\dfrac{dP}{d\rho}\nabla\rho=-\dfrac{1}{\rho_0+\delta\rho}\left.\dfrac{dP}{d\rho}\right|_{\rho_0+\delta\rho}\nabla(\rho_0+\delta\rho)\qquad\cdots\,\)(7)
ここで各項目の計算をまとめる。
\(\qquad\bullet\quad\dfrac{1}{\rho_0+\delta\rho}=\dfrac{1}{\rho_0}\left\{1-\dfrac{\delta\rho}{\rho_0}+\dfrac{1}{2}\left(\dfrac{\delta\rho}{\rho_0}\right)^2-\cdots\right\}\quad\simeq\quad\dfrac{1}{\rho_0}\)
\(\qquad\bullet\quad\left.\dfrac{dP}{d\rho}\right|_{\rho_0+\delta\rho}\sim\left.\dfrac{dP}{d\rho}\right|_{\rho_0}+\left.\dfrac{d^2P}{d\rho^2}\right|_{\rho_0}\delta\rho\quad\simeq\quad\left.\dfrac{dP}{d\rho}\right|_{\rho_0}\)
\(\qquad\bullet\quad\nabla(\rho_0+\delta\rho)\quad\simeq\quad\nabla\delta\rho\qquad(\quad\because\quad\rho_0\sim const.\quad)\)
\(\quad\pdr{\delta v}{t}\,\simeq\,-\dfrac{1}{\rho_0}\dfrac{dP}{d\rho}\nabla\delta\rho=-\dfrac{1}{\rho_0}\dfrac{dP}{d\rho}\pdr{\delta\rho}{x}\qquad\cdots\,\)(8)
\(\blacklozenge\,\,\,\)ここで、(4)式と(8)式を再掲する。
\(\quad\left\{\begin{array}{lr}\pdr{\delta\rho}{t}+\rho_0\pdr{\delta v}{x}=0&\cdots\,\mathrm{(4)}\\
\pdr{\delta v}{t}=-\dfrac{1}{\rho_0}\dfrac{dP}{d\rho}\pdr{\delta\rho}{x}&\cdots\,\mathrm{(8)}\end{array}\right.\)
\(\quad\)ここで、(4)式を時間微分する。
\(\qquad\ppdr{\delta\rho}{t^2}+\rho_0\pdr{}{t}\pdr{\delta v}{x}=\ppdr{\delta\rho}{t^2}+\rho_0\pdr{}{x}\pdr{\delta v}{t}=0\qquad\cdots\,\)(9)
\(\quad\)(9)式に(8)式を代入する。
\(\qquad\ppdr{\delta\rho}{t^2}+\rho_0\pdr{}{x}\left(-\dfrac{1}{\rho_0}\dfrac{dP}{d\rho}\pdr{\delta\rho}{x}\right)=\ppdr{\delta\rho}{t^2}-\dfrac{dP}{d\rho}\ppdr{\delta\rho}{x^2}=0\qquad\cdots\,\)(10)
ここで\(\quad\dfrac{dP}{d\rho}=c^2\quad\)とおくと、式(10)は次のようになる。
\(\qquad\ppdr{\delta\rho}{t^2}=c^2\ppdr{\delta\rho}{x^2}\qquad\cdots\,\)(11)
速度\(\,c\,\)の進行波となる。\(\qquad c=\sqrt{\dfrac{dP}{d\rho}}\)
1-3: 圧縮率で音速を表す
圧縮率\(\,\kappa\,\)とは系にかかる圧力に対して、系の体積がどの程度変化するかを表す状態量であり、以下の式で表される。
\(\qquad\kappa=-\dfrac{1}{V}\dfrac{dV}{dP}\qquad\cdots\,\)(12)\(\hspace{20mm}\)【参考】:個体弾性率\(\,\,\, K=-V\dfrac{dP}{dV}\qquad\)
音速の二乗は
\(\quad c^2=\dfrac{dP}{d\rho}=\dfrac{dV}{d\rho}\dfrac{dP}{dV}\quad\cdots\,\)(13)
ここで\(\quad V=\dfrac{M}{\rho}\quad\)なので、\(\quad \dfrac{dV}{d\rho}=\dfrac{d}{d\rho}\dfrac{M}{\rho}=-\dfrac{M}{\rho^2}=-\dfrac{V^2}{M}\qquad\)これを式(13)に代入すると
\(\qquad\dfrac{dP}{d\rho}=-\dfrac{V^2}{M}\dfrac{dP}{dV}=-\dfrac{V}{\rho}\dfrac{dP}{dV}=\dfrac{1}{\kappa\rho}\qquad\cdots\,\)(14)
よって音速は
\(\qquad c=\sqrt{\dfrac{1}{\kappa\rho}}\qquad\cdots\,\)(15)
1-4: 気体の音速
理想気体で考える。理想気体の状態方程式は、\(\mu\,\)を\(1\)モルの質量として、
\(\qquad P=RT\dfrac{n}{V}=RT\dfrac{M/V}{M/n}=RT\dfrac{\rho}{\mu}\qquad\cdots\,\)(20)
圧力\(\,P\,\)の全微分は
\(\qquad dP=\left(\pdr{P}{T}\right)_{\rho}dT+\left(\pdr{P}{\rho}\right)_Td\rho\qquad\cdots\,\)(21)
\(\qquad\left(\pdr{P}{T}\right)_{\rho}=\dfrac{R\rho}{\mu}=\dfrac{P}{T}\quad,\qquad\left(\pdr{P}{\rho}\right)_T=\dfrac{RT}{\mu}=\dfrac{P}{\mu}\)
これらを式(12)に代入すると
\(\qquad dP=\dfrac{P}{T}dT+\dfrac{P}{\rho}d\rho\qquad\Longrightarrow\quad \dfrac{dP}{P}-\dfrac{d\rho}{\rho}=\dfrac{dT}{T}\qquad\cdots\,\)(22)
\(\blacklozenge\,\,\,\)気体の内部エネルギー\(\,U\,\)は定積比熱より求められる。
\(\qquad dU=C_vdT\qquad\)・単原子分子:\(\,C_v=\dfrac{3}{2}R\qquad\)二原子分子:\(\,C_v=\dfrac{5}{2}R\qquad\cdots\,\)(23)
熱力学の第一法則に」より
\(\qquad dU=-PdV+dQ’\qquad\cdots\,\)(24)
音の振動は断熱的なので\(\,dQ’\,\)をゼロとすると
\(\qquad dU=-PdV=-Pd\dfrac{\mu}{\rho}=P\mu\dfrac{d\rho}{\rho^2}=RT\dfrac{\rho}{\mu}\mu\dfrac{d\rho}{\rho^2}=RT\dfrac{d\rho}{\rho}\qquad\cdots\,\)(25)
\(\blacklozenge\,\,\,\)2原子分子を考え、式(23)と式(25)の係数を比較すると、以下の関係が導かれる。
\(\qquad\dfrac{5}{2}\dfrac{dT}{T}=\dfrac{d\rho}{\rho}\qquad\cdots\,\)(26)
式(26)と式(22)から
\(\qquad\dfrac{dT}{T}=\dfrac{2}{5}\dfrac{d\rho}{\rho}=\dfrac{dP}{P}-\dfrac{d\rho}{\rho}\quad\Longrightarrow\quad\dfrac{dP}{P}=\dfrac{7}{5}\dfrac{d\rho}{\rho}\qquad\cdots\,\)(27)
音速は式(27)と式(20)より
\(\qquad c=\sqrt{\dfrac{dP}{d\rho}}=\sqrt{\dfrac{7}{5}\dfrac{P}{\rho}}=\sqrt{\dfrac{7}{5}\dfrac{RT}{\mu}}\qquad\cdots\,\)(28)
\(\blacklozenge\,\,\,\)空気( 平均分子量:29 温度:0℃ )で計算してみると
\(\qquad\sqrt{\dfrac{7}{5}\dfrac{RT}{\mu}}=\sqrt{\dfrac{7}{5}\dfrac{8.314\times\,273}{29.0\times 10^{-3}}}\simeq 331.0\,\,\)m/s
おおよそ再現出来ている。
◆ 参考資料
- 流体力学 講義ノート 第9回 流体の波動
:https://www2.yukawa.kyoto-u.ac.jp/~norihiro.tanahashi/pdf/hydrodynamics/note_hydrodynamics_9.pdf
ガス粒子の衝突
1.ガス粒子の衝突
宇宙での流体を考える際、ガスの力学的な時間に対して電子・分子の衝突時間が十分に短い場合、流体近似が成り立つ。 これは速度空間が十分に緩和して、連続体として見なせる状態である。気体を1種類の分子からなっており、半径\(\,r\,\)の剛体球で、ある分子1個が一定の速さ\(\,v\,\)で運動し、他の分子は静止しているという前提で衝突を考える。
分子間の距離が\(\,2r\,\)になると衝突するので、衝突断面積は\(\,2\pi(2r)^2\,\)であり、単位体積当たり\(\,n\,\)個の分子が分布する中を速さ\(\,\langle v\rangle\,\)で走る分子は平均時間間隔\(\quad\tau=\dfrac{1}{n\pi(2r)^2\langle v\rangle}\quad\)で衝突を繰り返すことになる。
この時間\(\,\tau\,\)のことを平均自由時間または緩和時間という。また、分子は衝突の間に平均\(\quad l_M=\langle v\rangle r=\dfrac{1}{n\pi(2r)^2}\quad\)の距離だけ動くことが出来る。この\(\,l_M\,\)を平均自由行程と呼ぶ。
\(\star\,\)まず、全部の分子が同じ速度で運動している場合を考える。
1つの分子の速度(速さではなく)を\(\,\bm{v}\,\)とし、もう一つの分子の速度を\(\,\bm{v}’\,\)とする。2つの分子の相対速度を\(\,\bm{v}_r\,\)とすると、これらの関係は左図のようになる。
余弦定理から\(\quad\bm{v}_r^2=\bm{v}^2+\bm{v}’^2-2\,\bm{vv}’\cos\theta\quad\)の関係を導く。ここで\(\,|\bm{v}|=|\bm{v}’|\,\)とすると、\(\,v_r\,\)は
\(\,v_r=\sqrt{2v^2-2v^2\cos\theta}=v\sqrt{4\dfrac{1-\cos\theta}{2}}=2v\sin\frac{\theta}{2}\) となる。
同じ\(\,\theta\,\)をもつ衝突について、紙面の内外での確立に差がないとして、立体角の要素を計算する。
\(\quad\di\int_0^{2\pi}\int_0^{\pi}\sin\theta\,d\theta\,d\omega=\int_0^{\pi}2\pi\sin\theta\,d\theta\quad\)そして、\(\,v_r\,\)の平均値は次のようになる。
\(\,\bar{v}_r=\dfrac{\di\int_0^{\pi}2v\sin\frac{\theta}{2}2\pi\sin\theta d\theta}{\di\int_0^{\pi}2\pi\sin\theta d\theta}=\dfrac{4\pi v\di\int_0^{\pi}\sin\frac{\theta}{2}\frac{1}{2}\sin\frac{\theta}{2}\cos\frac{\theta}{2}d\theta}{2\pi\left[-\cos\theta\right]_0^{\pi}}=2v\di\int_0^{\pi}\sin^2\frac{\theta}{2}\cos\frac{\theta}{2}d\theta\)
ここで\(\,\phi=\theta/2\,\)と置くと、\(d\theta=2d\phi\,\)となり、\(v_r\,\)は
\(\,v_r=4v\di\int_0^{\pi/2}\sin^2\phi\cos\phi\,d\phi=2v\di\int_0^{\pi/2}\sin 2\phi\cos\phi\,d\phi=2v\di\int_0^{\pi/2}\dfrac{1}{2}\left(\sin 3\phi\,\sin\phi\right)d\phi\)
よって、\(v_r=v\left[-\dfrac{1}{3}\cos 3\phi-\cos\phi\right]_0^{\pi/2}=\dfrac{4}{3}v\quad\)となる。
\(\star\,\)次に、速度が違う場合を考える。
同じく余弦定理\(\quad\bm{v}_r^2=\bm{v}^2+\bm{v}’^2-2\,\bm{vv}’\cos\theta\quad\)から同じように立体角と角度の積分の計算するが、ちょっと複雑なので結果だけを示す。\(\,v>v’\,\)の場合、相対速度は
\(\quad v_r=\dfrac{3v^2+v’^2}{3v}\quad\)となる。
\(\star\,\)続いて、全ての分子が Maxwell-Boltzmann 分布で運動している場合を考える。
速度空間において、速さが\(\,v\,\)と\(\,v+dv\,\)の間にある分子数を求める。そのためには\(F(v_x,v_y,v_z)dv_xdv_ydv_z\,\)を\(\,v\,\)と\(\,v+dv\,\)の間にある部分について\(\,v_x,v_y,v_z\,\)の積分をする。
この部分の球殻の体積は\(\,4\pi v^2dv\,\)であるので、この球殻に含まれる分子数は、全体の分子数を\(\,N\,\)、分子の質量を\(\,m\,,\,\alpha=\dfrac{m}{2k_BT}\,\)とおくと
\(F(v_x,v_y,v_z)dv=N\left(\dfrac{\alpha}{\pi}\right)^{3/2}\exp(-\alpha(v_x^2+v_y^2+v_z^2)4\pi v^2dv\,\)となり、
\(F(v)dv=4\pi N\left(\dfrac{\alpha}{\pi}\right)^{3/2}v^2\,e^{-\alpha\,v^2}dv\quad\)となる。
以降、\(\,F(v)\,\)は割合を表すととし全体の分子数を省略する。
\(\quad F(v)=\alpha^{3/2}\dfrac{4}{\sqrt{\pi}}\exp(-\alpha v^2)dv\qquad\)とする。
これは水素分子の各温度毎(K)の速度分布を計算したものである。
\(\,\bullet\,\)分布関数の最大値\(\,v_{max}\,\)を求める。
\(\quad\dfrac{dF}{dv}=4\pi\left(\dfrac{\alpha}{\pi}\right)^{3/2}(2v-2\alpha v^3)e^{-\alpha v^2}=0\quad\)より
\(\quad 1-\alpha v^2=0\quad\)なので、\(\,v_{max}=\dfrac{1}{\sqrt{\alpha}}\)
\(\,\bullet\,\)続いて分布関数の平均値\(\,\langle v\rangle\,\)を求める。
\(\quad\langle v\rangle=4\pi\left(\dfrac{\alpha}{\pi}\right)^{3/2}\di\int_0^{\infty}v^3e^{-\alpha v^2}dv\quad\)ここで、\(v^2=w\quad\)とおくと、\(\,dv=\dfrac{1}{2v}dw\quad\)より
\(\quad\di\int_0^{\infty}v^3e^{-\alpha v^2}dv=\di\int_0^{\infty}\dfrac{w}{2}e^{-\alpha w}dw=\left[\dfrac{w}{2}e^{-\alpha w}\dfrac{(-1)}{\alpha}\right]_0^{\infty}-\di\int_0^{\infty}\dfrac{1}{2}e^{-\alpha w}\dfrac{(-1)}{\alpha}dw\)
\(\qquad =\dfrac{(-1)}{2\alpha}\left[\dfrac{w}{e^{\alpha w}}\right]_0^{\infty}+\dfrac{1}{2\alpha}\di\int_0^{\infty}e^{-\alpha w}dw=0+\dfrac{1}{2\alpha}\dfrac{(-1)}{\alpha}\left[e^{-\alpha w}\right]_0^{\infty}=\dfrac{1}{2\alpha^2}\quad\)なので
\(\quad\langle v\rangle=4\pi\left(\dfrac{\alpha}{\pi}\right)^{3/2}\dfrac{1}{2\alpha^2}=\dfrac{2}{\sqrt{\pi\alpha}}\quad\)と計算できる。
\(\,\bullet\,\)続いて標的分子が分布関数に従う場合の相対速度\(\,v_r’\,\)を求める。
速度が違う場合の\(\,v_r=\dfrac{3v^2+v’^2}{3v}\,\)の式と分布関数\(\,F(v)=\alpha^{3/2}\dfrac{4}{\sqrt{\pi}}e^{-\alpha v^2}dv\,\)を使用する。
\(\quad v_r’=\alpha^{3/2}\dfrac{4}{\sqrt{\pi}}\left\{\di\int_0^v\dfrac{3v^2+v’^2}{3v}v’^2\,e^{-\alpha v’^2}dv’+\di\int_v^{\infty}\dfrac{3v’^2+v^2}{3v’}v’^2\,e^{-\alpha v’^2}dv’\right\}\)
この計算は最終的にGauss関数(誤差関数)を含む以下の形になる。
\(\quad v_r’=\dfrac{1}{\sqrt{\alpha\pi}}\left\{e^{-\alpha v^2}+\left(\dfrac{\alpha v}{2}+\dfrac{1}{v}\right)\di\int_0^ve^{-\alpha v’^2}dv’\right\}\)
\(\,\bullet\,\)前項の結果を用いて、注目分子の速さの分布関数での相対距離\(\,v_r\,\)を求める。
\(\quad v_r=\alpha^{3/2}\dfrac{4}{\sqrt{\pi}}\dfrac{1}{\sqrt{\alpha\pi}}\di\int_0^{\infty}v^2e^{-\alpha v^2}\left\{e^{-\alpha v^2}+\left(\dfrac{\alpha v}{2}+\dfrac{1}{v}\right)\di\int_0^ve^{-\alpha v’^2}dv’\right\}dv\)
\(\qquad =\dfrac{4\alpha}{\pi}\left\{\di\int_0^{\infty}v^2e^{-2\alpha v^2}dv+\di\int_0^{\infty}v^2e^{-\alpha v^2}\left(\dfrac{\alpha v}{2}+\dfrac{1}{v}\right)\di\int_0^ve^{-\alpha v’^2}dv’\right\}dv\)
第一項の積分はガウス積分だが、第二項の計算は相当複雑である。最終結果は以下のようになる。
\(\quad v_r=\sqrt{\dfrac{8}{\pi\alpha}}=4\sqrt{\dfrac{k_BT}{\pi m}}=\sqrt{2}\,\langle v\rangle\quad\)のように、平均速さの\(\,\sqrt{2}\,\)倍になる。
2.宇宙流体での衝突
\(\star\,\)水素分子
水素分子で考える。分子半径が\(\,\,0.145\,\)nm なので、衝突断面積(\(\,\sigma\,\))は、
\(\quad\sigma=4\pi r^2=26.45\times 10^{-16}\)cm\(^2\,\)となる。
分子質量\(\,\,m\,\,\)は\(\,\,m=\dfrac{2}{N_A}=\dfrac{2}{6.02\times 10^{23}}=3.32\times 10^{-24}\,\)g なので、速さ\(\,\langle v\rangle\,\)は、
\(\quad\langle v\rangle =\sqrt{\dfrac{8k_BT}{\pi m}}=\sqrt{\dfrac{8\times 1.38\times 10^{-16}T}{\pi\times 3.32\times 10^{-24}}}=1.03\times 10^4\sqrt{T}\,\)cm/s となる。
相対速度\(\,\,v_r\,\,\)は\(\quad v_r=\sqrt{2}\,\langle v\rangle=1.46\times 10^4\sqrt{T}\,\,\)cm/s である。
ガス密度を\(\,\rho\,\)(g/cm\(^3\))\(\,\)とすると、分子密度\(\,n\,\)は、\(\,n=\rho/(2/N_A)=\rho\times 3.01\times 10^{24}\,\)となる。
緩和時間\(\,\tau\,\)(秒)は、\(\,\tau=\dfrac{1}{n\sigma\,v_r}=\dfrac{1}{\rho\times 3.01\times 10^{24}\times 26.45\times 10^{-16}\times 1.46\times 10^4\sqrt{T}}\)
最終的に\(\quad \tau=\dfrac{1}{1.16\times 10^{14}\,\rho\sqrt{T}}=\dfrac{1}{3.86\times 10^{-11}n\sqrt{T}}\,\)(秒)
\(\qquad =\dfrac{1}{3.86\times 10^{-11}n\sqrt{T}\times 3600\times\ 24 \times 365.25}=\dfrac{1}{1.22\times 10^{-3}n\sqrt{T}}\,\,\)(年)\(\quad\)となる。
星間雲\(\quad \rho=10^{-22}\,\)g cm\(^{-3}\quad n=301\,\,\)分子 cm\(^{-3}\quad T=10\,\)Kで考えると
\(\quad\tau=\dfrac{1000}{1.22\times 301\times\sqrt{10}}=0.86\,\,\)(年) となる。
\(\star\,\)水素プラズマ
電子の衝突は、電子の運動エネルギーとクーロンポテンシャルが同程度になる距離なので
\(\quad \dfrac{1}{2}m_ev_e^2\,\,\sim\,\,\dfrac{e^2}{4\pi\epsilon_0r_0}\quad\Longrightarrow\quad r_0\,\,\sim\,\,\dfrac{e^2}{2\pi\epsilon_0m_ev_e^2}\quad\)となる。
衝突断面積\(\,\sigma\,\)は\(\quad\sigma=4\pi r_0^2\,\,\sim\,\,\dfrac{4\pi e^2}{4\pi^2\epsilon_0^4m_e^2v_e^4}\,\sim\,\dfrac{e^2}{\epsilon_0^2m_e^2v_e^4}\)
以降SI単位として計算する。\(\,k_B\,:\,1.38\times 10^{-23}\,\)JK\(^{-1}\quad m_e\,:\,9.11\times 10^{-31}\,\)kg
\(\qquad e\,:\,1.60 \times 10^{-19}\,\)C\(\qquad \epsilon_0=8.85\times 10^{-12}\)C\(^2\)N\(^{-2}\)m\(^{-2}\)
電子の速さを\(\,v_e=\sqrt{\dfrac{k_BT}{m_e}}\,\)とすれば\(\quad v_e=\sqrt{\dfrac{k_BT}{m_e}}=\sqrt{\dfrac{1.38\times 10^{-23}T}{9.11\times 10^{-31}}}=3.89\times 10^3\sqrt{T}\,\,\)
\(\quad r_0=\dfrac{e^2}{2\pi\epsilon_0k_BT}=\dfrac{(1.60\times 10^{-19})^2}{2\pi\times 8.85\times 10^{-12}\times 1.38\times 10^{-23}T}=3.33\times 10^{-5}\dfrac{1}{T}\,\)m となる。
\(\blacklozenge\,\,\)どうもプラズマは、デバイ長やそれ以外の要素が大きいようで、単純な計算ではできないみたい。この検討はここまでとする。
◆ 参考文献
- 宇宙流体力学の基礎改訂版:福江 純、和田 桂一、梅村 雅之(2022)、日本評論社
- 天体物理学の基礎 I:観山 正見、野本 憲一、二間瀬 敏史(2010)、日本評論社
- 京大「エネルギー理工学設計演習・実験2(別冊)」
http://www.nucleng.kyoto-u.ac.jp/people/ikuji/edu/vac/app-A/mfp.html - 「14. 衝突頻度と平均自由行程」広島大学 学術情報リポジトリ
https://ir.lib.hiroshima-u.ac.jp/14983/files/40230
![]() |
![]() |
![]() |
![]() |
ロッシュポテンシャル
\(\def\bm{\boldsymbol}\)\(\def\di{\displaystyle}\)\(\def\ve{\varepsilon_0}\)\(\newcommand{\pdr}[2]{\dfrac{\partial #1}{\partial #2}}\)\(\newcommand{\ppdr}[2]{\dfrac{\partial^2 #1}{\partial #2}}\)
ロッシュポテンシャルとラグランジュ点
連星で、2つの成分星を\(\,M_1,\,M_2\,\)、星の重心間の距離を\(\,a\,\)とし、\(M_1\,\)星の重心を原点、2つの星の重心を結ぶ線を\(\,x\,\)軸、軌道面内に\(\,y\,\)軸とし、軌道は円軌道を仮定する。ここで任意の点\(\,(x,y,z)\,\)に働く力は2つの星の万有引力と遠心力の全ポテンシャル\(\,\Psi\,\)で表される。
\(\quad\Psi=-G\dfrac{M_1}{r_1}-G\dfrac{M_2}{r_2}-\dfrac{\Omega^2}{2}\left\{\left(x-\dfrac{M_2a}{M_1+M_2}\right)^2+y^2\right\}\qquad\cdots\,\)(1)
ここで
\(\quad r_1^2=x^2+y^2+z^2\quad,\,\,r_2^2=(a-x)^2+y^2+z^2\qquad\cdots\,\)(2)
ここで、\(\psi(x,y,z)\,\)を次のように定義する。
\(\quad\dfrac{GM_1}{a}\psi(x,y,z)=\Psi\qquad\cdots\,\)(3)
質量比\(\,q=M_2/M_1\,\)を使い、軌道面に限定して考えると
\(\quad\psi=-\dfrac{1}{r_1}-\dfrac{q}{r_2}-\dfrac{(1+q)(x^2+y^2)}{2}-\dfrac{q^2}{2(1+q)}\qquad\cdots\,\)(4)
これを 3D サーフェス表示と等高線で表示する。
また、ラグランジュ点は式(4)を微分して、\(\partial\psi/\partial x=0\,\)とした以下の式から求められる
\(\quad -\dfrac{x}{r_1^3}+\left(\dfrac{1-x}{r_2^3}-1\right)q+(1+q)x=0\qquad\cdots\,\)(5)
\(L_1\,\)を求めるには、\(r_1=x\,,\quad r_2=1-x\,\,\)を代入して通分すると
\(\quad -(1+q)x^5+(2+3q)x^4-(1+3q)x^3+x^2-2x+1=0\qquad\cdots\,\)(6)
\(L_2\,\)を求めるには、\(r_1=x\,,\quad r_2=x-1\,\,\)を代入して通分すると
\(\quad (1+q)x^5-(2+3q)x^4+(1+3q)x^3-(1+2q)x^2+2x-1=0\qquad\cdots\,\)(7)
(6)、(7)式をニュートン法で解く。
Python の実行例は 地球-月系の計算結果である。
◆ 参考文献
- W Ser型食連星はくちょう座V367星の 測光ならびに分光観測:小木美奈子(2012)、岡山理科大学
- LagrangePoint詳解-2: c Hiromu Nakagawa
import numpy as np
from matplotlib import pyplot as plt
import japanize_matplotlib
import warnings
import sys
# 警告非表示
warnings.simplefilter('ignore')
mx = 40
my = 30
pot = np.zeros((my,mx)) # ポテンシャルの配列
aa = np.zeros(6) # 5次方程式の係数
xx = np.linspace(-4.,4.,40)
yy = np.linspace(-2.,4.,30)
X, Y = np.meshgrid(xx,yy)
def newton(x,a): # 5次方程式をニュートン法で解く
x1 = x # 初期値
for i in range(20):
y1 = a[5]*x1**5 + a[4]*x1**4 + a[3]*x1**3 + a[2]*x1**2 + a[1]*x1 + a[0]
y2 = 5.*a[5]*x1**4 + 4.*a[4]*x1**3 + 3.*a[3]*x1**2 + 2.*a[2]*x1 + a[1]
if abs(y2) < 1.e-4:
return 0
dx = y1/y2
if abs(dx) < 1.e-4: # 差分が一定値以下なら
return x1
x1 = x1 -dx
return x1
def poten(x,y,q): # ロッシュポテンシャルを計算する
r1 = x**2 + y**2
r2 = (1.-x)**2 + y**2
zz = 1./np.sqrt(r1) + q/np.sqrt(r2) + .5*(1.+q)*r1 -q*x + .5*q**2/(1+q)
return -zz
q = float(input( '質量比( M2/M1) : '))
if q < 1.e-2 or q > 0.95:
sys.exit() # 入力エラーで終了
# ラグランジュ点の計算
aa[0] = 1
aa[1] = -2
aa[2] = 1
aa[3] = -1 - 3 * q
aa[4] = 2 + 3 * q
aa[5] = -1 - q
x1 = 0.5 # L1 の初期値
x2 = newton(x1,aa)
print('\n ---> L1 の座標:','{:.3f}'.format(x2),' (伴星のX座標を1)')
aa[0] = -1
aa[1] = 2
aa[2] = -1 - 2 * q
aa[3] = 1 + 3 * q
aa[4] = -2 - 3 * q
aa[5] = 1 + q
x1 = 1.2 # L2 の初期値
x2 = newton(x1,aa)
print('\n ---> L2 の座標:','{:.3f}'.format(x2),' (伴星のX座標を1)')
# ロッシュポテンシャルの計算
pot = poten(X, Y, q)
fig=plt.figure(figsize=(15, 10))
ax = plt.subplot(111, projection='3d')
ax.plot_surface(X, Y, pot, cmap='viridis')
ax.set_zlim(-30,0)
cs = ax.contourf(X,Y,pot,zdir='z',offset=-30, cmap='autumn')
plt.colorbar(cs)
plt.show()
質量比( M2/M1) : 0.0123 ---> L1 の座標: 0.849 (伴星のX座標を1) ---> L2 の座標: 1.168 (伴星のX座標を1)
![]() |
![]() |
積分その1
積分その1
【1】\(\,\,\displaystyle\int_0^{\infty}xe^{-ax}dx\quad\)の計算
\(\quad\displaystyle\int_0^{\infty}xe^{-ax}dx\quad\cdots\,\)(1)
まず、\(\,\displaystyle\int xe^{-ax}dx\,\,\)の不定積分から計算する。
部分積分で考えると
\(\quad\displaystyle\int xe^{-ax}dx=x(-ae^{-ax})+\dfrac{1}{a}\displaystyle\int e^{-ax}dx\)
\(\qquad=-axe^{-ax}+\dfrac{1}{a}\displaystyle\int e^{-ax}dx\quad\cdots\,\)(2)
(1)式の定積分を考える時、まず次の極限を計算する。
\(\quad\displaystyle\lim_{x\to\infty}\dfrac{x}{e^x}\quad\cdots\,\)(3)
\(\quad x\ge 0\,\,\)のとき\(\,\,e^x\ge\dfrac{x^2}{2}\,\,\)を考える。
\(\quad f(x)=e^x-\dfrac{x^2}{2}\quad\cdots\,\)(4)
とすると、\(x\ge 0\quad\)のとき、\(f'(x)\gt 0\,\,,\quad f^{\prime\prime}\gt 0\quad\)なので、\(f(x)\gt 0\quad\)である。よって
\(\,\,e^x\ge\dfrac{x^2}{2}\,\,\)が成り立つ。そこで
\(\quad 0\ge \dfrac{x}{e^x}\ge \dfrac{x}{\frac{x^2}{2}}\quad\)つまり\(\quad 0\ge \dfrac{x}{e^x}\ge \dfrac{2}{x}\quad\)となり、
\(\quad\displaystyle\lim_{x\to\infty}\dfrac{2}{x}=0\quad\)なので\(\quad\displaystyle\lim_{x\to\infty}\dfrac{x}{e^x}=0\quad\cdots\,\)(5)
となる。
この結果 (2)式を\(\,0\,\)から\(\,\infty\,\)まで積分すると
\(\quad\displaystyle\int_0^{\infty}xe^{-ax}dx=\left[-axe^{-ax}\right]_0^{\infty}+\dfrac{1}{a}\displaystyle\int_0^{\infty}e^{-ax}dx\)
\(\qquad =0-0+\dfrac{1}{a}\left[-\dfrac{1}{a}e^{-ax}\right]_0^{\infty}=\dfrac{1}{a^2}\quad\cdots\,\)(6)
グラフは\(\,y=xe^{-x}\,\)である。\(x\,\to\,\infty\,\)まで積分すると斜線部分の面積は\(\,1\,\)となります。
![]() |
![]() |
![]() |
![]() |
ビリアル定理
1-1. クラウジウスのビリアル定理
ビリアル定理は 19世紀に Clausius によって考案された。
「系の平均活力は、その(平均)ビリアル(の大きさ)に等しい。」
ここで登場する「活力」(”vis viva”)は、今日の運動エネルギーに相当します。 そもそも「活力」は古典力学の草創期にライプニッツが導入した量で、今日の運動エネルギーの2倍に相当する量 \(\,mv^2\,\)ですが、クラウジウスは \(\,mv^2/2\,\)とちょうど運動エネルギーと同じ量として用いています。ビリアル virial はラテン語のvis(「力」)からクラウジウスが作った造語です。ここでは、(力)×(位置ベクトル)をビリアルと呼ぶことにします。
距離の逆2乗則に従う重力やクーロン力の中心力で相互作用しあっている多体系では、長時間平均した運動エネルギー\(\,\langle\,T\,\rangle\,\)と平均の全ポテンシャルエネルギー\(\,\langle\,V\,\rangle\,\)は次の関係式を満たす。
\(\quad 2\langle\,T\,\rangle+\langle\,V\,\rangle=0\quad\cdots\,\)(1)
1-2. ビリアル定理の古典力学的証明
ここで位置ベクトル\(\,\vec{r}\,\)と運動量\(\,\vec{p}\,\)の内積の総和を以下のように考える。
\(\quad G=\displaystyle\sum_i\vec{p}_i\,\cdot\,\vec{r}_i\quad\cdot\,\)(2)
(2)式を時間\(\,t\,\)で微分する。
\(\quad\dfrac{dG}{dt}=\displaystyle\sum_i\vec{p}_i\cdot\dfrac{d\vec{r}_i}{dt}+\displaystyle\sum_i\dfrac{d\vec{p}_i}{dt}\cdot\vec{r}_i\)
\(\qquad\quad=\displaystyle\sum_i\,m_i(\vec{v}_i)^2+\displaystyle\sum_i\vec{F}_i\cdot\vec{r}_i\)
\(\qquad\quad=\displaystyle\sum_i2\times\dfrac{1}{2}m_i(\vec{v}_i)^2+\displaystyle\sum_i\vec{F}_i\cdot\vec{r}_i\)
\(\qquad\quad=2T+\displaystyle\sum_i\vec{F}_i\cdot\vec{r}_i\quad\cdots\,\)(3)
ここでは、次の関係式を使用している。
\(\quad \vec{p}_i=m_i\vec{v}_i\,,\qquad\vec{v}_i=\dfrac{d\vec{r}_i}{dt}\,,\qquad\vec{F}_i=\dfrac{d\vec{p}_i}{dt}\quad\cdots\,\)(4)
(3)式の両辺を\(\,0\,\)から時間\(\,t\,\)の範囲で積分して\(\,t\,\)で割り、\(t\to\infty\,\)の極限をとって長時間平均する。すると粒子が動き得る範囲は有限なので\(\,G\,\)も有限だから、左辺は\(\,0\,\)に収束する。
\(\quad\displaystyle\lim_{t\to\infty}\dfrac{1}{t}\displaystyle\int_0^t\dfrac{dG}{dt}=\displaystyle\lim_{t\to\infty}\dfrac{G(t)-G(0)}{t}=0\quad\cdots\,\)(5)
したがって、
\(\quad 0=2\langle\,T\,\rangle+\left\langle\displaystyle\sum_i\vec{F}_i\cdot\vec{r}_i\right\rangle\quad\cdots\,\)(6)
つまり、ビリアル定理を得る。
\(\quad\langle\,T\,\rangle=-\dfrac{1}{2}\left\langle\displaystyle\sum_i\vec{F}_i\cdot\vec{r}_i\right\rangle\quad\cdots\,\)(7)
次にポテンシャルエネルギー\(\,V\,\)が中心力ポテンシャルで、粒子間の距離に反比例する形で、系のポテンシャル\(\,V\,\)が各粒子対の相互作用の和によって書き表される場合、以下のように表される。
\(\quad V(r_1,\cdots,r_N)=\displaystyle\sum_{i\lt j}\dfrac{a_{ij}}{|\vec{r}_i-\vec{r}_j|}\quad\cdots\,\,\)(8)
粒子\(\,i\,\)に作用する力の合計\(\,\vec{F}_i\,\)は、ポテンシャルを位置座標で微分することで以下のように表すことが出来る。
\(\quad\vec{F}_i=-\nabla_{r_i}V=\displaystyle\sum_{j\ne i}\dfrac{a_{ij}(\vec{r}_i-\vec{r}_j)}{|\vec{r}_i-\vec{r}_j|^3}=\displaystyle\sum_{j\ne i}\vec{F}_{ij}\quad\cdots\,\)(9)
ここで、\(\,\vec{F}_{ij}\,\)は、粒子\(\,j\,\)から粒子\(\,i\,\)に働く力である。
\(\quad\vec{F}_{ij}=a_{ij}\dfrac{\vec{r}_i-\vec{r}_j}{|\vec{r}_i-\vec{r}_j|^3}\quad\cdots\,\)(10)
この\(\,\vec{F}_{ij}\,\)を用いると、ビリアル(力と位置ベクトルの内積)は次のように表せる。
\(\quad\displaystyle\sum_i\vec{F}_i\cdot\vec{r}_i=\displaystyle\sum_{i,j(i\lt j)}\vec{F}_{ij}\cdot\vec{r}_i+\displaystyle\sum_{i,j(j\lt i)}\vec{F}_{ij}\cdot\vec{r}_i\)
\(\qquad\qquad\quad=\displaystyle\sum_{i,j(i\lt j)}\vec{F}_{ij}\cdot\vec{r}_i+\displaystyle\sum_{i,j(i\lt j)}\vec{F}_{ji}\cdot\vec{r}_j\qquad\because\,i\leftrightarrow j\)
\(\qquad\qquad\quad=\displaystyle\sum_{i,j(i\lt j)}\vec{F}_{ij}\cdot(\vec{r}_i-\vec{r}_j)\quad\because\,\vec{F}_{ji}\!=\!-\vec{F}_{ij}\,\,\cdots\,\)(11)
(10)式を(11)式に代入すると
\(\quad\displaystyle\sum_i\vec{F}_i\cdot\vec{r}_i=\!\displaystyle\sum_{i,j(i\lt j)}\dfrac{a_{ij}}{|\vec{r}_i-\vec{r}_j|}\)
\(\qquad\qquad\qquad=\displaystyle\sum_{i,j(i\lt j)}\!V_{ij}=\!\displaystyle\!\sum_{i,j(i\lt j)}\!V_{ji}\quad\cdots\,\)(12)
よって(6)式は以下のように表すことが出来る。
\(\quad 0=2\langle\,T\,\rangle+\left\langle\displaystyle\sum_i\vec{F}_i\cdot\vec{r}_i\right\rangle=2\langle\,T\,\rangle+\left\langle\displaystyle\sum_{i,j(i\lt j)}V_{ij}\right\rangle\)
\(\qquad\qquad=2\langle\,T\,\rangle+\langle\,V\,\rangle\quad\cdots\,\)(13)
![]() |
![]() |
![]() |
![]() |
ドメイン変更
tamami.tech というドメインを取得したけど、更新料が 3,400 円と高額なので、tamami8.net という 更新料 1,480 円のドメインに変更しました。
まあ変更と行っても、”*.tech” のドメインは 2021年の6月まで有効ですが、HTML のソースを順次書き換え中です。
最近、ロリポップではサーバーレンタルとドメイン取得をまとめてやると、更新料無料らしい。まあ、しょうがないと諦めています。
CO回転スペクトル
1. ミクロの世界の量子状態
1-1. 箱の中の粒子
一次元の粒子の波動方程式を考える。\(\,x\,\)軸上に、\(\,x=0\,\)から\(\,x=L\,\)までの長さが\(\,L\,\)の区間(\(\,0\ge\,x\ge\,L\,\))を考え、この区間の内部でだけ粒子が運動できるモノとする。このような状況の粒子を、1次元の箱の中の粒子と呼ぶ。
箱の外に出られないように、箱の外の位置エネルギーを\(\,\infty\,\)とする。また、箱の中では粒子に力が働かず自由に運動できるように、位置エネルギーを\(\,0\,\)とする。
解くべき方程式は、定常状態の波動方程式 \(\,\hat{H}\phi(x)=E\phi(x)\,\) であり、粒子の質量を\(\,m\,\)とするとハミルトニアン\(\,\hat{H}\,\)は以下のようになる。
\(\quad \hat{H}=-\dfrac{\hbar^2}{2m}\Delta+U(x)\quad\cdots\,\)(1.1)
変数は\(\,x\,\)のみの1次式なので、ラプラシアン\(\,\Delta\,\)は\(\,x\,\)の2次微分 \(\,d^2/dx^2\,\)に等しい。箱の外では位置エネルギーが無限大だから、\(\,\hat{H}\phi=E\phi\,\)が成り立つためには、\(\phi=0\,\)にならなければならない。
すなわち、箱の外の確率密度は\(\,0\,\)になり、粒子は箱の外には存在しない。次に箱の中では、位置エネルギーを\(\,U(x)=0\,\)としたので、波動方程式は次のような2次微分を含む方程式になる。
\(\quad -\dfrac{\hbar^2}{2m}\ \dfrac{d^2\phi}{dx^2}=E\phi\quad\cdots\,\)(1.2)
この式では\(\,\phi\,\)と\(\,x\,\)以外はすべて定数であり、次の形の基本的な微分方程式となっている
\(\quad \dfrac{d^2\phi}{dx^2}=-k^2\phi\quad\cdots\,\)(1.3)
ここで右辺の\(\,k\,\)は、求めるべきエネルギー\(\,E\,\)を含み、次式で与えられる定数である
\(\quad k=\left(\dfrac{2mE}{\hbar^2}\right)^{1/2}\quad\cdots\,\)(1.4)
式 (1.3)の一般解は指数関数または三角関数を含む次式の形に表される
\(\quad \phi(x)=ae^{ikx}+be^{-ikx}=A\sin kx+B\cos kx\quad\cdots\,\)(1.5)
波動関数の重要な性質として、連続性の条件があり、境界条件を考慮しなければならない。
箱の外では\(\,\phi=0\,\)であるから、箱の中(\(\,0\ge x\ge L\,)\)の\(\,\phi\,\)も箱の両端で値が\(\,0\,\)にならなくてはいけない。
そこで境界条件として \(\,\phi(0)=\phi(L)=0\,\) を使うと、\(\,x=0\,\)で\(\,\cos\,\)関数の部分は\(\,0\,\)にならないので、\(\sin\,\)の関数部分だけが残る。
また、\(\,x=L\,\)で\(\,\sin kx\,\)の値が\(\,0\,\)になるには \(\,kL=n\pi\,\)となり、その結果、波動関数とエネルギーがそれぞれ次のように決まる。
\(\quad\phi(x)=A\sin(n\pi x/L)\quad\cdots\,\)(1.6)
\(\quad E=\dfrac{(n\pi\hbar/L)^2}{2m}\quad\cdots\,\)(1.7)
ここで整数\(\,n\,\)は\(\,1,2,3\,\)などの自然数であり、これは量子数と呼ばれる。
(1.7) 式から、箱の中の粒子に許されるエネルギーは、図のように飛び飛びのエネルギー準位になる。
この結果は有限な空間に閉じ込められた粒子のエネルギーが飛び飛びになって量子化されることを示している。
ここで注目するべきことは、エネルギーが最低の基底状態、すなわち\(\,n=1\,\)の状態でも、エネルギーが\(\,0\,\)にならないことである。このエネルギーをゼロ点エネルギーという。
(1.7) 式で与えられる波動関数は図の中央に示したように、エネルギー準位を区別する量子数\(\,n\,\)が増えるにつれ上下の振動が激しくなる。波動関数の値が\(\,0\,\)になる位置を節(ふし)という。
節になる位置は、\(\phi\,\)でも\(\,\phi^2\,\)も同じであり、図の右に示した\(\,\phi^2\,\)の図から明らかなように、節の位置に粒子が観測される確率は\(\,0\,\)である。
1-2. 調和振動子
波動方程式の簡単な例として、左図のようなバネ(バネ定数を\(\,k\,\))につるされた質量\(\,m\,\)の錘(おもり)の振動を考える。位置エネルギーは、\(\,U(x)=\frac{1}{2}kx^2\,\) であり、釣り合った位置から\(\,x\,\)だけ錘がずれると、図の右に示した放物線に沿って上がって行き、常に中心に引き戻される力が働いて振動する。このような振動を調和振動子という。
この場合のハミルトニアンの運動エネルギーは箱の中の粒子の場合と同じになり、位置エネルギーは\(\,U(x)=\frac{1}{2}kx^2\,\) であるから、この調和振動子の波動方程式は次のように表される。
\(\left(-\dfrac{\hbar^2}{2m}\dfrac{d^2}{dx^2}+\dfrac{1}{2}kx^2\right)\Psi=E\Psi\quad\cdots\,\)(1.8)
この微分方程式は、次のように波動関数を仮定すると解くことができる。
\(\quad\Psi(x)=e^{-ax^2}\,\displaystyle\sum_{i=0}b_ix^i\quad\cdots\,\)(1.9)
解いた結果、次式に従うエネルギー準位が得られる。
\(\quad E_n=\left(n+\dfrac{1}{2}\right)h\nu\qquad(n=0,1,2,\cdots)\quad\cdots\,\)(1.10)
ここで、\(\,n=0,1,2,\cdots\,\) は振動の量子数(振動量子数)であり、\(h\,\)はプランク定数で、\(\nu\,\)はこの振動の固有振動数である。
調和振動子に許されるエネルギー準位は図の左に示すように等間隔になる。その間隔は、振動のエネルギー量子\(\,h\nu\,\)であり,固有振動数\(\,\nu\,\)はバネの力の定数が\(\,k\,\)で質量が\(\,m\,\)の古典力学の振動子と全く同じで、次式によって与えられる。
\(\quad\nu=\dfrac{1}{2\pi}\sqrt{\dfrac{k}{m}}\quad\cdots\,\)(1.11)
量子論の振動子には、古典力学と明らかに違った注目すべき特徴がある。エネルギーが最低の状態(基底状態)でも、そのエネルギーは\(\,0\,\)にならず、\(\,\dfrac{1}{2}h\nu\,\)というエネルギーをもって振動する。このエネルギーをゼロ点エネルギーといい、その時の振動をゼロ点振動という。
図の左に示したように、調和振動子の波動関数は中心の周りで振動子、遠くに行くと減衰する。値が\(\,0\,\)になる節の数は量子数\(\,n\,\)に等しく、エネルギーが高いほど多くなる。
1-3. 2粒子系の運動方程式
2個の粒子を含む場合すなわち2粒子系の波動方程式は、以下のように取り扱うと1個の粒子の問題になり簡単になる。図の左に示したように、質量\(\,m_1\,\)と\(\,m_2\,\)の2個の球が距離$r$で結ばれているとして、その運動を考える。
ここで重心を座標原点に固定すると、\(\,r_1\,\)と\(\,r_2\,\)の和が\(\,r\,\)になり、2個の粒子の相対運動は長さが\(\,r\,\)の棒の先端につけた1つの球の運動の問題に単純化される。そのとき先端の粒子の質量\(\,\mu\,\)は、2個の粒子の質量の逆数の和の逆数になる。
\(\quad\)換算質量\(\quad\mu=\dfrac{1}{\frac{1}{m_1}+\frac{1}{m_2}}\)
こうして2個の粒子の相対運動は、換算質量をもつ1つの粒子の運動として取り扱うことが出来る。したがって、その波動方程式は次のように、1つの粒子の方程式と同じになる。
\(\quad\left(-\dfrac{\hbar^2}{2\mu}\Delta+U\right)\Psi=E\Psi\quad\cdots\,\)(1.13)
このような問題では、原点からの距離\(\,r\,\)が特別に重要なので、直交座標系\(\,(x,y,z)\,\)の代わりに、図の右に示す3次元の極座標\(\,(r,\theta,\varphi\,)\,\)を変数にとる。すると体積要素\(\,dxdydz\,\)とラプラシアン\(\,\Delta\,\)は、以下のように表される。
\(\quad\)体積要素\(\quad dxdydz=r^2\sin\theta dr d\theta d\varphi\quad\cdots\,\)(1.14)
\(\quad\)ラプラシアン\(\quad \Delta=\dfrac{1}{r^2}\dfrac{\partial}{\partial r}\left(r^2\dfrac{\partial}{\partial r}\right)+\dfrac{1}{r^2}\Lambda\quad\cdots\,\)(1.15)
\(\quad\)ルジャンドリアン\(\quad \Lambda=\dfrac{1}{\sin\theta}\dfrac{\partial}{\partial\theta}\left(\sin\theta\dfrac{\partial}{\partial\theta}\right)+\dfrac{1}{\sin^2\theta}\dfrac{\partial^2}{\partial\varphi^2}\quad\)(1.16)
ここで、(1.15)式のラプラシアンに含まれる\(\,\Lambda\,\)は角度に依存する演算子であり、ルジャンドリアンと呼ばれる。
1-4. 剛体回転子
2量子系の波動方程式は、換算質量\(\,\mu\,\)を使うと原点から距離\(\,r\,\)のところに繋がれた1個の粒子の回転運動になる。距離\(\,r\,\)が一定なので\(\,r\,\)に関する微分は消えてしまい、また位置エネルギーも\(\,0\,\)とおいてよいので、その波動方程式は次のように単純な式になる。
\(\quad-\dfrac{\hbar^2}{2I}\Delta\Psi=E\Psi\quad\cdots\,\)(1.17)
ここで、次式で表される\(\,I\,\)は、重心を軸とする回転運動の慣性モーメントである。
\(\quad I=\mu r^2\quad\cdots\,\)(1.18)
(1.17)式の波動方程式は三角関数を組み合わせた式を仮定して解くことが出来、その結果、剛体回転子に許されるエネルギーは次式のように、整数値をとる回転量子数\(\,J\,\)によって飛び飛びの準位になる。
\(\quad E=\dfrac{J(J+1)\hbar^2}{2I}\qquad (J=0,1,2,\cdots)\quad\cdots\,\)(1.19)
剛体回転子の波動方程式は、数学でよく知られている球面調和関数\(\,Y_{J,m}(\theta,\varphi)\,\)というものになり、それは\(\,J\,\)のほかにもう一つの量子数\(\,m\,\)にも依存している。\(\,J\,\)が決まれば、それに従って\(\,m\,\)のとりうる範囲が決まる。\(\,m\,\)は\(\,J\,\)から\(\,J-1,J-2,\cdots,0,\cdots,-J\,\)まで、合計\(\,2J+1\,\)通りの値が可能である。
\(\quad Y_{00}(\theta,\varphi)=\left(\dfrac{1}{4\pi}\right)^{1/2}\quad\cdots\,\)(1.20)
\(\quad Y_{10}(\theta,\varphi)=\left(\dfrac{3}{4\pi}\right)^{1/2}\cos\theta\quad\cdots\,\)(1.21)
\(\quad Y_{1\pm1}(\theta,\varphi)=\mp\left(\dfrac{3}{8\pi}\right)^{1/2}\sin\theta\ e^{\pm i\varphi}\quad\cdots\,\)(1.22)
\(\quad Y_{20}(\theta,\varphi)=\left(\dfrac{5}{16\pi}\right)^{1/2}(3\cos^2\theta-1)(\quad\cdots\,\)(1.23)
\(\quad Y_{2\pm1}(\theta,\varphi)=\mp\left(\dfrac{15}{8\pi}\right)^{1/2}\sin\theta\ \cos\theta\ e^{\pm i\varphi}\quad\cdots\,\)(1.24)
\(\quad Y_{2\pm2}(\theta,\varphi)=\mp\left(\dfrac{15}{32\pi}\right)^{1/2}\sin^2\theta\ e^{\pm i2\varphi}\quad\)(1.25)
2.分子の回転
2-1. 回転運動のエネルギー
質量\(\,m\,\)、速度\(\,v\,\)の直線運動の大きさが運動量\(\,p=mv\,\)で、運動に伴うエネルギーは
\(\quad E=\dfrac{1}{2}mv^2=\dfrac{1}{2}\dfrac{p^2}{m}\quad\cdots\,\)(2.1)
である。一方、回転運動の大きさを表す量が角運動量であり、回転の軸から距離\(\,r\,\)にあって、運動量\(\,p\,\)で運動している物体の角運動量\(\,p_{\theta}\,\)は、
\(\quad p_{\theta}=r\times p\quad\cdots\,\)(2.2)
と表される。\(\,p_{\theta}\,\)は運動量を含むベクトルである。この角運動量を用いて、回転運動のエネルギーは
\(\quad E=\dfrac{1}{2}\dfrac{p_{\theta}^2}{I}\quad\cdots\,\)(2.3)
のように表される。式(2.1)と(2.3)を比較すると対応していることが分かる。
直線運動を行う物体で運動の持続性(慣性)に影響する量である質量に対応するのが、回転運動における\(\,I\,\)であり、慣性モーメントと呼ばれる。慣性モーメントは
\(\quad I=\sum_im_ir_i^2=\sum_im_i(x_i^2+y_i^2+z_i^2)\quad\)(2.4)
であるが、2原子分子では換算質量\(\,\mu\,\)と原子間の距離\(\,r\,\)を使うと、以下の簡単な式で表すことが出来る
\(\quad \mu=\dfrac{m_1\cdot m_2}{m_1+m_2}\quad\cdots\,\)(2.5)
\(\quad I=\mu r^2\quad\cdots\,\)(2.6)
式(2.3)には運動量\(\,p_{\theta}\,\)が含まれており、この式を量子の世界にもちこむには、その\(\,p_{\theta}\,\)を演算子に置き換えればよい。その結果
\(\quad p_{\theta}=\sqrt{J(J+1)}\hbar\quad\cdots\,\)(2.7)
と置き換えることが出来る。従って、式(2.3)の2原子分子の場合の量子論的な表現は式(1.19)のように
\(\quad E=\dfrac{J(J+1)\hbar^2}{2I}=BJ(J+1)\qquad (J=0,1,2,\cdots)\quad\cdots\,\)(2.8)
となる。ここで\(\,J\,\)は回転量子数であり、\(\,0\,\)を含む整数のみが許される。すなわち回転エネルギーは量子化されており、飛び飛びの準位のみが許されることになる。\(\,B\,\)は回転定数と呼ばれるもので次式で表される。
\(\quad B=\dfrac{h}{8\pi^2I}\quad\)周波数単位}\(\qquad \dfrac{h}{8\pi^2c_0I}\quad\)波数単位\(\quad\cdots\,\)(2.9)
2-2. 観測される回転スペクトル
回転運動だけに起因するスペクトルはマイクロ波領域から遠赤外線領域に現れる。図は期待状態のCO分子の吸収スペクトルである。規則的な間隔の吸収線からなる。式(2.8)より2原子分子の場合の回転エネルギー準位は下図のようになることが理論的に予測できる。
回転量子数の\(\,J=0,1,2,3,\cdots\,\)に対して、エネルギー準位は\(\, 0,2B,6B,12B,\cdots\,\)となる。従って隣接するエネルギー準位間の大きさは \(\,2B,4B,6B,\cdots\,\)である。これを横軸にエネルギーをとってプロットすると図の下に示したようにスペクトルの位置が\(\,2B\,\)の間隔で並ぶ。これは実測のスペクトルとぴったり一致する。
このようにして、光子を1つ吸収することによる(2原子分子の)回転遷移は、\(\,J\,\)の値が1つだけ変化する状態の間に起こることがわかる。吸収はエネルギーが増える方向への遷移であるが、エネルギーが減る方向の瀬には発光スペクトルとして観測され、実際同様なスペクトルを示すことがわかっている。すなわち回転遷移の選択則は\(\,\Delta J=\pm 1\,\)である。ここでは、この選択則を実測との対応から導いたが、理論的にも同じ結果が導かれる。
以上のようにして、実測スペクトルの間隔から回転定数\(\,B\,\)が決定できる。この\(\,B\,\)は(2.9)式に定義された量であるから、\(\,B\,\)がわかるということは慣性モーメント\(\,I\,\)が実験的に求められたことを意味する。これより\(\,B\to I\to r\,\)のようにして、分子の構造を与える情報\(\,r\,\)が得られることになる。
2原子分子であるNO分子についても同様な解析が可能であり、実測スペクトルから分子構造に関する情報を導くことができる。それらの結果を以下の表にまとめる。
\(\mu\) (kg) | 回転定数\(\,B\) (Hz) | \(B\,\)cm\(^{-1}\) | 結合距離 \(r\,\mathring{\mathrm{A}}\) | |
---|---|---|---|---|
\(^{12}\)CO | \(1.139\!\times\!10^{-26}\) | \(5.793\times 10^{10}\) | \(1.932\) | \(1.128\) |
\(^{13}\)CO | \(1.191\!\times\!10^{-26}\) | \(5.538\times 10^{10}\) | \(1.847\) | \(1.128\) |
NO | \(5.012\!\times\!10^{10}\) | \(1.672\) | \(1.151\) | |
COの場合、\(J\!=\!1\!-\!0\,\)の回転遷移では、\(E=BJ(J+1)\,\) より、\(115.9\)GHz となり、\(^{13}\)COでは、\(110.8\)GHzとなる。\(J\!=\!2\!-\!1\,\)の場合は \(\,6B-2B=4B\,\)なので、\(^{12}\)CO は\(231.7\)GHzで、\(^{13}\)COは\(\,221.5\,\)GHzになる。
2-3. 観測されたスペクトル
爆発的星形成銀河(starburst銀河)M82の中心数100pc領域における、赤外線から電波に至る波長領域のスペクトル。大別して、(1)周波数に対して、連続的に放射強度が変化している「連続波放射」と、(2)ある特定の周波数でのみ、強い放射を示す「スペクトル線放射」とに分けられる。
連続波放射は、主に波長の長い電波域におけるシンクロトロン放射(磁場の中を高速で運動している電子からの放射)と、サブミリ波~赤外線域にかけて、大きな山をつくっている成分、すなわち、星間ダスト(星間空間中に存在する、シリケイトなどの固体微粒子)からの熱放射の重ね合わせとしてよく説明でできる。
この他、ミリ波付近には自由・自由遷移放射(電離領域中を運動する自由電子からの放射)も見られる。一方、この波長帯のスペクトル線放射としては、原子からの再結合線(桃色)・微細構造線(赤)・超微細構造線、分子からの回転線(青)などがみられる。
![]() |
![]() |
![]() |
![]() |
ステファン・ボルツマンの法則-2
1.プランクの法則 (Planck’s law)
物理学における黒体から輻射(放射)される電磁波の分光放射輝度、もしくはエネルギー密度の波長分布に関する公式。プランクの公式とも呼ばれる。ある温度\(\,T\,\)における黒体からの電磁輻射の分光放射輝度を全波長領域において正しく説明することができる。1900年、ドイツの物理学者マックス・プランクによって導かれた。
プランクはこの法則の導出を考える中で、輻射場の振動子のエネルギーが、あるエネルギー素量(現在ではエネルギー量子と呼ばれている)\(\,\epsilon=h\nu\,\)の整数倍になっていると仮定した。このエネルギーの量子仮説(量子化)はその後の量子力学の幕開けに大きな影響を与えている。
プランクの放射式では、単位振動数\(\,\nu\,\)あたりの黒体輻射強度を以下のように表す。
\(\quad B_{\nu}(T)=\dfrac{2h}{c^2}\dfrac{\nu^3}{e^{h\nu/kT}-1}\quad\cdots\,\)(1)
なお、ここで\(\,c\,\)は光速、\(\,h\,\)はプランク定数、\(\,k\,\) はボルツマン定数である。このプランクの放射式に温度\(\,T\,\)に3000K\(\sim\)7000Kを代入してプロットしたのが下図である。
2.全黒体輻射強度 (Total blackbody radiation intensity)
このプランクの放射式を全振動数で積分したものを、全黒体輻射強度\(\,B(T)\,\) という。(1)式で、\(\,x=h\nu/kT\,\)とおくと、
\(\quad B_x(T)=\dfrac{2k^3T^3}{c^2h^2}\dfrac{x^3}{e^x-1}\quad\cdots\,\)(2)
となり、全振動数で積分すると
\(\quad \displaystyle\int_0^{\infty}B_{\nu}(T)d\nu=\displaystyle\int_0^{\infty}B_x(T)\,\dfrac{kT}{h}\,dx\quad\cdots\,\)(3)
\(\quad B(T)=\dfrac{2k^4T^4}{c^2h^3}\displaystyle\int_0^{\infty}\dfrac{x^3}{e^x-1}\,dx\quad\cdots\,\)(4)
(4)式の積分は\(\,\pi^4/15\,\)となるので、
全黒体輻射強度(放射輝度)[単位 Wm\(^{-2}\)sr\(^{-2}\)K\(^{-4}\)] は
\(\quad B(T)=\dfrac{2\pi^4k^4T^4}{15c^2h^3}=\dfrac{\sigma}{\pi}T^4\quad\cdots\,\)(5)
となる。
3.\(\,\,\,\,x^3/e^x-1\,\)の\(\,\,\)積分の計算
\(\quad\dfrac{x^3}{e^x-1}=\dfrac{x^3}{e^x}\dfrac{1}{1-e^{-x}}\quad\cdots\,\)(6)
(6)式の最後の項は公比\(\,e^{-x}\,\)の等比数列の和とすると
\(\quad\dfrac{1}{1-e^{-x}}=1+e^{-x}+e^{-2x}+\cdots\)
となるので
\(\quad\dfrac{x^3}{e^x-1}=x^3(e^{-x}+e^{-2x}+e^{-3x}+\cdots)\quad\cdots\,\)(7)
ここで\(\,x^pe^{-ax}\,\)の形の関数を積分するために、ガンマ関数を利用する。
ガンマ関数の定義から
\(\quad\varGamma(q)=\displaystyle\int_0^{\infty}y^{q-1}e^{-y}dy=(q-1)\,!\quad\cdots\,\)(8)
ここで、\(\,p=q-1\,\)とすると
\(\quad p\,!=\displaystyle\int_0^{\infty}y^pe^{-y}dy\quad\cdots\,\)(9)
\(\,\,y=ax\,\)とすると\(\,dy=a\,dx\,\)で、(9)式に代入すると
\(\quad p\,!=\displaystyle\int_0^{\infty}(ax)^pae^{-ax}dx=a^{p+1}\displaystyle\int_0^{\infty}x^pe^{-ax}dx\quad\cdots\,\)(10)
よって
\(\quad\displaystyle\int_0^{\infty}x^pe^{-ax}dx=\dfrac{p\,!}{a^{p+1}}\quad\cdots\,\)(11)
これを使って(7)式を積分すると
\(\quad\displaystyle\int_0^{\infty}\!\dfrac{x^3}{e^x-1}dx=\displaystyle\int_0^{\infty}\!x^3(e^{-x}+e^{-2x}+\cdots)dx=\dfrac{3!}{1^4}+\dfrac{3!}{2^4}+\cdots\quad\)
ここでゼータ関数より
\(\quad\zeta(4)=\displaystyle\sum_{n=1}^{\infty}\dfrac{1}{n^4}=\dfrac{1}{1^4}+\dfrac{1}{2^4}+\cdots=\dfrac{\pi^4}{90}\quad\cdots\,\)(12)
なので
\(\quad\displaystyle\int_0^{\infty}\dfrac{x^3}{e^x-1}dx=\dfrac{3!}{1^4}+\dfrac{3!}{2^4}+\cdots=3\,!\,\zeta(4)=\dfrac{\pi^4}{15}\quad\cdots\,\)(13)
4.放射エネルギー
全立体角に対する積分は
\(\quad\displaystyle\int\cos\theta\,d\Omega=\displaystyle\int_0^{2\pi}\displaystyle\int_0^{\frac{\pi}{2}}\cos\theta\,\sin\theta\,d\theta\,d\phi=\pi\quad\cdots\,\)(14)
よって黒体の単位面積当たりの放射エネルギー\(\,I\,\)は
\(\quad I=\pi\displaystyle\int B_{\nu}d\nu=\dfrac{2\pi^5k^4}{15c^2h^3}T^4=\sigma T^4\quad\cdots\,\)(15)
ここで、\(\,\sigma\,\)はステファン・ボルツマン係数である
\(\quad\sigma=\dfrac{2\pi^5k^4}{15c^2h^3}=5.67\times 10^{-8}\)W/m\(^2\)K\(^4\quad\cdots\,\)(16)
![]() |
![]() |
![]() |
![]() |
ステファン・ボルツマンの法則-1
1. あらゆる物体は光を放っている
電磁波というのは電荷を持った粒子が振動することで発生する。そして電荷を持った粒子は、どこからか飛んできた電磁波を受けると、その影響で振動する。そこらじゅうにある原子は内部に電荷を持っているのだから、何らかの方法で電磁波を放出し、また受け止めていると考えられる。
例えば太陽の光だって電磁波の一種だが、それを受け止めた物体は熱くなる。電磁波の形で遠くから届いたエネルギーが熱に変わるからだ。逆のことも起きている。キャンプファイヤーが終わった後で火を消すと辺りは真っ暗になるが、その真っ暗闇の中で、さっきまで焼かれていた石がぼんやりと赤黒く光っていることがある。そこらにある普通の石が、内部から光を放っているように見えて、少し不気味ではあるが、何だか綺麗だ。
他にも色んな例がある。ストーブの中で焼けた鉄板が赤く光っている。焚き火のときに炭が赤く光っている。ハロゲンヒーターが赤熱している。その光を浴びると、すぐに暖かく感じる。これは間にある空気を伝わって暖められるわけではない。光によって直接暖められるせいだ。このように、熱せられた物体から光や電磁波が出る現象を「熱放射」と呼ぶ。
では物体は何度くらいから光を出し始めるのだろう。実は目に見えないだけで、低い温度の物体でも電磁波を出しているらしい。体温の高い人が近くに来ると直接触らなくても暖かさを感じることがあるだろう。それはその人の体から、可視光線よりも低いエネルギーの光である赤外線が放射されているのだ。
2. ステファン・ボルツマンの法則
では物体は一体どれほどのエネルギーを熱放射の形で放出しているのだろうか?それを調べたのはヨーゼフ・ステファンという 19 世紀の科学者であり、実はあのボルツマンの師である。彼は他の学者たちが過去に行った幾つかの実験結果をまとめた上で、自らも白金の針金に電流を流して赤熱させる実験を行ったりして、熱放射のエネルギー\(\,K\,\)と物体の温度\(\,T\,\)との間に次のような関係があるのではないかという提案を行った。 (1879)
\(\quad K=\sigma T^4\quad\cdots\,\)(1)
この式は「ステファン・ボルツマンの法則」と呼ばれている。なぜボルツマンの名前が入っているかというと、後に、弟子のボルツマンが、当時の最新理論である電磁気学と熱力学を使って、この法則の根拠を説明することに成功したからである。 (1884)
にもかかわらず、この結果はしばらくの間認められることはなかった。多くの科学者がこの実験結果を確認しようとしたが、追試がうまくいかなかったのである。今ではその理由ははっきり分かっているのだが、この法則には少しばかり適用条件があって、当時はそのことがまだ良く分かっていなかったのである。ステファンの実験がうまく行ったのは偶然の幸運であろうと言われている。
3. 理論的導出
電磁波はエネルギーだけでなく運動量をも持つのだということを、電磁気学で学んだ。そこで、こんなことを考えてみよう。容器の壁が温度\(\,T\,\)の物質で出来ていて、そこから放出された電磁波で容器の中が一様に満たされており、そのエネルギー密度が\(\,u\,\)になっていたとする。
この仮定は次のような考えを根拠にしている。もし容器の壁から内側へ向けて電磁波が放出される一方だとすると、容器内部の電磁場のエネルギーは無限に増える一方となるだろう。しかしそうならないのは、容器の壁は電磁場を放出すると同時に吸収も行っていて、一定の状態で安定しているせいだと考えられる。容器の中の電磁場と容器の壁とは熱平衡に達しているというわけだ。電磁場に対しても、温度という概念が使えるということでもある。というわけで、\(\,u\,\)は\(\,T\,\)の関数だと言えるだろう。
さて、そのときの運動量密度\(\,w\,\)と\(\,u\,\)の間には\(\,u=cw\,\)という関係がある。
4. 電磁波のエネルギー密度
マックスウェルの方程式から次の波動方程式が導かれる
\(\quad\nabla^2\mathbf{B}-\dfrac{1}{c^2}\dfrac{\partial^2\mathbf{B}}{\partial t^2}=0\qquad\nabla^2\mathbf{E}-\dfrac{1}{c^2}\dfrac{\partial^2\mathbf{E}}{\partial t^2}=0\quad\cdots\,\)(2)
空間に領域Dをとり、その内部で電磁場が荷電粒子系に及ぼす単位体積当たりの力は
\(\quad\mathbf{f}=\rho\mathbf{E}+\mathbf{j}\times\mathbf{B}\quad\cdots\,\)(3)
であるので荷電粒子系の運動方程式は
\(\quad\dfrac{d\mathbf{P}_{\text{粒子}}}{dt}=\displaystyle\int_{\mathrm{D}}\mathbf{f}\,dv\quad\cdots\,\)(4)
となる。\(\mathbf{P}_{\text{粒子}}\,\)は粒子系の全運動量である。
ここで電場についてのガウスの法則と拡張されたアンペールの法則を使って\(\,\mathbf{f}\,\)を表すと
\(\mathbf{f}=\epsilon_0(\nabla\cdot\mathbf{E})\mathbf{E}-\epsilon_0\dfrac{\partial\mathbf{E}}{\partial t}\times\mathbf{B}+\dfrac{1}{\mu_0}(\nabla\times\mathbf{B})\times\mathbf{B}\quad\cdots\,\)(5)
これよりD内の電磁場の運動量は単位体積当たり以下の値を持つ
\(\quad\mathbf{p}_{\text{場}}\equiv\epsilon_0\mathbf{E}\times\mathbf{B}\quad\cdots\,\)(6)
ここで電束密度\(\,\mathbf{D}=\epsilon_0\mathbf{E}\)、磁束密度\(\,\mathbf{B}=\mu_0\mathbf{H}\,\)であり、真空中を伝搬する電磁場の電場\(\,\mathbf{E}\)と磁束密度\(\,\mathbf{B}\)の大きさの間には\(\,|\mathbf{E}|=c|\mathbf{B}|\)なる関係が成立する。
電場の大きさを\(\,E=|\mathbf{E}|\)とすると電磁場のエネルギー密度は
\(\quad u=\dfrac{1}{2}(\mathbf{E}\cdot\mathbf{D}+\mathbf{B}\cdot\mathbf{H})=\dfrac{1}{2}\left(\epsilon_0E^2+\dfrac{E^2}{c^2\mu_0}\right)\)
\(\qquad =\dfrac{1}{2}(\epsilon_0E^2+\epsilon_0E^2)=\epsilon_0E^2\quad\cdots\,\)(7)
一方、運動量密度\(\,w\,\)は
\(\quad w=|\mathbf{D}\times\mathbf{B}|=\epsilon_0E\cdot\dfrac{E}{c}=\dfrac{\epsilon_0E}{c}=\dfrac{u}{c}\quad\cdots\,\)(8)
5. 電磁波による圧力
ここで、容器の内部の電磁波を考える。電磁波は容器の壁に当たるたびに壁に運動量を与えるのだから、容器の壁を圧力\(\,p\,\)で押すことだろう。まずは、その\(\,p\,\)と\(\,u\,\)の関係を導く。
説明が簡単に済むように、一辺の長さが\(\,L\,\)の立方体の容器を考えてみる。そしてどの壁を考えても同じことではあるのだが、極座標の積分計算を分かりやすくする都合で、\(\,z\,\)軸の正の方向にある壁に与える圧力を考える。
容器内の全ての電磁波が一つの方向を向いているとする。その全運動量\(\,P\,\)は運動量密度\(\,w\,\)と体積\(\,V=L^3\,\)を掛けて、\(P=wL^3\,\)と表せる。
全ての電磁波の向かう一つの方向というのが、\(\,z\,\)軸に対して\(\,\theta\,\)の角度を持った方向だとすると、壁に当たって反射するたびに、合計で\(\,2P\cos\theta\,\)の運動量を与えることになる。
その頻度であるが、電磁波の速度の\(\,z\,\)軸方向成分が\(\,c\,\cos\theta\,\)であり、往復\(\,2L\,\)の距離を進むたびに再び同じ壁に当たるのだから、1秒間に\(\,c\,\cos\theta/2L\,\)回の衝突を行うことになる。
力というのは1 秒あたりの運動量変化のことであり、それを面積\(\,L^2\,\)で割れば圧力\(\,p\,\)が求まることになる。
\(\quad p=2P\cos\theta\times\dfrac{c\,\cos\theta}{2L}\div L^2\)
\(\qquad =2wL^3\cos\theta\times\dfrac{c\cos\theta}{2L^3}\)
\(\qquad =wc\,\cos^2\theta=u\,\cos^2\theta\quad\cdots\,\)(9)
という形になる。しかしこの結果は、全ての電磁波が一つの方向を向いていると仮定した上での話であり、実際にはあらゆる方向を向いているはずだから、全方向について平等に平均をしてやらないといけないだろう。その為には立体角の考えを使えば良さそうだ。自分が半径1の球の中心にいて、周囲をぐるっと見回していると想像してみよう。どの方向も等しいのだから、この球の表面積にわたって今の結果を積分してやればいい。平均を取るためにはそれを表面積\(\,4\pi\,\)で割ってやればいい。
\(\quad p=\dfrac{u}{4\pi}\displaystyle\int_0^{\pi}\cos^2\theta\sin\theta\,d\theta\displaystyle\int_0^{2\pi}d\phi\)
\(\qquad =\dfrac{u}{4\pi}\left[-\dfrac{1}{3}\cos^3\theta\right]_0^{\pi}\times 2\pi\)
\(\qquad =\dfrac{u}{3}\quad\cdots\)(10)
6. 熱力学的考察
可逆過程での内部エネルギー(\(\,U\,\))変化において以下の関係式がある。
\(\quad dU=T\,dS-p\,dV\quad\) これを変形して
\(\quad\left(\dfrac{\partial U}{\partial V}\right)_T=T\left(\dfrac{\partial p}{\partial T}\right)_V-p\quad\cdots\,\)(11)
という式を作ることが出来る。これに\(\,T\,\)一定の条件下で\(\,V\,\)を変化させ、マクスウェルの関係式を当てはめる。電磁波の場合、この左辺はエネルギー密度\(\,u\,\)になる。
圧力\(\,p\,\)に式(10)を代入すると式(11)は
\(\quad u=\dfrac{1}{3}\,T\,\dfrac{du}{dT}-\dfrac{u}{3}\\[4pt]\therefore\,4u=T\,\dfrac{du}{dT}\qquad\dfrac{du}{u}=4\,\dfrac{dT}{T}\\[4pt]\quad\log u=4\log T+C\qquad\therefore\,\,u=aT^4\quad\cdots\text{(12)}\)となり、形の上ではステファン・ボルツマンの法則と同じものが導かれたことになる。
しかし\(\,u\,\)というのは、物体と電磁場とが熱平衡にあるときの電磁場のエネルギー密度を示しているのであって、物体から放出されるエネルギー\(\,K\,\)とは少し違う概念である。
7. 概念の整理
\(\,K\,\)の意味を「物体の単位表面積から単位時間あたりに放出されるエネルギー」だと言えるようにまとめることが出来れば使いやすい法則になりそうである。
微小面積\(\,dS\,\)から放出される電磁波をイメージする。
その微小面の法線方向から\(\,\theta\,\)だけずれた方向へ向かって微小時間\(\,dt\,\)に出て行った電磁波は、図のような筒状体積の内部にあることになるだろう。その体積は\(\,c\,\cos\theta\,dtdS\,\)である。
しかしこの筒状体積とエネルギー密度\(\,u\,\)をかけるだけで、そのエネルギー量が求まるほど単純ではない。なぜなら、このエネルギー密度\(\,u\,\)というのは、あらゆる方向に向かう電磁波が均等に重なり合った合計のエネルギー密度を表しているからである。
今知りたいのは、\(\,\theta\,\)方向へ流れる電磁波のエネルギーだけなのだが、その割合は全方向に向かう電磁波のエネルギーの中でどれくらいあると言えるだろうか。いや、厳密に\(\,\theta\,\)方向に向かう電磁波だけに限定してしまうとそれは無に等しいのである。だからある程度の微小幅を許容して考えないといけない。
こういう時は立体角の考えを使う。半径1の球殻を考え、その表面積の比で表すのである。角度\(\,d\theta\,\)と\(\,d\phi\,\)によって張られる球殻上の微小表面積は\(\,sin\theta\,d\theta\,d\phi\,\)と表されるから、全球殻の面積\(\,4\pi\,\)と比べた割合は\(\,1/4\pi\,\sin\theta\,d\theta\,d\phi\,\)となるだろう。
先ほど考えた筒状領域を、厳密に一方向として考えるのではなく、角度\(\,d\theta\,\)と\(\,d\phi\,\)の微小角度の幅を許してやれば、その範囲を向いている電磁波のエネルギーは
\(\quad\dfrac{1}{4\pi}uc\,\cos\theta\sin\theta\,d\theta\,d\phi\,dt\,dS\quad\cdots\,\)(13)
と表せることになる。この式(13)を微小面より上にある全ての方向について積分してやれば、微小面\(\,dS\,\)から微小時間\(\,dt\,\)内にあらゆる方向へ出て行くエネルギーを表せる。
\(\quad\dfrac{1}{4\pi}uc\,dt\,dS\displaystyle\int_0^{\pi/2}\cos\theta\sin\theta\,d\theta\displaystyle\int_0^{2\pi}d\phi\\[3pt]\qquad=\dfrac{1}{2}uc\,dt\,dS\displaystyle\int_0^{\pi/2}\dfrac{1}{2}\sin 2\theta\,d\theta\\[3pt]\quad =\dfrac{1}{2}uc\,dt\,dS\left[-\dfrac{1}{4}\cos 2\theta\right]_0^{\pi/2}\\[3pt]\quad =\dfrac{1}{4}uc\,dt\,dS\quad\cdots\,\text{(14)}\)
意味を考えれば今計算したのは\(\,K\,dt\,dS\,\)のことでもある。つまり\(\,K\,\)というのは、
\(\quad K=\dfrac{c}{4}u\quad\cdots\,\)(15)
に相当することが言えるのである。ここに先ほどの\(\,u=aT^4\,\)を代入すれば、
\(\quad K=\dfrac{ca}{4}\,T^4=\,\sigma T^4\quad\cdots\,\)(16)
という形が成り立っていることが言えるのである。普通はステファン・ボルツマンの法則と言えば、ここで使ったような意味での\(\,K\,\)と温度\(\,T\,\)の関係のことであり、「ステファン・ボルツマンの定数」として良く知られた\(\,\sigma\,\)もその意味に合うような数値が書かれていることが多い。
![]() |
![]() |
![]() |
![]() |
レーン=エムデン方程式(3)
レーン=エムデン方程式は以下の形で表される
\(\quad\dfrac{1}{\xi^2}\dfrac{d}{d\xi}\left(\xi^2\dfrac{d\theta}{d\xi}\right)=-\theta^N\qquad\cdots\,\)(1)
境界条件は、\(\xi\to 0\,\)のとき\(\,\theta=1\,,\,\dfrac{d\theta}{d\xi}=0\,\)である。
[1]中心近傍での扱い
この方程式は\(\xi=0\,\)が確定特異点になっている。この場合、特異点近傍でローラン級数展開を考える。\(\xi=0\,\)の周りの形式解を以下のようにする。
\(\quad\theta(\xi)=1+a_2\xi^2+a_4\xi^4+\cdots\qquad\)(2)
式(2)を式(1)の左辺に代入する。まず、\(\xi\,\)で微分すると
\(\quad\dfrac{d\theta}{d\xi}=2a_2\xi+4a_4\xi^3+\cdots\qquad\)(3)
さらに\(\xi^2\,\)をかけて微分すると
\(\quad\dfrac{d}{d\xi}\left(\xi^2\dfrac{d\theta}{d\xi}\right)=6a_2\xi^2+20a_4\xi^4+\cdots\qquad\)(4)
\(\therefore\,\,\)左辺\(\,\,=6a_2+20a_4\xi^2+\cdots\qquad\)(5)
右辺は、まず\(B_2=a_2\xi^2\,,\,B_4=a_4\xi^4\,,\cdots\,\)として、
\(\quad(1+B_2+B_4+\cdots)^N\,\)の展開を考える。
最初の\(\,1\,\)の項は\(\,1^NB_2^0B_4^0\cdots\,\)であるので、
係数は\(\,\dfrac{N!}{N!0!0!\cdots}\,\)となる。
次の項は\(1^{N-1}B_2^1B_4^0\cdots\,\)で、
係数は\(\,\dfrac{N!}{(N-1)!1!0!\cdots}\,\)となり、
次の\(1^{N-1}B_2^0B_4^1B_6^0\cdots\,\)の係数は
\(\,\dfrac{N!}{(N-1)!0!1!0!\cdots}\,\)となる。因って
\(\quad(1+B_2+B_4+\cdots)^N=1+NB_2+NB_4+\cdots\qquad\)(6)
右辺\(\,=-\theta^N=-1-Na_2\xi^2-Na_4\xi^4-\cdots\,\)となるので
\(\,\xi^m\,\)の係数を比較すると
\(\quad 6a_2=-1\,,\qquad 20a_2=-Na_4\qquad\)となり
\(\quad a_2=-\dfrac{1}{6}\,,\qquad a_4=\dfrac{N}{120}\quad\)を得る。
よって中心近傍での解は
\(\quad \theta(\xi)\simeq 1-\dfrac{1}{6}\xi^2+\dfrac{N}{120}\xi^4\quad\cdots\,\)(7)
と表すことが出来る。
[2]数値計算
この方程式は(その2)で求めた解析解以外の解は数値計算に因らざる得ない。
\(\bigstar\,\)常微分方程式の解法
\(\quad\dfrac{dy_i}{dx}=f_i[x,y_i(x)]\quad i=1,2,\dots\quad\)(8)
を初期値問題として数値的に解く。
まず、\(x=x_0\,\)において初期値\(\,y_i(x_0)=y_{i0}\,\)が与えられたとき、\(x=x_0+h\,\)での\(\,y_i\,\)の値\(\,y_i(x_0+h)\,\)をもとめ、次のステップでその\(\,(x,y_i)\,\)を新たな初期値として処理を繰り返す。
\(x_0\,\)と\(\,x_0+h\,\)の間の中間の\(\,x\,\)における導関数、\(f(x,y)\,\)の関数値を何回か計算し、それらの適当な平均を用いることによって、高次のTaylor展開を達成するのが Runge-Kutta 法である。
\(\bigstar\,\)Runge-Kutta 法
\(\quad k_1=hf(x_0,y_0)\qquad\qquad\qquad\cdots\,\)(9)
\(\quad k_2=hf(x_0+\frac{h}{2},y_0+\frac{k_1}{2})\qquad\cdots\,\)(10)
\(\quad k_3=hf(x_0+\frac{h}{2},y_0+\frac{k_2}{2})\qquad\cdots\,\)(11)
\(\quad k_4=hf(x_0+h,y_0+k_3)\qquad\cdots\,\)(12)
\(\quad y(x_0+h)=y(x_0)+\frac{1}{6}[K_1+2(K_2+k_3)+k_4]\,\cdots\,\)(13)
まず点(\(\,x_i,y_i\,\))での微係数から\(\,y\,\)の増分\(\,k_1\,\)を推定する(式9)。次に(\(\,x_i,y_i\,\))と今得た点(\(\,x_i+\Delta x,y_i+k_1\,)\)の中間点での微係数から、これに平行な破線で示す直線によって増分\(\,k_2\,\)の推定を行う(式10)。これより第3の推定点\((\,x_i+\Delta x,y_i+k_3\,)\)を得る。第4の推定点まで求めたところで、これら4つの増分推定を\(\,1:2:2:1\,\)の重みを付けて平均化し、これを最終的な増分推定として\(\,y_{i+1}\,\)を求める(式13)。この式は Taylor展開の4次の精度を持つ。
\(\bigstar\,\) 初期値とプログラミング
解くべき式は(1)式だが、\(\xi=x\,,\,\theta=y\,\)と書き換えて見やすくする。
\(\quad\dfrac{1}{x^2}\dfrac{d}{dx}\left(x^2\dfrac{dy}{dx}\right)=-y^N\quad x\!=\!0\quad y\!=\!0\,,\,\dfrac{dy}{dx}=0\,\cdots\,\)(14)
この常微分方程式を連立1階微分方程式に変換する。
\(y1=y\,,\quad y2=x^2\dfrac{dy}{dx}\,\)とおくと
\(\quad \dfrac{dy1}{dx}=\dfrac{y2}{x^2}\qquad\cdots\,\)(15)
\(\quad \dfrac{dy2}{dx}=-x^2y1^N\qquad\cdots\,\)(16)
初期値としては\(\,xi=x=0\,\)でなく、演算の刻み\(\,\xi=x=h\,\)とし、\(y1\,\)は式(7)を使い
\(\quad y1=1-\dfrac{1}{6}h^2+\dfrac{N}{120}h^4\quad\cdots\,\)(17)
\(y2\,\)は、\(\xi^2\dfrac{d\theta}{d\xi}=y2\quad\)としたので
\(\quad y2=h^2\left(-\dfrac{1}{3}h+\dfrac{N}{30}h^3\right)=-\dfrac{1}{3}h^3+\dfrac{N}{120}h^4\quad\cdots\,\)(18)
で設定する。
計算は中心\(\,\,\xi=x=h\,\,\)から始めて、密度係数がゼロ\(\,\,\,\theta=y1=0\,\,\)となる表面までを積分する。また\(\,\,h=1.0\times10^{-6}\,\,\)とした。
[3]数値計算による結果
\(N\,\)を\(0\sim 4.5\,\)まで計算した結果をグラフと表に示す。
この表は、\(N=n\,\)でガス球の表面(圧力ゼロ\(\,\theta=0\,\))の時の\(\,\xi_n\,\)とそこでの\(\,-1/\xi\,(d\theta/d\xi)\,\)で、中心密度\(\rho_c\,\)を計算するためのパラメータである。
Nの値 | \(\xi_n\) | \(\left(-\frac{1}{\xi}\frac{d\theta}{d\xi}\right)_{\xi=\xi_n}\) |
\(0.0\) | \(2.549\) | \(3.291\times10^{-1}\) |
\(0.5\) | \(2.753\) | \(1.816\times10^{-1}\) |
\(1.0\) | \(3.142\) | \(1.013\times10^{-1}\) |
\(1.5\) | \(3.654\) | \(5.564\times10^{-2}\) |
\(2.0\) | \(4.353\) | \(2.923\times10^{-2}\) |
\(2.5\) | \(5.355\) | \(1.424\times10^{-2}\) |
\(3.0\) | \(6.897\) | \(6.152\times10^{-3}\) |
\(3.5\) | \(9.536\) | \(2.180\times10^{-3}\) |
\(4.0\) | \(14.972\) | \(5.356\times10^{-4}\) |
\(4.5\) | \(31.836\) | \(5.386\times10^{-5}\) |
\(\bigstar\,\)結果の見方
\(\xi\,\)は \(\,R=\alpha\xi=\left[\dfrac{(N+1)K}{4\pi G}\rho_c^{\frac{1}{N}-1}\right]^{1/2}\xi\,\cdots\)(21) である。
また中心から\(\,\xi\,\)より内側に含まれる質量\(\,m(\xi)\,\)は
\(\quad m(\xi)=\displaystyle\int_0^{\alpha\xi}4\pi\rho r^2dr\)
\(\qquad\quad =4\pi\alpha^3\rho_c\displaystyle\int_0^{\xi}\xi^2\theta^Nd\xi\)
\(\qquad\quad =-4\pi\alpha^3\rho_c\displaystyle\int_0^{\xi}\dfrac{d}{d\xi}\left(\xi^2\dfrac{d\theta}{d\xi}\right)d\xi\)
\(\qquad\quad =-4\pi\alpha^3\rho_c\xi^2\dfrac{d\theta}{d\xi}\quad\cdots\,\)(22)
で与えられる。したがって星の全質量\(\,M\,\)は
\(\quad M=-4\pi\alpha^3\rho_c\left(\xi^2\dfrac{d\theta}{d\xi}\right)_{\xi=\xi_n}\)
\(\qquad =-4\pi\left[\dfrac{(N+1)K}{4\pi G}\rho_c^{\frac{1}{N}-1}\right]^{3/2}\left(\xi^2\dfrac{d\theta}{d\xi}\right)_{\xi=\xi_n}\,\cdots\,\)(23)
ここで星の平均密度\(\,\bar{\rho}\,\)を考える。
\(\quad\bar{\rho}=\dfrac{M}{4\pi R^3/3}=\dfrac{3M}{4\pi R^3}\)
\(\qquad =-\dfrac{3}{4\pi\alpha^3\xi_n^3}\cdot4\pi\alpha^3\rho_c\left(\xi^2\dfrac{d\theta}{d\xi}\right)_{\xi=\xi_n}\)
\(\qquad =-\dfrac{3}{\xi_n}\rho_c\left(\dfrac{d\theta}{d\xi}\right)_{\xi=\xi_n}\quad\cdots\,\)(24)
\(\bigstar\,\)太陽の場合
これらの関係を太陽の場合で計算してみる。太陽の半径\(\,R=6.960\times10^8\)m 、質量\(\,M=1.9884\times10^{30}\)kg と、恒星の場合に適応できるポリトロープ指数\(\,N=3\,\,(\Gamma=\frac{4}{3}\,)\,\)の計算結果を用いる。
\(\quad\rho_c=\dfrac{M}{4\pi R^3}\left(-\dfrac{1}{\xi}\dfrac{d\theta}{d/xi}\right)_{\xi=\xi_n}^{-1}\quad\cdots\,\)(25)
\(\qquad=\dfrac{1.9884\times10^{30}}{4\pi\cdot (6.960\times10^8)^3\cdot6.152\times10^{-3}}\)
\(\qquad =7.629\times10^4\,\mathrm{kg}\cdot\mathrm{m}^{-3}\)
理科年表によると\(\,\rho_c=1.56\times10^5\,\mathrm{kg}\cdot\mathrm{m}^{-3}\,\)*となっており、計算値は約半分である。まあ、天文のざっくりモデルでの計算としては、まあ合っているといえる。
*[参考]理科年表の太陽の構造の部分
![]() |
![]() |
![]() |
![]() |
レーン=エムデン方程式(2)
その1で導出したレーン=エムデン方程式は
\(\quad\dfrac{1}{\xi^2}\dfrac{d}{d\xi}\left(\xi^2\dfrac{d\theta}{d\xi}\right)=-\theta^N\quad\cdots\,\)(1)
で、\(\,\theta\,\)は無次元化した密度、\(\,\xi\,\)は無次元化した半径である。中心での境界条件は、\(\xi\to 0\,\)のとき\(\,\theta=1\,,\,\dfrac{d\theta}{d\xi}=0\,\)である。
[1]\(N=0\,\)の時の解析解
式(1)を一般の\(\,N\,\)について解析解は得られないが、\(N=0\,,\,1\,,\,5\,\)の時は得ることが出来る。
なお、以降は見やすくするために\(\,\xi\,,\,\theta\,\)を\(\,x\,,\,y\,\)と記述する。式(1)は
\(\quad\dfrac{1}{x^2}\dfrac{d}{dx}\left(x^2\dfrac{dy}{dx}\right)=-x^N\quad\cdots\,\)(2)
となる。ここに\(\,N=0\,\)を代入すると
\(\quad\dfrac{1}{x^2}\dfrac{d}{dx}\left(x^2\dfrac{dy}{dx}\right)=-1\quad\cdots\,\)(3)
両辺に\(\,x^2\,\)をかけると
\(\quad\dfrac{d}{dx}\left(x^2\dfrac{dy}{dx}\right)=-x^2\quad\cdots\,\)(4)
この式の両辺を\(\,x\,\)で積分すると
\(\quad x^2\dfrac{dy}{dx}=-\dfrac{1}{3}x^3+C_1\,\cdots\,\)(5)
となる。\(C_1\,\)は積分定数。両辺を\(\,x^2\,\)で割ると
\(\quad\dfrac{dy}{dx}=-\dfrac{1}{3}x+\dfrac{C_1}{x^2}\quad\cdots\,\)(6)
再度\(\,x\,\)で積分すると
\(\quad y=-\dfrac{1}{6}x^2-\dfrac{C_1}{x}+C_2\quad\cdots\,\)(7)
ここで、\(\,x=0\,\)のとき\(\,y=1\,\)の条件を満たすためには、\(C_1=0\,\)で\(\,C_2=1\,\)とする必要がある。よって、\(N=0\,\)の時は
\(\quad y=1-\dfrac{1}{6}x^2\quad\cdots\,\)(8)
となる。
[2]\(N=1\,\)の時の解析解
式(2)に\(\,N=1\,\)を代入すると
\(\quad\dfrac{1}{x^2}\dfrac{d}{dx}\left(x^2\dfrac{dy}{dx}\right)+y=0\quad\cdots\,\)(9)
ここで、\(z=xy\,\)とすると
\(\quad\dfrac{dy}{dx}=\dfrac{d}{dx}\left(\dfrac{z}{x}\right)=\dfrac{1}{x}\dfrac{dz}{dx}-\dfrac{z}{x^2}\quad\cdots\,\)(10)
式(10)の両辺に\(\,x^2\,\)をかけると
\(\quad x^2\dfrac{dy}{dx}=x\dfrac{dz}{dx}-z\quad\cdots\,\)(11)
となり、これを\(\,x\,\)で微分すると
\(\quad\dfrac{d}{dx}\left(x^2\dfrac{dy}{dx}\right)=\dfrac{d}{dx}\left(x\dfrac{dz}{dx}-z\right)=x\dfrac{d^2z}{dx^2}\quad\cdots\,\)(12)
式(12)を式(9)に代入すると
\(\quad\dfrac{1}{x}\dfrac{d^2z}{dx^2}+y=\dfrac{d^2z}{dx^2}+z=0\quad\cdots\,\)(13)
という2階微分方程式になる。
この方程式の解は
\(\quad z=A\cos x+B\sin x\)
で、\(z=xy\,\)と変数を元に戻すと
\(\quad y=\dfrac{A\cos x}{x}+\dfrac{B\sin x}{x}\quad\cdots\,\)(14)
ここで、\(x\to 0\,\)の時、\(y\to 1\,\)なので、\(A=0\,,B=1\,\)となり
\(\quad y=\dfrac{\sin x}{x}\quad\cdots\,\)(15)
となる。
[3]\(N=5\,\)の時の解析解
式(2)に\(\,N=5\,\)を代入すると
\(\quad\dfrac{1}{x^2}\dfrac{d}{dx}\left(x^2\dfrac{dy}{dx}\right)=-y^5\quad\cdots\,\)(16)
ここで、特殊な変数変換を行う。
\(\quad x=e^{-t}\quad\cdots\,\)(17)
\(\quad y=\sqrt{\dfrac{1}{2}e^t}\,Z\quad\cdots\,\)(18)
式(17)より
\(\quad\dfrac{dx}{dt}=-x\quad\cdots\,\)(17)
また、\(Z\,\)の\(\,x\,\)微分は
\(\quad\dfrac{dZ}{dx}=\dfrac{dZ}{dt}\dfrac{dt}{dx}=-\dfrac{1}{x}\dfrac{dZ}{dt}\quad\cdots\,\)(18)
微分演算子には以下の関係がある
\(\quad\dfrac{d}{dx}=-\dfrac{1}{x}\dfrac{d}{dt}\quad\cdots\,\)(19)
まず、式(18)を\(\,x\,\)で微分すると
\(\quad\dfrac{dy}{dx}=-\dfrac{1}{2\sqrt{2}}x^{-\frac{3}{2}}Z+\dfrac{1}{\sqrt{2}}x^{-\frac{1}{2}}\dfrac{dZ}{dx}\)
\(\qquad=-\dfrac{1}{2\sqrt{2}}x^{-\frac{3}{2}}Z-\dfrac{1}{\sqrt{2}}x^{-\frac{3}{2}}\dfrac{dZ}{dt}\quad\cdots\,\)(20)
式(20)の両辺に\(\,x^2\,\)をかけると
\(\quad x^2\dfrac{dy}{dx}=-\dfrac{1}{2\sqrt{2}}x^{\frac{1}{2}}Z-\dfrac{1}{\sqrt{2}}x^{\frac{1}{2}}\dfrac{dZ}{dt}\quad\cdots\,\)(21)
この式(21)を\(\,x\,\)で微分する
\(\quad\dfrac{d}{dx}\left(x^2\dfrac{dy}{dx}\right)=-\dfrac{1}{4\sqrt{2}}x^{-\frac{1}{2}}Z-\dfrac{1}{2\sqrt{2}}x^{\frac{1}{2}}\dfrac{dZ}{dt}\)
\(\qquad\qquad-\dfrac{1}{2\sqrt{2}}x^{-\frac{1}{2}}\dfrac{dZ}{dt}-\dfrac{1}{\sqrt{2}}x^{\frac{1}{2}}\dfrac{d}{dx}\left(\dfrac{dZ}{dt}\right)\)
\(\qquad=-\dfrac{1}{4\sqrt{2}}x^{-\frac{1}{2}}Z+\dfrac{1}{\sqrt{2}}x^{-\frac{1}{2}}\dfrac{d^2Z}{dt^2}\quad\cdots\,\)(22)
ここで
\(\quad -x^2y^5=-x^2\dfrac{1}{4\sqrt{2}}x^{-\frac{5}{2}}Z^5=-\dfrac{1}{4\sqrt{2}}x^{-\frac{1}{2}}Z^5\quad\cdots\,\)(23)
なので式(22)(23)を(16)に代入すると
\(\quad-\dfrac{1}{4}Z+\dfrac{d^2Z}{dt}=-\dfrac{1}{4}Z^5\quad\)より
\(\quad\dfrac{d^2Z}{dt^2}=\dfrac{1}{4}Z(1-Z^4)\quad\cdots\,\)(24)
という形に書き直すことができる。式(24)の両辺に\(\dfrac{dZ}{dt}\,\)をかけると
\(\quad\left(\dfrac{dZ}{dt}\right)\left(\dfrac{d^2Z}{dt^2}\right)=\dfrac{1}{4}Z\left(\dfrac{dZ}{dt}\right)-\dfrac{1}{4}Z^5\left(\dfrac{dZ}{dt}\right)\,\cdots\,\)(25)
ここで左辺は以下のように変形でき
\(\quad\left(\dfrac{dZ}{dt}\right)\left(\dfrac{d^2Z}{dt^2}\right)=\dfrac{1}{2}\dfrac{d}{dt}\left[\left(\dfrac{dZ}{dt}\right)^2\right]\quad\cdots\,\)(26)
右辺の計算も
\(\quad Z\left(\dfrac{dZ}{dt}\right)=\dfrac{1}{2}\dfrac{d}{dt}\left(Z^2\right)\quad\cdots\,\)(27)
と変形できるので、式(25)は
\(\quad\dfrac{d}{dt}\left[\left(\dfrac{dZ}{dt}\right)^2\right]=\dfrac{1}{4}\dfrac{d}{dt}\left(Z^2\right)-\dfrac{1}{12}\dfrac{d}{dt}\left(Z^6\right)\quad\cdots\,\)(28)
\(\,t\,\)で積分すると
\(\quad\left(\dfrac{dZ}{dt}\right)^2=\dfrac{Z^2}{4}-\dfrac{Z^6}{12}+D\quad\cdots\,\)(29)
となる。\(\,D\,\)は積分定数である。\(x\to 0\,\)のとき、\(Z=0\,dZ/dt=0\,\)なので、\(D=0\,\)である。よって次の式が得られる。
\(\quad\dfrac{dZ}{Z\sqrt{1-\frac{1}{3}Z^4}}=\pm\dfrac{1}{2}dt\quad\cdots\,\)(30)
左辺の正負の符号は境界条件より負を選択する。
\(\quad\dfrac{dZ}{Z\sqrt{1-\frac{1}{3}Z^4}}=-\dfrac{1}{2}dt\quad\cdots\,\)(31)
この方程式を解くために、さらなる変数変換を行う。
\(\quad\dfrac{1}{3}Z^4=\sin^2\phi\quad\cdots\,\)(32)
この式の微分形は以下のようになる。
\(\quad\dfrac{4}{3}Z^3\,dZ=2\sin\phi\,\cos\phi\,d\phi\quad\cdots\,\)(33)
式(32)(33)を(31)の左辺に代入すると
\(\quad\dfrac{dZ}{Z\sqrt{1-\frac{1}{3}Z^4}}=\dfrac{1}{Z\sqrt{1-\sin^2\phi}}\dfrac{3}{4Z^3}2\sin\phi\,\cos\phi\,d\phi\)
\(\qquad=\dfrac{3}{2}\dfrac{1}{Z^4}\dfrac{\sin\phi\,\cos\phi\,d\phi}{\sqrt{1-\sin^2\phi}}=\dfrac{1}{2\sin^2\phi}\dfrac{\sin\phi\,\cos\phi}{\cos\phi}d\phi\)
\(\qquad=\dfrac{1}{2}\dfrac{1}{\sin\phi}\,d\phi=\dfrac{1}{2}\mathrm{cosec}\,\phi\,d\phi\quad\cdots\,\)(34)
より(31)の方程式は以下の形に変形できる。
\(\quad\mathrm{cosec}\,\phi=-dt\quad\cdots\,\)(35)
両辺を積分すると
\(\quad\ln\left|\tan\dfrac{\phi}{2}\right|=-t+C’\quad\cdots\,\)(36)
より
\(\quad\tan\dfrac{\phi}{2}=Ce^{-t}\quad\cdots\,\)(37)
ここで、式(32)から
\(\quad Z^4=3\sin^2\phi=\dfrac{12\tan^2\frac{\phi}{2}}{\left(1+\tan^2\frac{\phi}{2}\right)^2}\quad\cdots\,\)(38)
この式(38)に式(37)を代入すると
\(\quad Z^4=\dfrac{12C^2e^{-2t}}{(1+C^2e^{-2t})^2}\)
\(\quad Z=\pm\left[\dfrac{12C^2e^{-2t}}{(1+C^2e^{-2t})^2}\right]^{\frac{1}{4}}\quad\cdots\,\)(39)
ここで、\(Z\,\)は式(18)より正で、式(17)(18)を代入すると
\(\quad y=\left[\dfrac{3C^2}{(1+C^2e^{-2t})^2}\right]^{\frac{1}{4}}=\left[\dfrac{3C^2}{(1+C^2x^2)^2}\right]^{\frac{1}{4}}\quad\cdots\,\)(40)
\(x=0\,\)で\(\,y=0\,\)という条件なので、
\(\quad 1=(3C^2)^{\frac{1}{4}}\quad\)より\(\quad C^2=\dfrac{1}{3}\quad\)が得られ
\(\quad y=\left[\dfrac{1}{(1+\frac{1}{3}x^2)^2}\right]^{\frac{1}{4}}=\dfrac{1}{\sqrt{1+\frac{1}{3}x^2}}\quad\cdots\,\)(41)
となる。
![]() |
![]() |
![]() |
![]() |
これまでの経緯
やってみたかった事
・高校では『天文研』に所属して、太陽の黒点の観測や夏の奥多摩での観測合宿などやってそれなりに、天文してました。完成はしませんでしたが、放課後交代で反射鏡を磨いていました。その後、大学と生業は『化学』という天文とは多少距離がある分野で長らく過ごしていました。
・やってみたかったと言うほどではないですが、いつかは...的な思いがあったのは、確かです。
放送大学で
・一応、定年(60歳)後を見据えて50代から、放送大学の専科履修生として勉強していました。定年後に放送大学の「自然と環境コース」に全科履修生として入学しました。(2013年度)
・それまで、専科履修生で取得した単位があったので、3年からの編入としました。昔の大学の単位を認定してくれる制度はあったようですが、面倒だと考えて語学も含めて全ての単位を放送大学で取得することにしました。その夏には卒業研究の申請を提出しています。
・それまで学習センターは試験や面接授業でしか利用していませんでしたが、いろいろなサークルなどがあり、先輩諸兄が指導してくれる環境がありがたかったです。卒論はどうせなら\(\,\mathrm{\TeX}\,\)を使いたいと思っていたところ、サークルにて\(\,\mathrm{\TeX}\,\)の講座を開催してくれたので、渡りに船と勉強させてもらいました。
・自然科学系の卒研に必要な科目と言ったら、まず数学です。一応、工学部だったのですが、昔取った杵柄が全く役に立ちませんでした。その数学を親切に教えてくださる先輩がいて、初歩から(再?)勉強させてもらっています。
天文での卒研
・卒業研究は天文を選択しました。報告書(卒論)はもちろん\(\,\mathrm{\TeX}\,\)で記述しました。放送大学のWAKABAの(資料)の「卒業研究資料_2014年度_06_自然と環境」に掲載されています。テーマはKeplar望遠鏡の観測データを使用した「連星の黒点挙動」です。
\(\,\mathrm{\TeX}\,\)について
・そんな感じで、備考録では\(\,\mathrm{\TeX}\,\)を使える環境をということで、ロリポップのレンタルサーバーの上で WordPress を使っています。
参考までに MathJaxの表示設定はこちらです。
<script type=”text/javascript” async=”” src=”https://cdnjs.cloudflare.com/ajax/libs/mathjax/2.7.1/MathJax.js?config=TeX-AMS_CHTML”></script>
この方式でなくともプラグインして[mathjax](括弧は半角)と文頭で宣言すればよいらしいですが...今のところ、うまく行かないので、文頭にスクリプトを書いています。
ポリトロープについて
・恒星は天文の勉強の基本です。ポリトロープは恒星の入門コースとなっています。これから導かれるレーン=エムデン方程式は、恒星物理学の基本です。約6年前勉強した内容を忘れないようにということで、ここに残します。
・振り子は最近までちゃんとした計算を知りませんでした。そんなんで、勉強した楕円積分とともに、備考録の最初に残しておきます。
レーン=エムデン方程式(1)
[1]Lane Emden 方程式とは
・宇宙物理学や流体力学において、レーン=エムデン方程式(Lane-Emden equation)は、球対称な密度分布を示す力学平衡にある自己重力流体を記述する微分方程式である。名称は宇宙物理学者のジョナサン・ホーマー・レーンとロバート・エムデンに由来する。
・Jonathan Homer Lane(1819-1880,ニューヨーク州のゲネシー生まれ)エール大学でオルムステッド(Denison Olmsted)に師事し、数学と自然科学を学ぶ。1857 年にワシントン市内に特許事務所を開き、特許関係の顧問や特許弁護士などの業務を始めた。1860 年に南北戦争が勃発すると、ペンシルバニア州北部の兄弟の住む村に移り、冷却用ガス膨張装置を用いた極低温装置の開発や天文学にも興味をもち、ひとりで文献に当たりながら太陽内部構造の研究に思いを深めていた。
・太陽の進化についてはジョン・ハーシェルやケルビン卿の収縮説に興味をもっていたが、内部構造を数学的に扱っていない点がレーンには不満であった。そこで彼は太陽をガス球と仮定し、重力平衡と熱力学過程を基本にして内部構造に関する基本方程式を書き始めた。
「太陽の理論的温度」(Lane,J.H.1870,Amer.J.Sci.&Arts,50,57-74) 自己重力のもとで球状に保たれるガス体の全体的平衡状態を取り扱った歴史的論文である。彼は太陽内部を完全気体と仮定したが、これは大胆な発想であった。完全気体を考えると\(\,p=R\rho T\)、静力学的平衡は\(\,dp=-\dfrac{Gm(r)}{r^2}\rho\,dr\,\)、内部は断熱的に対流平衡となるので、内部構造は\(\,1-\left(\dfrac{\rho}{\rho_c}\right)^{\gamma-1}=\displaystyle\int_0^x\dfrac{\mu(x)}{x^2}dx\,\)にて与えられる。
[2]アウグスト・リターのガス球論
・August Ritter(1826-1908)1853年にゲッチンゲン大学卒業。アーヘン工業大学(現RWTH-Aachen)の構造力学と一般力学の担当教授。「大気の高さに関する考察とガス状天体の内部構造について」ガス球の平衡・振動・安定性の問題に取り組み、自己重力にあるポリトロープガス球の平衡状態についての先駆的研究。
・ポリトロープとは圧力と密度に\(\,pV^n=K\,\) または\(\,p=K\rho^{1+1/n}\,\)の関係を満たす過程のこと。
\(\,T/T_c=y\,,\quad r/a=x\,\)とおいて、重力平衡式と組み合わせると以下の基本微分方程式が得られる。
\(\quad\dfrac{d^2y}{dx^2}+\dfrac{2}{x}\dfrac{dy}{dx}+y^n=0\,\quad\cdots\,\)(1)
この式は ポリトロープガス球の内部構造に関する基本方程式として、エムデンやエディントンに大きな影響を与えた。
[3]エムデンのガス球論
・ロベルト・エムデン(Jacob Robert Emden,1862-1940)スイスの天体物理学および気象学の専門家。1907 年にミュンヘン工科大学の教授。「ガス球論」出版。リターの基本方程式もそのままの形で再度導かれており、通常、エムデンの微分方程式と呼ばれている。リターらと異なり、エムデンはガス球半径を有限の場合と無限の場合を併行して考察している。エムデンの基本的考え方はレーンやリターと変わっていないが、理論が整備され、応用例が多いことからガス球理論の集大成と呼ばれ、4 版を重ねて、星の内部構造論の基本的テキストとなって広く知られるようになった。なお、エムデンの甥、マーチン・シュワルツシルドは 1958 に著書「星の構造と進化」を著わしているが、この書で歴史的なマイルストーンの研究者として、順にエムデン、エディントン、チャンドラセカールの 3 人の名前を挙げている。
[4]星の力学的安定
・連続の式:
半径\(\,r\,\)、ガス密度は半径の関数とし、\(\,\rho(r)\,\) と表すと、その半径より内側に含まれる質量\(\,M_r\,\)は
\(\quad M_r=\displaystyle\int_0^r4\pi r’^2\rho(r)dr’\quad\cdots\,\)(2)
と書くことができる。
この式の両辺を\(\,r\,\)で微分して微分形式で書き表すと
\(\quad\dfrac{dM_r}{dr}=4\pi r^2\rho(r)\quad\cdots\,\)(3)
となる。この式は流体力学の連続の式(continuity equation)の一種である。
・静水圧平衡:
星は自分自身の重力による収縮しようとする力とガス圧による対抗する力がバランスして、力学平衡を保っている。
\(\,r\,\)を中心からの距離、\(\,P\,\)をそこでの圧力、\(\,G\,\)を重力定数として
\(\quad\dfrac{dP}{dr}=-\dfrac{GM_r}{r^2}\rho\quad\cdots\,\)(4)
と表せる。
連続の式(3)、静水圧平衡の式(4)を解くことで、天体の内部構造を明らかにすることができる。 しかしこの段階では、独立変数が半径\(\,r\,\)で、従属変数は密度\(\,\rho\)、圧力\(\,P\)、質量\(\,M_r\,\)の3つで方程式は2つしかない。ためもうひとつ変数を結びつける式が必要となる。
・ポリトロープ関係
ここで、前出のポリトロープ関係と呼ばれる関係式を導入する。
\(\quad P=K\rho^{1+\frac{1}{N}}\quad\cdots\,\)(5)
\(\,K\,,\,\,N\,\)は定数のパラメータであり、\(\,N\,\)のことをポリトロープ指数と呼ぶ。
[5]レーン=エムデン方程式の導出
天体内部の構造を決めるためには、式(3)、(4)、(5)を解けば良いことになる。まず式(4)から
\(\quad\dfrac{r^2}{\rho}\dfrac{dP}{dr}=-GM_r\quad\cdots\,\)(6)
となり、両辺を\(\,r\,\)で微分すると
\(\quad\dfrac{d}{dr}\left(\dfrac{r^2}{\rho}\dfrac{dP}{dr}\right)=-G\dfrac{dM_r}{dr}\quad\cdots\,\)(7)
となる。これを式(3)に代入すると
\(\quad\dfrac{d}{dr}\left(\dfrac{r^2}{\rho}\dfrac{dP}{dr}\right)=-4\pi Gr^2\rho\quad\cdots\,\)(8)
となる。また、式(5)を\(\,r\,\)で微分すると
\(\quad\dfrac{dP}{dr}=K\left(1+\dfrac{1}{N}\right)\rho^{\frac{1}{N}}\dfrac{d\rho}{dr}\quad\cdots\,\)(9)
が得られる。式(9)を式(8)に代入すると
\(\quad K\left(1+\dfrac{1}{N}\right)\dfrac{d}{dr}\left(r^2\rho^{\frac{1}{N}-1}\dfrac{d\rho}{dr}\right)=-4\pi Gr^2\rho\quad\cdots\,\)(10)
ここで密度\(\,\rho(r)\,\)を星の中心密度\(\,\rho_c\,\)として
\(\quad\rho=\rho_c\theta^N\quad\cdots\,\)(11)
の形で表すこととする。\(\theta\,\)は無次元化した密度。
両辺を\(\,r\,\)で微分すると
\(\quad\dfrac{d\rho}{dr}=\dfrac{d}{dr}\left(\rho_c\theta^N\right)=\rho_cN\theta^{N-1}\dfrac{d\theta}{dr}\quad\cdots\,\)(12)
また、式(11)から以下の関係が導かれる
\(\quad\rho^{\frac{1}{N}-1}=\rho_c^{\frac{1}{N}-1}\,\theta^{1-N}\quad\cdots\,\)(13)
式(10)に式(12)(13)を代入すると
\(\quad K\left(1+\dfrac{1}{N}\right)\dfrac{d}{dr}\left(r^2\rho_c^{\frac{1}{N}-1}\theta^{1-N}\rho_cN\theta^{N-1}\dfrac{d\theta}{dr}\right)\!=\!-4\pi Gr^2\rho_c\theta^N\quad\)
\(\quad\dfrac{K(N+1)\rho_c^{\frac{1}{N}-1}}{4\pi G}\dfrac{1}{r^2}\dfrac{d}{dr}\left(r^2\dfrac{d\theta}{dr}\right)=-\theta^N\quad\cdots\,\)(14)
また式(5)を星の中心にあてはめると
\(\quad P_c=K\rho_c^{1+\frac{1}{N}}\quad\cdots\,\)(15)
この式(15)を式(14)に代入して\(\,K\,\)を消去すると
\(\quad\dfrac{(N+1)P_c}{4\pi G\rho_c^2}\dfrac{1}{r^2}\dfrac{d}{dr}\left(r^2\dfrac{d\theta}{dr}\right)=-\theta^N\quad\cdots\,\)(16)
左辺の係数部分を以下のように定義すると
\(\quad\alpha\equiv\left[\dfrac{(N+1)K}{4\pi G}\rho_c^{\frac{1}{N}-1}\right]^{1/2}\quad\cdots\,\)(17)
\(\quad\dfrac{\alpha^2}{r^2}\dfrac{d}{dr}\left(r^2\dfrac{d\theta}{dr}\right)=-\theta^N\quad\quad\cdots\,\)(18)
ここで、半径を以下の形に無次元化する。
\(\quad\xi=\dfrac{r}{\alpha}\quad\cdots\,\)(19)
これより無次元化した密度\(\,\theta\,\)の\(\,r\,\)の微分は以下のようになる
\(\quad\dfrac{d\theta}{dr}=\dfrac{d\theta}{d\xi}\dfrac{d\xi}{dr}=\dfrac{1}{\alpha}\dfrac{d\theta}{d\xi}\quad\cdots\,\)(20)
微分演算子は以下の形になる
\(\quad\dfrac{d}{dr}=\dfrac{1}{\alpha}\dfrac{d}{d\xi}\quad\cdots\,\)(21)
式(18)にこれらを用いると
\(\quad\dfrac{1}{\xi^2}\dfrac{d}{d\xi}\left(\xi^2\dfrac{d\theta}{d\xi}\right)=-\theta^N\quad\cdots\,\)(22)
この式は無次元化した密度\(\,\theta\)、無次元化した半径\(\,\xi\,\)を結びつけるもので、
レーン=エムデン方程式と呼ばれている。
![]() |
![]() |
![]() |
![]() |
ポリトロープその2
天文においてポリトロープとは、圧力と密度がポリトロピック関係を満たし力学平衡にある球対称な流体である。
一般に、圧力\(\,P\,\)と密度\(\,\rho\,\)の間に
\(\quad P=K\rho^{1+\frac{1}{N}}\quad\cdots\,\) (1)
の関係があるとき、この関係式をポリトロピックと呼ぶ。ここで\(\,N\,\)はポリトロピック指数、\(\,K\,\)は比例定数である。
熱力学では \(\,1+\dfrac{1}{N}=\gamma\,\) であるが、天文では後述の Lane-Emden 方程式の解法のために、この形の表記を行う。
また、断熱圧縮率(\(\,\Gamma_1\,\))を以下のような形で表す。
\(\quad\Gamma_1\equiv\,\left(\dfrac{\partial\ln P}{\partial\ln\rho}\right)_S=\left(\dfrac{\partial\ln P}{\partial\ln\rho}\right)_T\gamma\quad\cdots\,\) (2)
理想気体の場合\(\, (\partial\ln P/\partial\ln\rho)_T=1\,\)となるので、
\(\quad\Gamma_1=\gamma\equiv\,\dfrac{C_P}{C_V}\quad\cdots\,\) (3)
ポリトロピック指数と星の種類
・\(N=0\,\):密度一定のガス
・\(N=0.5\sim 1\,\):中性子星
・\(N=\frac{3}{2}\,\):非相対論的ガス球 白色矮星、ガス惑星 RGの核
・\(N=3\sim3.5\,\):相対論的ガス 放射優勢の太陽などの恒星
・\(N=5\,\):半径が無限大の星
・\(N=\infty\,\):等温のガス球
[1]状態方程式
粒子がエネルギー\(\,\epsilon\,\)、化学ポテンシャル\(\,\mu\,\)に存在する確率は分布関数
\(\quad f(\epsilon)=\dfrac{1}{\exp[(\epsilon-\mu)/k_BT]\pm1}\quad\cdots\)(4)
符号\(\,\pm\,\)は+がフェルミ粒子、-がボース粒子に対応する。
高温あるいは低密度で粒子の存在期待値が小さい場合、量子力学の効果が無視出来て
\(\quad f(\epsilon)=e^{\mu/k_BT}\,e^{-\epsilon/k_BT}\,\)と近似出来る。
マックスウェル・ボルツマン分布が成り立つ条件は\(\,e^{\mu/k_BT}\ll1\,(\mu<0\,)\,\)であり、理想気体の状態方程式が数密度が以下を充たす場合である。\(A_m\,\)を粒子当たりの平均分子量とすると
\(\quad n=\dfrac{\rho}{A_mM_u}\ll\dfrac{g(2\pi mk_BT)^{3/2}}{h^3}\quad\cdots\,\)(5)
[2]光子ガス
光子はスピン1のボース粒子でボース-アインシュタイン統計に従う。エネルギーを\(\epsilon\)とすると、\(\epsilon=pc\,,\quad v=c\,\)で統計的重みは\(\,g=2\)、光子数は保存されないので\(\mu=0\)である。これより放射の圧力とエネルギー密度は
\(\quad U_r=3P_r=aT^4\quad\cdots\,\) (6) ステファンボルツマンの式
放射の圧力(\(P_r\))の寄与の比を \(\,1-\beta\equiv P_r/P\,\)と表し、気体の比熱比を\(\,\gamma\,\)とすると、断熱指数は
\(\quad\Gamma_1=\dfrac{4}{3}+\dfrac{\beta(4-3\beta)(3\gamma-4)}{3\beta+36(\gamma-1)(1-\beta)}\quad\cdots\,\)(7)
と表せる。ガス圧が優勢の\(\,\beta=1\,\)の時は\(\,\Gamma_1=\gamma\,\)で、
放射優勢の極限\(\,\beta\to 0\,\)で\(\,\Gamma_1=4/3\,\)となる。
[3]電子縮退気体
化学ポテンシャル\(\,\mu\,\)と熱エネルギー\(\,k_BT\,\)との比として縮退パラメーター\(\,\Phi\,\)を定義する
\(\quad\Phi=\mu/k_BT\quad\cdots\,\)(8)
(5)式より縮退による理想気体からのズレは\(\,\Phi\simeq 0\,\)で現れるが、低温・高密度で\(\,\Phi\to\infty\,\)の極限で分布関数が階段関数で近似できる場合を完全縮退という。このときの化学ポテンシャルの最大値をフェルミエネルギー\(\,\epsilon_F\,\)、運動量をフェルミ運動量\(\,p_F\,\)という。
ここでフェルミモーメント\(\,p_F/m_ec=x_F\,\)とおくと、電子の数密度、圧力、運動エネルギー密度はそれぞれ
\(\quad n_e=\dfrac{8\pi}{h^3}\displaystyle\int_0^{p_F}\!p^2dp=\dfrac{8\pi}{3}\left(\dfrac{mc}{h}\right)^3x_F^3=Bx_F^3\quad\cdots\,\)(9)
\(\quad P_e=\dfrac{8\pi c}{3h^3}\displaystyle\int_0^{p_F}\!\dfrac{p^4dp}{\sqrt{p^2+m_e^2c^2}}=8A\!\displaystyle\int_0^{x_F}\!\dfrac{x^4\,dx}{\sqrt{1+x^2}}\cdots\,\)(10)
\(\quad U_e=\dfrac{8\pi}{h^3}\displaystyle\int_0^{p_F}(\sqrt{p^2c^2+m_e^2c^4}-m_ec^2)p^2dp\)
\(\qquad=\dfrac{8\pi m_e^4c^5}{h^3}\displaystyle\int_0^{x_F}\!(\sqrt{1+x^2}-1)x^2dx=Ag(x_F)\quad\cdots\,\)(11)
ここで、\(\,A\,,B\,\)はそれぞれ
\(\quad A\equiv\dfrac{\pi m_e^4c^5}{3h^3}=6.002\times10^{22}\)dyn cm\(^{-2}\)
\(\quad B\equiv\dfrac{8\pi}{3}\left(\dfrac{mc}{h}\right)^3=9.739\times10^5\)g cm\(^{-3}\)
また
\(\quad f(x)=8\displaystyle\int_0^x\!\dfrac{y^4\,dy}{\sqrt{1+y^2}}=\!x(2x^2-3)\sqrt{1+x^2}+3\sinh^{-1}x\)
\(\quad g(x)=24\displaystyle\int_0^x(\sqrt{1+y^2}-1)y^2\,dy\)
\(\qquad=3\{2x(x^2+1)^{3/2}-x\sqrt{x^2+1}-\sinh^{-1}x\}-8x^3\)
\(\qquad=8x^3(\sqrt{1+x^2}-1)-f(x)\)
のように定義されます。
\(\quad U_e=Ag(x_F)\,,\quad P_e=Af(x_F)\quad\)また、
\(\quad P=(\gamma-1)U\quad\) なので、電子ガスにおいて
\(\quad \gamma=\dfrac{f(x_F)}{g(x_F)}+1\quad\cdots\,\)(12)
と表すことが出来る。
図のように、\(p_F/m_ec=x_F\ll 1\,\)では、\(\gamma\to 5/3\,\)となり、
相対論的極限の \(x_F\gg 1\,\)では、\(\gamma\to 4/3\,\)となる。
![]() |
![]() |
![]() |
![]() |
ポリトロープその1
ポリトロープとは、熱力学で使うときと天文で使うときで多少感じが違う。まずは、熱力から
● ポリトロープ変化
圧力\( P\)と比容積\( V\)が以下の関係で結ばれる変化をポリトロープ変化という(この変化は可逆変化)
\(\quad PV^n=C\quad n\): ポリトロープ指数
\(\qquad n=0\qquad P=C\quad\,\,\): 等圧変化
\(\qquad n=1\qquad PV=C\quad\): 等温変化
\(\qquad n=\gamma\qquad PV^{\gamma}=C\quad\): 断熱変化
\(\qquad n=\infty\quad\,\, V=C\qquad\): 等積変化
【1】気体の比熱
熱力学の第1法則は
\(\,U\,\)を内部エネルギー、\(\,Q\,\)を熱量、\(\,W\,\)を仕事とすると
\(\quad dU=d’W+d’Q\quad \cdots\)(1)
と表せる。ここで \( d’\,\)は準静的過程を表す。
\(\quad d’W=-PdV\quad \cdots\)(2)
\(\quad dU=-PdV+d’Q\quad \cdots\)(3)
体積が一定として、\(dT\,\)で割ると
\(\quad C_V=\dfrac{d’Q}{dT}=\left(\dfrac{\partial U}{\partial T}\right)_V\quad \cdots\)(4)
同様に、圧力が一定とすれば
\(\quad C_P=\left(\dfrac{\partial U}{\partial T}\right)_P+P\left(\dfrac{\partial V}{\partial T}\right)_P\quad \cdots\)(5)
また、(4)式を \(T\)で積分すると\(\,\,U=C_VT\quad \cdots\)(6)
1モルでは \(\,\,PV=\mathrm{R}T\quad \cdots\)(7)\(\,\,\)で
\(\quad P\left(\dfrac{\partial V}{\partial T}\right)_P=\mathrm{R}\quad \cdots\)(8) なので、(5)式は
\(\quad C_P=C_V+\mathrm{R}\quad \cdots \)(9)
【2】断熱変化
\(dU\,\)の微分を考える
\(\quad dU=\left(\dfrac{\partial U}{\partial T}\right)_VdT+\left(\dfrac{\partial U}{\partial V}\right)_TdV\quad \cdots\)(10)
これを(3)式に代入すると
\(\quad d’Q=\left(\dfrac{\partial U}{\partial T}\right)_VdT+\left[P+\left(\dfrac{\partial U}{\partial V}\right)_T\right]dV\quad \cdots\)(11)
\(P=\)一定として(11)式を\(dT\,\)で割ると
\(\quad C_P=C_V+\left[P+\left(\dfrac{\partial U}{\partial V}\right)_T\right]\left(\dfrac{\partial V}{\partial T}\right)_P\quad \cdots\)(12)
断熱変化では \(d’Q=0\,\)なので(11)式は
\(\quad \left(\dfrac{\partial U}{\partial T}\right)_VdT+\left[P+\left(\dfrac{\partial U}{\partial V}\right)_T\right]dV=0\quad \cdots\)(13)
理想気体では\(\,(U\,\)は温度だけの関数なので \(\left(\dfrac{\partial U}{\partial V}\right)_T=0\)
(13)式に(6)式と(7)式を代入すると
\(\quad C_VdT+\dfrac{\mathrm{R}T}{V}dV=0\quad \cdots\)(14)
比熱比を\(\,\gamma=\dfrac{C_P}{C_V}\quad \cdots\)(15)\(\quad\)として(14)式を\(C_VT\,\)で割ると
\(\quad \dfrac{dT}{T}+(\gamma-1)\dfrac{dV}{V}=0\quad \cdots\)(16)\(\quad\)なので
\(\quad \ln T+(\gamma-1)\ln V= \text{一定} \cdots\)(17)\(\quad\)よって
\(\quad TV^{\gamma-1}= \text{一定} \cdots\)(18)
式(7)より\(\quad T=\dfrac{PV}{\mathrm{R}}\quad\cdots\)(19)を代入すると
\(\quad PV^{\gamma}=\text{一定}\quad \cdots\)(20)「ポアソンの法則」が得られる。
常温での\(\,C_V\,\)は並進運動と回転運動の自由度の\(\frac{\mathrm{R}}{2}\,\)倍となる。
単原子分子では並進運動の自由度が3なので、\(\,C_V=\frac{3}{2}\mathrm{R}\)。
2原子分子では\(\theta\,,\,\phi\,\)の回転が加わり、自由度が5になるので、\(\,C_V=\frac{5}{2}\mathrm{R}\)となる。
3原子分子(直線形)では、並進の自由度と回転の自由度は、2原子分子と同様だが振動の自由度が4となる。この振動は常温での比熱に寄与しないので、\(\,C_V\,\)は同じとなる。
よって、単原子分子は \(\gamma=\frac{5}{3}\,\)となり、2原子分子と直線形の3原子分子は \(\,\gamma=\frac{7}{5}\,\)となる。
![]() |
![]() |
![]() |
最初のテスト
TeX test
$$ x^2=c$$
$$\nabla=\dfrac{\partial}{\partial x} $$