投稿者「takeshibaoh」のアーカイブ

電荷密度

$\def\bm{\boldsymbol}$$\def\di{\displaystyle}$$\def\ve{\varepsilon}$$\def\vd{\varDelta}$$\newcommand{\pdr}[2]{\dfrac{\partial #1}{\partial #2}}$$\newcommand{\ppdr}[2]{\dfrac{\partial^2 #1}{\partial #2}}$

1.クーロンの法則

電荷$\,q\,$と$\,Q\,$を持つ２個の点電荷の間の距離を$\,R\,$,それらの間に作用する力の大きさを$\,F\,$と書くと

$F=\dfrac{1}{4\pi\varepsilon_0}\dfrac{q\cdot Q}{R^2}\qquad\cdots\,$(1)

点電荷$\,q\,$の位置ベクトルを$\,\bm{z}_0\,$、$Q\,$を$\,\bm{z}\,$とすると
$R=|\bm{z}_0-\bm{z}|\,$なので、電荷$\,Q\,$が点電荷$\,q\,$に作用する力は

$F(\bm{z}_0)=\dfrac{q}{4\pi\varepsilon_0}\dfrac{Q(\bm{z}_0-\bm{z})}{|\bm{z}_0-\bm{z}|^3}\qquad\cdots\,$(2)

次に、$n\,$個の点電荷$\,Q_1,Q_2,\cdots,Q_n\,$がそれぞれ$\,\bm{z}_1,\bm{z}_2,\cdots,\bm{z}_n\,$
の位置にあるとき、これらが点電荷$\,q\,$に作用する力は、
それぞれの点電荷が作用する力のベクトルの和で与えられる。

$\bm{F}(\bm{z}_0)=\dfrac{q}{4\pi\ve_0}\cdot\di\sum_{i=1}^n\dfrac{Q_i(\bm{z}_0-\bm{z}_i)}{|\bm{z}_0-\bm{z}_i|^3}\qquad\cdots\,$(3)

これがオリジナルのクーロンの法則であるが、$q\,$がゼロの場合、力が働かない事になる。
マックスウェルは近接作用である場の概念を導入し、電場というものを考えた。
これは、電荷がその周りに電場を作り、その電場に電荷が置かれると力が働くと言うことである
つまり、電場$\,\bm{E}\,$(ベクトル量を導入すると (2)式は以下のように分割される。

$\bm{F}(\bm{z}_0)=q\bm{E}(\bm{x}=\bm{z}_0)\quad\cdots\,$(4)$\qquad\bm{E}(\bm{x})=\dfrac{Q}{4\pi\ve_0}\dfrac{\bm{x}-\bm{z}}{|\bm{x}-\bm{z}|^3}\quad\cdots\,$(5)

2.ガウスの法則

(1) 積分形

点電荷$\,Q\,$が閉曲面$\,S\,$に囲まれている。
閉曲面$\,S\,$上にとった微小曲面$\,dS\,$とその上の電磁場$\,\bm{E}(\bm{x})\,$を考える。
$\,dS\,$に立てた単位法線ベクトルを$\,\bm{n}(\bm{x})\,$とすると、
$\bm{E}(\bm{x})\,$の$\,\bm{n}(\bm{x})\,$方向の成分は

$\bm{E}(\bm{x})\cdot\bm{n}(\bm{x})=\dfrac{Q}{4\pi\ve_0R^2}\cos\theta\qquad\cdots\,$(6)

ここで$\,R\,$は点電荷$\,Q\,$から微小面$\,dS\,$までの距離であり、$\theta\,$は$\,\bm{E}(\bm{x})\,$と$\,\bm{n}(\bm{x})\,$
とのなす角である。微小面$\,dS\,$を点電荷$\,Q\,$からみた立体角$\,d\Omega\,$は

$\dfrac{dS\cos\theta}{d\Omega}=\dfrac{R^2}{1^2}\quad$より$\quad d\Omega=\dfrac{\cos\theta\,dS}{R^2}\qquad\cdots\,$(7)

で与えられる。これを式(6)に代入すると$\quad\bm{E}(\bm{x})\cdot\bm{n}(\bm{x})\,dS=\dfrac{Q}{4\pi\ve_0R^2}\cos\theta\,dS=\dfrac{Q}{4\pi\ve_0}d\Omega$
これを閉曲面$\,S\,$上にわたって面積分すると

$\di\int_S\bm{E}(\bm{x})\cdot\bm{n}(\bm{x})\,dS=\dfrac{Q}{4\pi\ve_0}\di\int d\Omega=\dfrac{Q}{4\pi\ve_0}\cdot 4\pi=\dfrac{Q}{\ve_0}\qquad\cdots\,$(8)

点電荷が閉曲面の外側にある場合、$Q\,$からみた同一立体角内に二つの相対する微小曲面$\,dS\,,\,dS’\,$があることになる。
閉曲面の分割は任意であるから$\,dS\ne dS’\,$でよいが、$Q\,$からみた立体角$\,d\Omega\,$が同じなので、相殺される。
$S\,$内に、点電荷$\,Q_1,Q_2,\cdots,Q_n\,$がある場合、全体の電場$\,\bm{E}(\bm{x})\,$の面積分は、
含まれる点電荷の和を$\,\varepsilon_0\,$で割ったものになる。

$\di\int_S\bm{E}(\bm{x})\cdot\bm{n}(\bm{x})\,dS=\dfrac{1}{\ve_0}(Q_1+Q_2+\cdots+Q_n)\qquad\cdots\,$(9)

これを積分型のガウスの法則という。

(2) 微分形

　体積が$\,\varDelta^3x\equiv\varDelta x\varDelta y\varDelta x\,$である立方体の箱を考え
この箱の中に正電荷があるとする。この電荷が作る電場は箱の
中から外向きに出ている。そこで、箱の表面$\,S\,$上の面積分を考える。
まず、$x\,$軸に直交する二つの面 A と B を考える。A 面上では

$\di\int_A\bm{E}\cdot\bm{n}\,dS=E_x(x+\varDelta x,y,z)\varDelta y\varDelta z$
ここで変数$\,y\,,\,z\,$では２次以上の微少量を無視している。B 面上では

$\di\int_A\bm{E}\cdot\bm{n}\,dS=-E_x(x,y,z)\varDelta y\varDelta z$
ここで負号がつくのは、B面上に立てた単位法線ベクトル$\,\bm{n}\,$の方向が$\,x\,$軸の負の
方向を向いているからである。これらを加え合わせると

$\di\int_{A+B}\bm{E}\cdot\bm{n}\,dS=\dfrac{E_x(x+\varDelta x,y,z)-E_x(x,y,z)}{\varDelta x}\cdot\varDelta x\varDelta y\varDelta z=\pdr{E_x(\bm(x)}{x}\varDelta^3 x$

となる。$y\,$軸と$\,z\,$軸に垂直な面についても同様の計算をし、これらを加えると

$\di\int_S\bm{E}\cdot\bm{n}\,dS=\left(\pdr{E_x(\bm{x})}{x}+\pdr{E_y(\bm{x})}{y}+\pdr{E_z(\bm{x})}{z}\right)\varDelta^3x\qquad\cdots\,$(10)

ここで$\,\mathrm{div}\,\bm{E}(\bm{x})\,$という量を定義する。

$\mathrm{div}\,\bm{E}(\bm{x})=\di\lim_{\varDelta^3x\to 0}\dfrac{1}{\varDelta^3 x}\di\int_S\bm{E}\cdot\bm{n}\,dS\qquad\cdots\,$(11)　より

$\mathrm{div}\,\bm{E}(\bm{x})=\pdr{E_x(\bm{x})}{x}+\pdr{E_y(\bm{x})}{y}+\pdr{E_z(\bm{x})}{z}\qquad\cdots\,$(12)

これをベクトル$\,\bm{E}(\bm{x})\,$の発散 (divergence) という。

ここで、閉曲面$\,S\,$によって囲まれる空間の体積を$\,V\,$とすると、ガウスの定理より

$\di\int_S\bm{E}(\bm{x})\cdot\bm{n}\,dS=\di\int_V\mathrm{div}\,\bm{E}(\bm{x})\,d^3x\qquad\cdots\,$(13)

ここで、式(9)の左辺は式(13)のように変形されるから、式(9)は

$\di\int_V\mathrm{div}\,\bm{E}(\bm{x})\,d^3x=\dfrac{1}{\ve_0}\di\int_V\rho(\bm{x})d^3x\quad$となる。

ここで積分領域$\,V\,$を任意の点$\,\bm{x}\,$の周りの微小領域$\,\varDelta^3x\,$にとれば
$\mathrm{div}\,\bm{E}(\bm{x})\varDelta^3x=\rho(\bm{x})\varDelta^3x/\ve_0\,$となり、任意の場所で

$\mathrm{div}\,\bm{E}(\bm{x})=\dfrac{1}{\ve_0}\rho(\bm{x})\qquad\cdots\,$(14)

これを微分形のガウスの法則という。

3.静電ポテンシャル

(1) ストークスの定理

任意のベクトル場を考えたとき、以下の式が成立するのがストークスの定理

$\di\int_C\bm{X}(\bm{x})\cdot d\bm{x}=\di\int_S\mathrm{rot}\,\bm{X}(\bm{x})\cdot\bm{n}(\bm{x})\,dS\qquad\cdots\,$(15)

左辺の線積分は図のような任意の閉曲線$\,C\,$について、
図のような方向への１周にわたる積分であり、右辺の積分は
閉曲線$\,C\,$によって囲まれた任意の曲面$\,S\,$上での面積分である。
また、$\mathrm{rot}\,\bm{X}(\bm{x})\,$とは、その成分が

$\mathrm{rot}\,\bm{X}(\bm{x})\equiv\left(\pdr{X_z}{y}-\pdr{X_y}{z},\,\pdr{X_x}{z}-\pdr{X_z}{x},\,\pdr{X_y}{x}-\pdr{X_x}{y}\right)\quad\cdots\,$(16)

ここで、まず２次元で閉曲面$\,C’\,$を考える。$\bm{X}(\bm{x})\,$を$\,\bm{Y}(x,y)\,$と書き換える。
すると線積分は、２次以上の微少量を無視すると

$\di\int_{C’}\bm{Y}(x,y)\cdot d\bm{x}=Y_x(x,y)\varDelta x+Y_y(x+\vd,y)\vd y+Y_y(x,y)(-\vd y)+Y_x(x,y+\vd y)(-\vd x)$
$\hspace{30mm}=\left(\pdr{Y_y(x,y)}{x}-\pdr{Y_x(x,y)}{y}\right)\vd x\vd y$
$\hspace{30mm}=(\mathrm{rot}\,\bm{Y}(x,y))_z\vd x\vd y\qquad\cdots\,$(17)

式(17)より

$\left(\pdr{Y_y(x,y)}{x}-\pdr{Y_x(x,y)}{y}\right)=\dfrac{1}{\vd x\vd y}\di\int_{C’}\bm{Y}(x,y)\cdot d\bm{x}$
このような線積分を各微小長方形に分割して行う。隣り合わせの線積分の方向が反対向きなので
打ち消し合うことになる。最終的に隣り合う相手のいない周辺上の線積分だけが残るので、

$\di\int_{C’}\bm{Y}(x,y)\cdot d\bm{x}=\di\sum_i(\mathrm{tot}\,\bm{Y}(x_i,y_i))_z\vd x_i\vd y_i$
$\hspace{32mm}=\di\int_{S’}(\mathrm{rot}\,\bm{Y}(x,y))_zdS’\qquad\cdots\,$(18)

これが２次元のストークスの定理である。
この定理から$\,\mathrm{rot}\,\bm{X}_x(\bm{x})=0\,$のとき、以下の式がなりたつ。

$quad\di\int_C\bm{X}(\bm{x})\cdot d\bm{x}=0\qquad\cdots\,$(19)
閉曲線$\,C\,$の周りの１周の線積分を図のように
P点から出発してＱ点にいたる２つの曲線$\,C_1\,$と$\,C_2\,$上の線積分に分解すると

$\di\int_{\mathrm{P},C_1}^{\mathrm{Q}}\bm{X}(\bm{x})\cdot d\bm{x}+\di\int_{\mathrm{Q},C_2}^{\mathrm{P}}\bm{X}(\bm{x})\cdot d\bm{x}$
$\hspace{10mm}=\di\int_{\mathrm{P},C_1}^{\mathrm{Q}}\bm{X}(\bm{x})\cdot d\bm{x}-\di\int_{\mathrm{P},C_2}^{\mathrm{Q}}\bm{X}(\bm{x})\cdot d\bm{x}=0$

従って、点Pから点Qまでの線積分
$\quad\di\int_{\mathrm{P}}^{\mathrm{Q}}\bm{X}(\bm{x})\cdot d\bm{x}\qquad\cdots\,$(20)$\qquad$は途中の道筋によらず一定の値を取る。

すなわち、２点間の線積分が途中の道筋によらず一定の値を取るための必要かつ十分な条件は
$\quad\mathrm{rot}\,\bm{X}(\bm{x})=0\quad\cdots\,$(21)$\quad$が成立することである。

(2) 静電ポテンシャル

静電場が式(21)を満たすか、検証する。
まず、単位体積中の電荷量、つまり電荷密度を$\,\rho(\bm{x}_i)\,$とすると、電場$\,\bm{E}(\bm{x})\,$は

$\bm{E}(\bm{x})=\dfrac{1}{4\pi\ve}\di\int_{-\infty}^{\infty}d^3x’\dfrac{(\bm{x}-\bm{x}’)\rho(\bm{x}’)}{|\bm{x}-\bm{x}’|^3}\qquad\cdots\,$(22)

ここで、まず次の微分の計算をする。

$\pdr{}{x}\dfrac{1}{|\bm{x}-\bm{x}’|^n}=-n\dfrac{x-x’}{|\bm{x}-\bm{x}’|^{n+2}}\qquad\cdots\,$(23)

次に$\,\mathrm{rot}\,\bm{E}(\bm{x})\,$の$\,x\,$成分を計算する。式(22) から、

$\pdr{E_x(\bm{x})}{y}-\pdr{E_y(\bm{x})}{x}$
$\hspace{10mm}=\dfrac{1}{4\pi\ve_0}\di\int_{-\infty}^{\infty}d^3x’\rho(\bm{x}’)\left[(z-z’)\pdr{}{y}\dfrac{1}{|\bm{x}-\bm{x}’|^3}-(y-y’)\pdr{}{x}\dfrac{1}{|\bm{x}-\bm{x}’|^3}\right]$
$\hspace{10mm}=-\dfrac{3}{4\pi\ve_0}\di\int_{-\infty}^{\infty}d^3x’\rho(\bm{x}’)\left[\dfrac{(z-z’)(y-y’)}{|\bm{x}-\bm{x}’|^3}-\dfrac{(y-y’)(z-z’)}{|\bm{x}-\bm{x}’|^3}\right]=0\qquad\cdots\,$(24)

$y\,$成分、$z\,$成分についても同様である。よって静電場において$\quad\mathrm{rot}\,\bm{E}(\bm{x})=0\qquad\cdots\,$(25)$\quad$が成り立つ。

式(14)と式(25)を満たしている静電場$\quad\bm{E}(\bm{x})\,$は、ポテンシャル関数$\,\phi(\bm{x})\,$の勾配(gradient)で表す事が出来る。

$\bm{E}(\bm{x})=-\mathrm{grad}\,\phi(\bm{x})\qquad\cdots\,$(26)

この関数$\,\phi(\bm{x})\,$を 静電ポテンシャル という。
式(22)と(23)から、以下の式が得られる。

$\bm{E}(\bm{x})=-\dfrac{1}{4\pi\ve_0}\,\mathrm{grad}_x\di\int_{-\infty}^{\infty}d^3x’\dfrac{\rho(\bm{x}’)}{|\bm{x}-\bm{x}’|}\qquad\cdots\,$(27)

ここで、$\mathrm{grad}_x\,$は変数$\,x,y,z\,$による微分を取ることを意味する。
式(26)と式(27)から

$\phi(\bm{x})=\dfrac{1}{4\pi\ve_0}\di\int_{-\infty}^{\infty}d^3x’\dfrac{\rho(\bm{x}’)}{|\bm{x}-\bm{x}’|}\qquad\cdots\,$(28)

が成り立つ。
$\star\,$ガウスの法則より、電場の発散は式(14)から$\quad\nabla\cdot\bm{E}(\bm{x})=\rho(\bm{x})/\ve_0\qquad\cdots\,$(29)$\quad$となる。

4.電荷のテスト

２次元のベクトル場を設定する

適当な座標に正負の電荷をおく

numpy.gradient 関数を使ってベクトル場を計算する

ベクトル場の流線を表示する（正電荷から負電荷へと向かうよう符号を逆にする）

ベクトル場から発散を計算し電荷密度分布を表示する

参考資料

電磁気学　砂川重信　（岩倉書店）　物理テキストシリーズ

演習しよう　電磁気学　（数理工学社）　ライブラリ　物理の演習しよう２

In [85]:

import numpy as np
import matplotlib.pyplot as plt
import japanize_matplotlib   # 日本語表示用

#グラフ範囲を指定
LX, LY = 2.2, 2.2
# ベクトル場のベクトルの間隔の変数を gw とする
gw = 0.05
nn=int(2*LX/gw)   # 本数
coords,ds = np.linspace(-LX,LX,nn,retstep=True) # 座標と刻み幅の等差数列
# k = 4Πε0
k = 1.
x, y= np.meshgrid(coords, coords)    # x,y メッシュグリッド作成

# 点電荷の座標 (x1,y1), 電荷 q1
x_1,y_1,q_1 = -1.1, 0.2, 1
# (x1,y1) にある電荷から、点 (x,y) までの距離
r_1=np.sqrt((x-x_1)**2+(y-y_1)**2) 
# 点電荷の座標 (x2,y2), 電荷 q2
x_2,y_2,q_2 = 1, -0.9, -1   
# (x2,y2) にある電荷から、点 (x,y) までの距離
r_2=np.sqrt((x-x_2)**2+(y-y_2)**2) 

# 静電ポテンシャル phi は
phi=1/k*(q_1/r_1+q_2/r_2)
# 電場ベクトル Ex , Ey
E_y,E_x =np.gradient(phi, ds)

# x,yの描画範囲の設定
fig=plt.figure()
plt.xlim([-LX,LX])
plt.ylim([-LY,LY])
# グラフ描画
plt.grid()
plt.plot(x_1,y_1,'o',color='b')
plt.text(x_1+0.1,y_1+0.1,r'$q_1=+1$',color='b')
plt.plot(x_2,y_2,'o',color='r')
plt.text(x_2+0.1,y_2+0.1,r'$q_2=-1$',color='r')

# 電場ベクトルをプロット　 (+) から (-) に向くようにベクトルの向きを反対に
plt.streamplot(x,y,-E_x,-E_y,linewidth=0.3)
# タイトル
plt.text(x_1-0.2,1.5,"電場ベクトル表示",fontsize=25,color='m')
# ファイル保存
fig.savefig("electric_field.png")
plt.show()

# 電荷密度の計算
div_E=E_x-E_y

# 電荷密度の等高線をプロット
# x,yの描画範囲の設定
fig, ax = plt.subplots()
plt.xlim([-LX,LX])
plt.ylim([-LY,LY])
# グラフ描画
plt.grid()
#plt.plot(x_1,y_1,'o',color='b')
plt.text(x_1+0.1,y_1+0.1,r'$q_1=+1$',color='b')
#plt.plot(x_2,y_2,'o',color='r')
plt.text(x_2+0.1,y_2+0.1,r'$q_2=-1$',color='r')

# 電荷分布
CS = ax.contourf(x, y, div_E)
#fig.colorbar(CS)
# タイトル
plt.text(x_1,1.5,"電荷密度分布",fontsize=25,color='m')
# ファイル保存
fig.savefig("charge_distribution-2.png")
plt.show()

No description has been provided for this image

In [ ]:

傾きの計算

Gradient-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}}$

グラジエントのテスト

2次元でのテスト

grad(z) $\equiv\left(\dfrac{\partial z}{\partial x},\dfrac{\partial z}{\partial y}\right)$

python の numpy の gradient関数は、中心差分であり、始点側は後退差分、終点側は前進差分としている。
source と destnation の要素数が同じである。

diff関数は前進差分で、要素数が１つ減る。

前進差分

$\dfrac{du}{dx}=\displaystyle\lim_{h\to\ 0}\dfrac{u(x+h)-u(x)}{h}$

中心進差分

$\dfrac{du}{dx}=\displaystyle\lim_{h\to\ 0}\dfrac{u(x+h)-u(x-h)}{2h}$

後退差分

$\dfrac{du}{dx}=\displaystyle\lim_{h\to\ 0}\dfrac{u(x)-u(x-h)}{h}$

テスト

平面座標から山の高さを返す関数を作り、 x, y, z 軸での標高線を描く

表面をワイヤーフレームで表示する

x , y のメッシュグリッドを作成し、そのメッシュグリッドで山の高さを計算する

ｚ軸からの等高線グラフに、微分データ(Gradient)のベクトルを表示する

山データから高さを y 軸に平行な線で区切ったデータでプロットする

微分データから高さを y 軸に平行な線で区切ったデータでプロットする

In?[82]:

import numpy as np
import pylab as plt
from mpl_toolkits.mplot3d import Axes3D

def yama(x1, y1):
    h = 15 - np.sqrt((x1+10)**4 +(y1+6)**2)*0.015-np.sqrt((x1-20)**2+(y1-15)**2)*0.12
    return h
MAX_R=25
MIN_R=-25
MESH=2.5
fig = plt.figure()
ax = fig.add_subplot(projection='3d')

x = np.arange(MIN_R, MAX_R, MESH)    # 1次元のリスト
y = np.arange(MIN_R, MAX_R, MESH)

print('◆　データ数    x : ',len(x),'   y : ',len(y))

X, Y = np.meshgrid(x, y)             # 2次元のリスト
Z = yama(X,Y)                        # 2次元のリストからの変数なので Z は２次元
Z[np.less(Z, 0)] = 0
dZ_dX, dZ_dY = np.gradient(Z)  # Z が２次元あので、微分すると２次元変数が２つとなる

ax.contour(X, Y, Z, zdir='z', offset=0, cmap=plt.cm.Reds)
ax.contour(X, Y, Z, zdir='x', offset=MIN_R, cmap=plt.cm.OrRd)
ax.contour(X, Y, Z, zdir='y', offset=MAX_R, cmap=plt.cm.Greens)
ax.set_zlim(0, 16)
ax.plot_wireframe(X, Y, Z ,linewidth=0.3)
plt.show()

print('\n 傾斜のベクトルは下から上ではなく、上から下になるように向きを逆にしている\n')

fig, ax = plt.subplots(figsize=(5,5))
ax.contour(X, Y, Z,cmap=plt.cm.Reds)
ax.quiver(X, Y, -dZ_dY, -dZ_dX, angles="xy", scale_units='xy', scale=0.7,color='b')
plt.show()

print('\n 地形をグラフとして表わす\n')

fig, ax = plt.subplots(figsize=(5,5))
plt.plot(X[0],Z[0],linewidth=0.3,color='r',label='x=-20')
plt.plot(X[2],Z[2],linewidth=0.3,color='g',label='x=-16')
plt.plot(X[4],Z[4],linewidth=0.3,color='b',label='x=-12')
plt.plot(X[6],Z[6],linewidth=0.3,color='y',label='x=-8')
plt.plot(X[8],Z[8],linewidth=0.3,color='m',label='x=-4')
plt.plot(X[10],Z[10],linewidth=0.3,color='c',label='x=0')
plt.plot(X[15],Z[15],linewidth=0.3,color='k',label='x=10')
plt.legend()
plt.show()

print('\n 傾きをグラフとして表わす\n')

fig, ax = plt.subplots(figsize=(5,5))
plt.ylim(-0.6,0.6)
plt.plot(X[0],dZ_dX[0],linewidth=0.4,color='r',label='x=-20')
plt.plot(X[8],dZ_dX[8],linewidth=0.4,color='g',label='x=-4')
plt.plot(X[12],dZ_dX[12],linewidth=0.4,color='b',label='x=4')
plt.plot(X[15],dZ_dX[15],linewidth=0.4,color='c',label='x=10')
plt.plot(X[19],dZ_dX[19],linewidth=0.4,color='m',label='x=18')
plt.legend()
plt.show()

◆　データ数    x :  20    y :  20

 傾斜のベクトルは下から上ではなく、上から下になるように向きを逆にしている

 地形をグラフとして表わす

 傾きをグラフとして表わす

In?[?]:

ああそうなんだ倶楽部・ひよこ句会（私の句掲載分）

私が参加させてもらっている「ああそうなんだ倶楽部・ひよこ句会」での、
私の句が選句されたものの一覧です。

全くのド素人句ですが、仲間内で選んでくれたものなの(一部は先生の選句です)で、ちょっと嬉しいです。

年月日	句
2025-01-19	顔合わせ火鉢囲んで餅膨る
2024-12-22	ひとり旅車窓に霞む遠枯野
2024-11-16_1	世を嘆くゆるりと見上ぐ冬北斗
2024-11-16_2	秋深し四駆も喘ぐ滝の音
2024-10-20	ため息と喧騒の中残暑なり
2024-09-29	いにしえの風に浸りて秋彼岸
2024-06-15_1	巣鴨駅伸びる線路に風薫る
2024-06-15_2	AIと珍問答の梅雨曇り
2024-06-15	ため息と喧騒の中残暑なり
2024-10-20	庭先で懐旧のとき風薫る
2024-03-16	客帰り後片付けに花の冷え
2023-12-16_1	冬晴れに地学巡検輝石あり
2023-12-16_2	立ち止まる冬三角の荘厳さ
2023-11-18	大あくびご祈祷待ちの七五三
2023-10-21	降りて行く上弦の月緩やかに
々	々
2023-09-16	天上にそびえる塔の曼珠沙華
2023-06-17	佳き仲間杯交わす初鰹
2023-04-16	ゆったりと登る階段復活祭
2023-02-18	春北斗夜の温もり汲み取るか
2022-11-19	流れ星掬い取れるか冬北斗
2022-10-15_1	神無月忘れぬ香り出雲蕎麦
2022-10-15_2	嗅覚の底に響くや焼き秋刀魚
2022-09-17_1	蟋蟀と煌めく星の奏かな
2022-09-17_2	世のしくみ思い巡らせ夜の長し
2022-05-21_1	迷い路ふと何処から薫風か
2022-05-21_2	風薫る緑のターフいざダービー
2022-05-21_3	風薫る漕ぐ脚軽く土手の花
2022-03-19	春雷か平和を希うウクライナ
2022-01-15	寅の尾を避けて行くのか去年今年
2021-12-5	澄み切って冬の三角荘厳な
2021-11-05	オリオンの彼方にかすむギャラクシー
2021-09-05_1	新蕎麦や香り楽しむ江戸の粋
2021-09-05_2	懐かしい仙石原の芒かな
2021-07-10	願い込め梅雨入り前の七福神
々	々
2021-05-31	風青し日本ダービー空駆ける
2021-04-30	パート先巣立ちの鳩の後始末
2021-03-27	花冷えに五輪の期待儚くも
2021-01-16	今年こそ皆で花見を待ち焦がれ
2020-12-19	オリオン座冷たく光るシリウスと
2020-10-04	マスク越し爽やか匂う秋の風

面の角度

1. 面の式

　点 $A, B, C$ がある。この３点で構成される面の式を求める。

ベクトル $\overrightarrow{AB}$ と　$\overrightarrow{AC}$ を考える。点 $A, B, C$ の座標をそれぞれ $(x_0,y_0,z_0)\,,\,(x_1,y_1,z_1)\,,\,(x_2,y_2,z_2)$ とすると、各々の成分は

$\overrightarrow{AB}=(x_1-x_0,y_1-y_0,z_1-z_0)^T=(x_{AB},y_{AB},z_{AB})^T$とおく。
$\overrightarrow{AC}=(x_2-x_0,y_2-y_0,z_2-z_0)^T=(x_{AC},y_{AC},z_{AC})^T$とおく。

ここでベクトル $\overrightarrow{AB}$ と $\overrightarrow{AC}$ の外積を考える。
$\overrightarrow{AB}\times\overrightarrow{AC}=\begin{pmatrix}
y_{AB}x_{AC}-z_{AB}y_{AC}\\
z_{AB}x_{AC}-x_{AB}z_{AC}\\
x_{AB}y_{AC}-y_{AB}x_{AC}\end{pmatrix}=\begin{pmatrix}P_x\\P_y\\P_z\end{pmatrix}\quad\cdots\,$(1)　法線ベクトル$\,\vec{\bm{n}}\,$と平行となる。

この法線ベクトルが点$\,A\,$を通るとすると
$P_x(x-x_0)+P_y(y-y_0)+P_z(z-z_0)=0\quad\cdots\,$(2)　が平面の式である

2. 面の最大傾斜

式(2)より
　$P_xx+P_yy+P_zz-(P_xx_0+P_yy_0+P_zz_0)=0$
　この( )の部分を$\,D\,$とおけば

　$P_xx+P_yy+P_zz-D=0\quad\cdots\,$(3)
　が得られる。
　各軸の切片は$\,P_x’=C/P_x\,,\,P_y’=C/P_y\,,\,P_z’=C/P_Z\,$
　となる。

　この式を変形すると　$\quad f(x,y)=z=P_x/P_z\,x+P_y/P_z\,y-D/P_z\,\quad\cdots\,$(4)

　また、この平面上の傾き$\,\mathrm{grad}\,f\,$は、一定で
$\quad \mathrm{grad}\,f=\left(\pdr{z}{x},\pdr{z}{y}\right)=(P_x/P_z\,,\,P_y/P_z)\quad\cdots\,$(5)　　となる。

2-1. 最大傾斜

最大傾斜角度は
$\,\theta=\arctan\sqrt{\left(\left|\dfrac{dz}{dx}\right|^2+\left|\dfrac{dz}{dy}\right|^2\right)}\quad\cdots\,$(6)
　であるので、実際の計算は
$\,180/\pi\times\arctan\sqrt{\left(\dfrac{P_x}{P_z}\right)^2+\left(\dfrac{P_y}{P_z}\right)^2}\quad\cdots\,$(7)　となる。

2-2. 傾斜方向

傾斜方向は
$\,aspect=\arctan 2\left(\dfrac{dz}{dy}\,,\,-\dfrac{dz}{dx}\right)\quad\cdots\,$(8)　となっているので、実際の計算は
$\,180/\pi\times\arctan 2\left(\dfrac{P_y}{P_z}\,,\,-\dfrac{P_x}{P_z}\right)\quad\cdots\,$(9)　となる。

なお、$\,\arctan 2(x,y)=y/x\,$であるが、$\,x\ne 0\,$や符号によって値が変わってくるので、注意のこと。

3. 面同士の角度

交差する面の法線ベクトルを$\,(Q_x,Q_y,Q_z)^T\,$とする。法線ベクトル同士の交差する角度が面の公差角度になるので、
$\quad\cos\theta=\dfrac{|P_xQ_x+P_yQ_y+P_zQ_z|}{\sqrt{P_x^2+P_y^2+P_z}\sqrt{Q_x^2+Q_y^2+Q_z^2}}\quad\cdots\,$(3)

4. 参考資料

Esri 社の Web ページを参考にしました。

傾斜角 (Slope) の仕組み

傾斜方向の仕組み

制限3体 program

制限3体

円制限3体をルンゲクッタ法で解く
近接連星系でラグランジュ点から放たれた点の軌跡を積分する
jupyter-notebook にて program したもの。

In?[9]:

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>

In?[?]:

0.

制限3体

1.制限３体問題

制限３体問題とは、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)}\,$を採用する。
時間の単位として、以下の値を採用する
ケプラーの第３法則から$\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.方程式の数値積分

　方程式の積分で差分法によるものは、最初に差分表の作成が必要となる。それで、差分表を作らずに数値積分する方法として、ルンゲクッタ法がある。
　ルンゲクッタ法はどのような常微分方程式にも使用できる。ただ、そのためには、方程式の系を１階の常微分方程式の系に直さなければならない。そのために、次のように連立させて解法する形とする。また変数は$\,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$が導かれる。

楕円の性質

1.楕円とは

「2定点$\,F，F’\,$からの距離の和が一定である点$\,P\,$の軌跡」を楕円という。

2.楕円の方程式

２定点の座標を、$F\,(c,0)\,\,F’\,(-c,0)\,$とし、点$\,P\,(x,y)\,$とする。距離の和が一定なので、

$\quad \sqrt{(x-c)^2+y^2}+\sqrt{(x+c)^2+y^2}=2a\qquad\cdots\,$(1)

この(1)式を展開変形すると

$\quad(x+c)^2+y^2=4a^2\sqrt{(x-c)^2+y^2}+(x-c)^2+y^2$

$\quad a\sqrt{(x-c)^2+y^2}=a^2-cx$

$\quad a^2(x-c)^2+a^2y^2=a^4-2a^2cx+c^2x^2$

$\quad a^2x^2=c^2x^2+a^2y^2=a^4-a^2c^2$

$\quad x^2(a^2-c^2)+a^2y^2=a^2(a^2-c^2)\qquad\cdots\,$(2)

ここで、$a>c>0\,$より　$\sqrt{a^2-c^2}=b\,$とおくと

$\quad x^2b^2+a^2y^2=a^2b^2\quad\longrightarrow\quad\dfrac{x^2}{a^2}+\dfrac{y^2}{b^2}=1\quad\cdots\,$(3)

この(3)式は通常現れる楕円の式である。

3.極座標における楕円

三角形 $FF’P\,$を考える。この三角形に余弦定理を使う。

$\quad \overline{PF’}^2=\overline{FF’}^2+\overline{FP}^2-2\overline{FF’}\,\overline{FP}\,\cos\angle FF’P$

ここで、$\overline{FF’}=2c\,\,,\,\,\overline{PF}=r\,\,,\,\,\overline{PF’}=2a-r\,\,,\,\,\angle FF’P=\pi-\theta\,$を代入すると

$\quad (2a-r)^2=r^2+(2c)^2-2r\cdot2\cos(\pi-\theta)$

$\quad 4a^2-4ar+r^2=r^2+4c^2+4r\cos\theta$

$\quad a^2-ar=c^2+rc\cos\theta$

$\quad a^2-c^2=r(a+c\cos\theta)$

$\quad r=\dfrac{a^2-c^2}{a+c \cos\theta}=\dfrac{\dfrac{a^2-c^2}{a}}{1+\dfrac{c}{a}\cos\theta}$

ここで、$\dfrac{c}{a}=\ve\,\,,\,\,\dfrac{a^2-c^2}{a}=\ell\quad$とおくと

$\quad r=\dfrac{\ell}{1+\ve\cos\theta}\qquad\cdots\,$(4)

となり、ケプラーの第1法則での式が現れる。

4.楕円の性質

点$\,P$が$\,(c,0)\,\,or\,\,(-c,0)\,$にあるときは定義により$\quad R_1+R_2=2a\quad$が成り立つ。

点$\,P$が$\,(c,0)\,$にあるときは$\quad 2c=R_2-R_1\quad$となる。

また図のように点$\,P$が$\,(0,b)\,$にあるときは、$\,R_1=R_2=a\,\quad$なので$\quad b^2=a^2-c^2\quad\cdots\,$(5)　となる。

これらを合わせると
$\quad b^2=a^2-c^2$

$\hspace{9mm} =\dfrac{1}{4}\left(R_1^2+2R_1R_2+R_2^2-R_1^2+2R_1R_2-R_2^2\right)$

$\hspace{9mm} =R_1R_2\qquad\cdots\,$(6)

また、式(5)より離心率$\,\ve\,$は次のように定義される。

$\quad\ve=\dfrac{\sqrt{a^2-b^2}}{a}\qquad\cdots\,$(7)

これより焦点の座標は $\,c=\ve a\,\quad\Longrightarrow\quad(\pm \ve a\,,\,0)\,$と表す事が出来る。

5.楕円の面積

面積は楕円の方程式を積分して求めることする。計算は第1象限のみを行い結果を4倍して面積とする。

$\quad S=4\di\int_0^ab\sqrt{1-\dfrac{x^2}{a^2}}\,dx\qquad\cdots\,$(8)

ここで、$\,x=a\cos\theta\,\,\left(0\ge\theta\ge\frac{\pi}{2}\right)\,\,$と置換すると

$\quad\dfrac{dx}{d\theta}=-a\sin\theta\quad$となるので、式(8)は

$\quad S=4\di\int_{\frac{\pi}{2}}^0b\sqrt{1-\cos^2\theta}\,(-a)\sin\theta\,d\theta$

$\hspace{9mm}=4\di\int_0^{\frac{\pi}{2}}ab\sin^2\theta\,d\theta\notag$

$\hspace{9mm}=4\di\int_0^{\frac{\pi}{2}}ab\,\dfrac{1-\cos 2\theta}{2}\,d\theta$

$\hspace{9mm}=\left[2ab\right]_0^{\frac{\pi}{2}}-\left[\sin 2\theta\right]_0^{\frac{\pi}{2}}$

$\hspace{9mm}=\pi ab\qquad\cdots\,$(10)$\qquad$となる。

音速

音速を流体力学的に導出する。中央大学棚橋先生の「流体力学」を参考させてもらいました。

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\,\,\,$２原子分子を考え、式(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

電離層

$\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. 地球電離圏

　地球電離圏は，太陽からの極端紫外線(Extra Ultra Violet: EUV)放射によって部分電離された地球の超高層大気領域である．$60\sim 800$km程度の高度範囲に存在し，国際宇宙ステーションや低高度を周回する人工衛星は，電離圏の中を飛翔している．電子密度および中性大気密度の鉛直構造を次に示す．

　$300$km付近に現れる電離圏のピーク高度では電子密度が $10^{11}$から$10^{12}$m$^{-3}$程度にまで達するが，中性大気の密度は同高度で$10^{15}$m$^{-3}$程度であるため，地球電離圏はもっとも濃い部分ですら$0.1\%$程度が電離しているに過ぎない弱電離プラズマであると言える．このため，地球電離圏では，電離大気が背景に存在する中性大気と運動量，熱量の交換を行うことにより，完全電離プラズマ中では見られないような物理過程が生じる．

また，地球電離圏は，特に高緯度領域において，その上部に広がる地球磁気圏と磁力線を介して結合しており，磁気圏からの電場の投影，および荷電粒子の降下によって大きく影響を受けることになる．

2. 酸素の電離

酸素分子は以下の工程でイオン化する．

$(1)\,\,\mathrm{O}_2\,\to\,2\mathrm{O}\hspace{18mm}5.12\,\,$[eV]
$(2)\,\,\mathrm{O}\,\to\,\mathrm{O}^++\,e\hspace{10mm}13.62\,\,$[eV]
$(3)\,\,\mathrm{O}_2\,\to\,\mathrm{O}_2^++\,e\hspace{8mm}12.07\,\,$[eV]

ここで、$\quad h\nu=h\dfrac{c}{\lambda}=E\,$
から光反応の波長を計算する．

$\,h=6.626\times 10^{-34}\,\,$J･s$\qquad c=2.998\times 10^8\,\,$m/s$\qquad e=1.602\times 10^{-19}\,\,$C$\qquad$より
$\,\lambda\,$[nm]$\,\times\,$E [eV]$\,=1240\quad$となる．よって各反応の光の波長は

$(1)\,\,\mathrm{O}_2\,\to\,2\mathrm{O}\hspace{18mm}5.12\,\,$[eV]$\qquad 242\,$nm
$(2)\,\,\mathrm{O}\,\to\,\mathrm{O}^++\,e\hspace{10mm}13.62\,\,$[eV]$\qquad 91\,$nm
$(3)\,\,\mathrm{O}_2\,\to\,\mathrm{O}_2^++\,e\hspace{8mm}12.07\,\,$[eV]$\qquad 103\,$nm

　この図は大気圏の高度と酸素の反応と光の波長を示したものです。ここでは触れないが、オゾンの生成・消長が高度によって変わる関係を表している。

　右端は高度による温度変化を表したものである。
　
　

　
　
　
　酸素プラズマの温度と各粒子の組成の関係を表したものである。

　左端は$\,5000\,$K からであるが、この温度では酸素原子状態が大半である。$8000\,$K 位から、$\mathrm{O}^+\,$が増加してくる。
　
　$20000\,$K 位から、$\mathrm{O}^{++}\,$イオンが増加してくる。

　熱圏での温度は$\,1500\,$K 程度なので、$\mathrm{O}\,$と$\,\mathrm{O}^+\,$が卓越している状態。
　
　
　
　
　
　

　酸素分子とオゾンの紫外吸収スペクトルを示す。

　前述したように、酸素分子は紫外域で励起され、解離＆電離などの反応で、オゾンを生成する。
　生成されたオゾンは紫外域全体に強い吸収を示し、熱圏の温度を高めている。

3. プラズマ振動

　電子とイオンが混ざり合った一様な状態があるとしよう。そこである電子のクラスター（塊）が少しだけ空間的変位すると、それらは他の電子、イオンにより反対方向の力を受けるであろう。すなわちバネの運動と同じように振動することが期待される。これはいろいろなプラズマで観測される典型的な振動であり、プラズマ振動と呼ばれる。

　平衡状態ではイオンと電子が一様に分布している。ところが何らかの原因で、それらの成分の分布に乱れが生じると、ある場所の電子の分布密度が変化する。このときイオンは電子に比べて重いので動かないものとする。

　さて電子が集合すると、電子間にはクーロンの斥力が作用して、もとの平衡点に戻ろうとする。しかし、電子には慣性があるため、その平衡点から行き過ぎて別の場所の密度分布が大きくなってしまう。こうしてプラズマ内の分子は振動する。

　問題を簡単にするため、一次元で考える。左図の$\,x\,$点にあった電子が、その平衡点から$\,s\,$だけ変位し、また同じ時刻に$\,x+\mathit{\Delta}t\,$の点にあった電子は$\,s+\mathit{\Delta}s\,$だけ移動したとする。するとはじめ$\,\mathit{\Delta}x\,$の間にあった電子群は、$\mathit{\Delta}x+\mathit{\Delta}s\,$の間に移動することになる。電子の数は保存するから以下の式が成り立つ。

$\qquad n_0\mathit{\Delta}x=n(\mathit{\Delta}x+\mathit{\Delta}s\,)\qquad\cdots\,$(1)

ここで$\,n_0\,$は平衡の状態の電子の平均数密度、$n\,$は移動後の平均数密度である。(1)から

$\qquad n=\dfrac{n_0}{1+\frac{\mathit{\Delta}s}{\mathit{\Delta}x}}\simeq n_0\left(1-\dfrac{ds}{dx}\right)\qquad\cdots\,$(2)

　空間内には、電子の他に動かない正イオンが平均数密度$\,n_0\,$で一様に分布しているから、空間内の電荷密度のズレ$\,\rho\,$は

$\quad\rho=-en+en_0=n_0e\dfrac{ds}{dx}\qquad\cdots\,$(3)

　ガウスの法則は以下のようになっている。
$\qquad \mathrm{div}\,\bm{E}(\bm{x})=\dfrac{\rho(\bm{x})}{\ve}\qquad\cdots\,$(4)$\qquad$なので

$\qquad \dfrac{dE}{dx}=-\dfrac{n_0e}{\ve}\dfrac{ds}{dx}\qquad\cdots\,$(5)

　これを$\,x\,$について積分すると、場所$\,x\,$に生じる電場は

$\qquad E(x)=\dfrac{n_0e}{\ve}s-C\qquad\cdots\,$(5)
　ここで、電子の平衡点からのズレ$\,s\,$が$\,0\,$のとき、電場は$\,0\,$であるから$\,C=0\,$である。
このとき電子の運動方程式は

$\qquad m\dfrac{d^2s}{dt^2}=-eE=-\dfrac{n_0e^2}{\ve}s\qquad\cdots\,$(6)

　つまり電子は角振動数$\,\,\omega_p=\sqrt{\dfrac{n_0e^2}{\ve m}}\,\,$で振動する。この$\,\omega_p\,$をプラズマの特性振動数という。
　電子の定数を入れると

$\qquad \omega_p=\sqrt{\dfrac{n_0e^2}{\ve m}}=5.64\times10\sqrt{n_0[\mathrm{m}^{-3}]}\quad[\mathrm{s}^{-1}]\qquad\cdots\,$(7)

　電離圏の電子密度は$\,\,10^{11}\sim 10^{12}\,[\mathrm{m}^{-3}]\,\,$なので、プラズマ振動数は数Mhz～数十Mhz となり、短波(HF)が反射されることで、地球の裏側への通信が可能となる。
　
　
　
　

◆ 参考資料

宇宙天気ミニ講座・電離圏編：https://sw-forum.nict.go.jp/forum/2019/pdf/kouza_3.pdf
古典的プラズマ振動について：加藤雄介研究室　https://webpark1378.sakura.ne.jp
電磁気学演習：砂川重信　岩波書店　物理テキストシリーズ-5　2020年第35刷

ガス粒子の衝突

$\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.ガス粒子の衝突

　宇宙での流体を考える際、ガスの力学的な時間に対して電子・分子の衝突時間が十分に短い場合、流体近似が成り立つ。
　これは速度空間が十分に緩和して、連続体として見なせる状態である。気体を１種類の分子からなっており、半径$\,r\,$の剛体球で、ある分子１個が一定の速さ$\,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\,$まず、全部の分子が同じ速度で運動している場合を考える。
１つの分子の速度(速さではなく)を$\,\bm{v}\,$とし、もう一つの分子の速度を$\,\bm{v}’\,$とする。２つの分子の相対速度を$\,\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)、日本評論社
京大「エネルギー理工学設計演習・実験２(別冊)」
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

In [16]:

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)

In [ ]:

角柱周りの流れ

Body-stream-4

$\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}}$

角柱周りの流れ(SMAC法)

流路にブロックが置かれた条件で、流れの時間発展を計算します。

非圧縮性流体の解析で使用する方程式は以下の2式である。($\,\bm{f_B}\,$は重力や電磁力などの体積力である)

$\left\{\begin{array}{l}\nabla\cdot\bm{V}=0\hspace{20mm}\cdots\,(1)\\
\pdr{\bm{V}}{t}+(\bm{V}\cdot\nabla)\bm{V}=-\nabla p+\dfrac{1}{\mathrm{Re}}\nabla^2\bm{V}+\bm{f_B}\qquad\cdots\,(2)\end{array}\right.$

出口条件の$\,U_c\,$は出口断面の平均速度とし、入口速度$\,U_{in}\,$に等しい。また流れ方向の圧力勾配は$\,0\,$とする。この条件を前進オイラー法と1次風上差分により、離散化したのが次の式である。

$\left\{\begin{array}{l}u_{N_i\,,\,j}^{n+1}=u_{N_i\,,\,j}^n-U_c\dfrac{u_{N_i\,,\,j}^n-u_{N_i-1\,,\,j}^n}{\varDelta x}\varDelta t\hspace{20mm}\cdots\,(3)\\
v_{N_i+1\,,\,j}^{n+1}=v_{N_i+1\,,\,j}^n-U_c\dfrac{v_{N_i+1\,,\,j}^n-v_{N_i-1\,,\,j}^n}{\varDelta x}\varDelta t\hspace{20mm}\cdots\,(4)\\
p_{N_i\,,\,j}^{\,\prime}=p_{N_i\,,\,j}^{\,\prime}\hspace{20mm}\cdots\,(5)\end{array}\right.$

圧力に関し、流出境界で法線方向の勾配$\,0\,$とするため、圧力補正に関する境界条件が計算領域全周で勾配$\,0\,$のノイマン型となる。このため、計算領域全体での流入流量と流出流量が一致していないと、圧力補正に関するポアソン方程式

$\nabla^2p^{\,\prime}=\dfrac{1}{\mathrm{Re}\varDelta t}\nabla\cdot\bm{V}^*\hspace{20mm}\cdots\,(6)$

が破綻する。そこで断面平均速度$\,\bar{u}_i\,$を次の式で定義し、

$\bar{u}_i=\di\sum_{j=1}^{N_i}u_{i\,,\,j}\,\,/\,N_j\hspace{20mm}$出口速度を次のように修正する。$\hspace{20mm}u_{N_i\,,\,j}=\,u_{N_i\,,\,j}\,\dfrac{U_{in}}{\bar{u}_{N_i}}\hspace{20mm}\cdots\,(7)$

Re 数を $40$ と低めで設定し、時間刻み係数$\,C_n\,$を$\,0.1\,$とし、圧力補正のポアソン方程式は SOR法にて解法した。この緩和係数を 1.85 とした。静止状態から流入して時間 120 までの時間発展を計算した。ポアソン方程式の解法に、ほとんどの計算を費やすことになった。Cでコンパイルしたものは、20秒で計算終了したが、Python では 21時間必要であった。

◆ 参考文献

数値流体解析の基礎：肖鋒、長崎孝夫(2020)、コロナ社
流体力学第２版：杉山弘、遠藤剛、新井隆景(2014)、森北出版

In [5]:

import numpy as np
from matplotlib import pyplot as plt
import japanize_matplotlib
import warnings
import time
# 警告非表示
warnings.simplefilter('ignore')

# -------  条件設定
ni = 250        # 分割数    X
nj = 50         #           Y
Lx = 25         # 計算領域のサイズ X
Ly = 5          #                  Y
Lb = 1          # ブロックの幅
Lx1 = 4         # 入口からブロックまでの距離
scheme =1       # 中心差分  cn=0.1   2: 風上 cn=0.5   3: ハイブリッド cn=0.5
cn = 0.1        # クーラン数 CFL 条件 < 1 scheme=1 の場合
re = 100.       # レイノルズ数
uin = 1.        # 流入速度
endtime = 120   # 処理終了時間
visc = 1. / re  # 粘度
flg = 1         # 処理フラッグ

dx = Lx / float(ni)       # セルの大きさ X
dy = Ly / float(nj)       #              Y
dxi = 1. / dx
dyi = 1. / dy
dxi2 = dxi * dxi
dyi2 = dyi * dyi
dt = cn / (uin * dxi + 2 * visc * (dxi2 + dyi2))   # 時間刻み

ni1 = int(Lx1 * dxi + 1e-6)               # ブロックまでのXセル数
ni2 = int((Lx1 + Lb) * dxi + 1e-6)          # ブロックの後端までのXセル数
nj1 = int((Ly / 2 - Lb / 2) * dyi + 1e-6)       # ブロック部分の通り道の上部セル数
nj2 = int((Ly / 2 + Lb / 2) * dyi + 1e-6)       #                       下部セル数

alpha = 1.85    # SOR法の緩和係数
eps = 1.e-6     # 繰り返し演算打ち切り誤差
vdx = visc * dxi
vdy = visc * dyi

# -------- 配列
u   = np.zeros((ni+2,nj+2))   # 流速    u
v   = np.zeros((ni+2,nj+2))   #         v
un  = np.zeros((ni+2,nj+2))   # 新流速  u
vn  = np.zeros((ni+2,nj+2))   #         v
p   = np.zeros((ni+2,nj+2))   # 圧力
pc  = np.zeros((ni+2,nj+2))   # 圧力補正
fe  = np.zeros((ni+2,nj+2))   # 検査体積右界面 流速
fn  = np.zeros((ni+2,nj+2))   #         上界面 流速
a_e  = np.zeros((ni+2,nj+2))  # ポアソン方程式係数東
a_w  = np.zeros((ni+2,nj+2))  #                   西
a_n  = np.zeros((ni+2,nj+2))  #                   北
a_s  = np.zeros((ni+2,nj+2))  #                   南
ap  = np.zeros((ni+2,nj+2))   #                   圧力
b   = np.zeros((ni+2,nj+2))   #                   div(u)
iu = np.zeros((ni+2,nj+2),dtype=np.int32)   # 処理フラッグ  u
iv = np.zeros((ni+2,nj+2),dtype=np.int32)   #               v
ip = np.zeros((ni+2,nj+2),dtype=np.int32)   #               p

# グラフ条件設定
nxh = ni // 5
nyh = nj // 5
uc = np.zeros((nyh,nxh))     # ベクトル表示用 u
vc = np.zeros((nyh,nxh))     # ベクトル表示用 v
x, y = np.meshgrid(np.linspace(0, Lx, nxh), np.linspace(0, Ly, nyh))
pd = np.zeros((nyh,nxh))     # 圧力表示用
xp, yp = np.meshgrid(np.linspace(0, Lx, nxh), np.linspace(0, Ly, nyh))

# -------- 初期設定
#
# 流れ u
for i in range(1,ni):
    for j in range(1,nj+1):
        iu[i][j] = flg  # 全体に flg 設定
for i in range(ni1,ni2+1):
    for j in range(nj1+1,nj2+1):
        iu[i][j] = 0    # ブロック部分は設定なし
# 流れ v
for i in range(1,ni+1):
    for j in range(1,nj):
        iv[i][j] = flg  # 全体に flg 設定
for i in range(ni1+1,ni2+1):
    for j in range(nj1,nj2+1):
        iv[i][j] = 0    # ブロック部分は設定なし
# 圧力 p
for i in range(1,ni+1):
    for j in range(1,nj+1):
        ip[i][j] = flg  # 全体にフラッグ設定
for i in range(ni1+1,ni2+1):
    for j in range(nj1+1,nj2+1):
        ip[i][j] = 0    # ブロック部分は設定なし

# 圧力ポアソン方程式の係数設定
for i in range(1,ni+1):
    for j in range(1,nj+1):
        a_w[i][j] = dxi2  # 西
        a_e[i][j] = dxi2  # 東
        a_s[i][j] = dyi2  # 南
        a_n[i][j] = dyi2  # 北
        ap[i][j]= (dxi2 + dyi2) * 2
        if j==1:  # 北側
            ap[i][j] = ap[i][j] - a_s[i][j]
            a_s[i][j] = 0
        if j==nj: # ブロックの上
            ap[i][j] = ap[i][j] - a_n[i][j]
            a_n[i][j] = 0
        if i==1:  # 入口 東側
            ap[i][j] = ap[i][j] - a_w[i][j]
            a_w[i][j] = 0
        if i==ni: # ブロックの西側
            ap[i][j] = ap[i][j] - a_e[i][j]
            a_e[i][j] = 0

# ブロックに隣接するセルの係数変更
for i in range(ni1+1,ni2+1):
    j = nj1
    ap[i][j] = ap[i][j] - a_n[i][j]
    a_n[i][j] = 0
    j = nj2+1
    ap[i][j] = ap[i][j] - a_s[i][j]
    a_s[i][j] = 0

for j in range(nj1+1,nj2+1):
    i = ni1
    ap[i][j] = ap[i][j] - a_e[i][j]
    a_e[i][j] = 0
    i = ni2+1
    ap[i][j] = ap[i][j] - a_w[i][j]
    a_w[i][j] = 0

# 流入境界の速度設定
for j in range(1,nj+1):
    u[0][j] = uin
    un[0][j] = uin

# 初期速度場
for i in range(1,ni):
    for j in range(1,nj+1):
        if iu[i][j] == flg:
            if i >= ni1 and i<= ni2:
                u[i][j] = uin * Ly / (Ly - Lb)
            else:
                u[i][j] = uin
for j in range(1,nj+1):
    u[ni][j] = uin  #  出口の流速
print("* 条件")
print("ni1 :",ni1,"  ni2: ",ni2,"  nj1: ",nj1,"  nj2: ",nj2)

st_time=time.time()
old_time=st_time

# ----------- 時間前進 計算ループ  -------------------------
n = 0
lp_time = 0

while lp_time <= endtime:  # 制限時間内か

    n = n + 1
    lp_time = lp_time + dt

# -- 境界条件
    for i in range(1,ni+1):
        u[i][0] = u[i][1]       # 下境界 du/dy=0
        u[i][nj+1] = u[i][nj]   # 上境界 du/dy=0

    for j in range(1,nj+1):
        v[0][j]= -v[1][j]       # 左境界 v=0
        p[ni+1][j] = p[ni][j]   # 右境界 dp/dx=0

# -- NS方程式： n+1 の仮速度算出

# un
    for i in range(ni+1):       # 流れ U
        for j in range(nj+1):
            u0 = (u[i][j] + u[i+1][j]) / 2
            v0 = (v[i][j] + v[i+1][j]) / 2

            if scheme == 1:

            # 中心差分
                fe[i][j] = u0 * (u[i+1][j] + u[i][j]) / 2 \
                            - vdx * (u[i+1][j] - u[i][j])
                fn[i][j] = v0 * (u[i][j+1] + u[i][j]) / 2 \
                            - vdy * (u[i][j+1] - u[i][j])

            elif schema == 2:
            # 一次風上差分
                fe[i][j] = u0 * (u[i+1][j] + u[i][j]) / 2 \
                            - (vdx + abs(u0)/2) * (u[i+1][j] - u[i][j])
                fn[i][j] = v0 * (u[i][j+1] + u[i][j]) / 2 \
                            - (vdy + abs(v0)/2) * (u[i][j+1] - u[i][j])

            else:
            # ハイブリッド
                fe[i][j] = u0 * (u[i+1][j] + u[i][j]) / 2 \
                            - max(vdx, abs(u0)/2) * (u[i+1][j] - u[i][j])
                fn[i][j] = v0 * (u[i][j+1] + u[i][j]) / 2 \
                            - max(vdy, abs(v0)/2) * (u[i][j+1] - u[i][j])


# -- 角柱の上下面の処理
    for i in range(ni1,ni2+1):
        fn[i][nj1] =  vdy * 2 * u[i][nj1]
        fn[i][nj2] = -vdy * 2 * u[i][nj2+1]

# -- 次の時刻の流れ un
    for i in range(1,ni+1):
        for j in range(1,nj+1):
            if iu[i][j] == flg:
                dudt = (fe[i-1][j] - fe[i][j]) * dxi \
                            + (fn[i][j-1] - fn[i][j]) * dyi
                un[i][j] = u[i][j] + dudt * dt \
                            + (p[i][j] - p[i+1][j]) * dxi * dt

# vn
    for i in range(ni+1):       # 流れ V
        for j in range(nj):
            u0 = (u[i][j] + u[i][j+1]) / 2
            v0 = (v[i][j] + v[i][j+1]) / 2

            if scheme == 1:

            # 中心差分
                fe[i][j] = u0 * (v[i+1][j] + v[i][j]) / 2 \
                            - vdx * (v[i+1][j] - v[i][j])
                fn[i][j] = v0 * (v[i][j+1] + v[i][j]) / 2 \
                            - vdy * (v[i][j+1] - v[i][j])

            elif schema == 2:
            # 一次風上差分
                fe[i][j] = u0 * (v[i+1][j] + v[i][j]) / 2 \
                            - (vdx + abs(u0)/2) * (v[i+1][j] - v[i][j])
                fn[i][j] = v0 * (v[i][j+1] + v[i][j]) / 2 \
                            - (vdy + abs(v0)/2) * (v[i][j+1] - v[i][j])

            else:
            # ハイブリッド
                fe[i][j] = u0 * (v[i+1][j] + v[i][j]) / 2 \
                            - max(vdx, abs(u0)/2) * (v[i+1][j] - v[i][j])
                fn[i][j] = v0 * (v[i][j+1] + v[i][j]) / 2 \
                            - max(vdy, abs(v0)/2) * (v[i][j+1] - v[i][j])

# -- 角柱の上下面の処理
    for j in range(nj1,nj2+1):
        fe[ni1][j] =  vdx * 2 * v[ni1][j]
        fe[ni2][j] = -vdx * 2 * v[ni2+1][j]

# -- 次の時刻の値 vn
    for i in range(1,ni+1):
        for j in range(1,nj):
            if iv[i][j] == flg:
                dvdt = (fe[i-1][j] - fe[i][j]) * dxi \
                            + (fn[i][j-1] - fn[i][j]) * dyi
                vn[i][j] = v[i][j] + dvdt * dt \
                            + (p[i][j] - p[i][j+1]) * dyi * dt

# -- 流出境界の速度 u
    umean = 0
    for j in range(1,nj+1):
        un[ni][j] = u[ni][j] - uin * (u[ni][j] - u[ni-1][j]) * dxi * dt
        umean += un[ni][j]

    umean = umean / nj
    for j in range(1,nj+1):
        if umean != 0:
            un[ni][j] = un[ni][j] / umean * uin
        else:
            un[ni][j] += uin

# -- 速度の発散 div
    for i in range(1,ni+1):
        for j in range(1,nj+1):
            b[i][j] = -((un[i][j] - un[i-1][j]) * dxi \
                            + (vn[i][j] - vn[i][j-1]) * dyi) / dt
            pc[i][j] = 0    # 圧力補正クリア


# -- 圧力補正の計算 ( SOR 法 )

    for iter in range(1,10001):
        err = 0
        for i in range(1,ni+1):
            for j in range(1,nj+1):
                if ip[i][j] == flg:
                    pcn = (a_w[i][j] * pc[i-1][j] + a_e[i][j] * pc[i+1][j] \
                            + a_s[i][j] * pc[i][j-1] + a_n[i][j] * pc[i][j+1] + b[i][j]) / ap[i][j]
                    err = max(err, abs(pcn - pc[i][j]))
                    pc[i][j] = pc[i][j] + alpha * (pcn - pc[i][j])

        if err < eps:
            break


# -- 速度、圧力を修正し、前の値を更新する
    dvmax = 0.
    for i in range(1,ni+1):
        for j in range(1,nj+1):
            vold = v[i][j]
            if iu[i][j] == flg:
                u[i][j] = un[i][j] + (pc[i][j] - pc[i+1][j]) * dxi * dt
            if iv[i][j] == flg:
                v[i][j] = vn[i][j] + (pc[i][j] - pc[i][j+1]) * dyi * dt
            p[i][j] = p[i][j] + pc[i][j]
            dvmax = max(dvmax, abs(v[i][j]-vold))

# -- 流出境界速度
    for j in range(1,nj):
        v[ni+1][j] = v[ni+1][j] - uin * (v[ni+1][j] - v[ni][j]) * dxi *dt

    for j in range(1,nj+1):
        u[ni][j] = un[ni][j]
        
# -- 途中結果の表示
    if n % 200 == 0 or lp_time >= endtime:
#    if n < 5:
        now_time=time.time()
        cal_time=now_time-st_time
        d_time=now_time-old_time
        old_time=now_time
        print("n=",n,"  time=",'{:.2f}'.format(lp_time),"  cal_time(sec)=",'{:.1f}'.format(cal_time),"  sprit(sec)=",'{:.1f}'.format(d_time),"  iter=",iter,"  err=",'{:.3e}'.format(err), "  dvmax=",'{:.3e}'.format(err) )

# -- グラフ表示
    if n % 1000 == 0 or lp_time >= endtime:
#    if n < 5:
        for j in range(nyh):      # V 速度ベクトル
            my = j * 5
            for i in range(nxh):  # U
                mx = i * 5
                uc[j][i]=(u[mx+1][my+1]+u[mx+2][my+1])/2.  # セルの中心
                vc[j][i]=(v[mx+1][my+1]+v[mx+1][my+2])/2.
                pd[j][i]=(p[mx+1][my+1]+p[mx+2][my+2])/2.
         
        plt.figure(figsize=(15,3))
        plt.xlabel("x",fontsize=15)
        plt.ylabel("y",fontsize=15)
        plt.title("流れ",fontsize=20)
        plt.quiver(x,y,uc,vc,color='blue',scale=80,headwidth=4,headaxislength=4,width=0.0015)
        plt.show()
  
        plt.figure(figsize=(14.3,3))
        plt.xlabel("x",fontsize=15)
        plt.ylabel("y",fontsize=15)
        plt.title("流線",fontsize=20)
        plt.streamplot(x, y, uc, vc, density=3, color='k', arrowstyle='-', linewidth=0.6)
        plt.show()

        plt.figure(figsize=(14.3,3.8))
        plt.xlabel("x",fontsize=15)   
        plt.ylabel("y",fontsize=15)
        plt.title("'圧力",fontsize=20)
        cont=plt.contourf(xp, yp, pd, cmap='jet')
#        plt.colorbar(cont, orientation='horizontal')
        plt.colorbar(aspect=40, pad=0.08, shrink=1.0,orientation='horizontal', extend='both')
        plt.show()
 

拡大流れ_SMAC法

$\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}}$

拡大流れ　(SMAC法)

非圧縮性流体の解析で使用する方程式は以下の2式である。($\,\bm{f_B}\,$は重力や電磁力などの体積力である)

ここで(2)式を時刻について差分化して、仮の速度場$\,\bm{V}^{\,*}\,$を求める。

$\quad\dfrac{\bm{V}^{\,*}-\bm{V}}{\Delta t}=-(\bm{V}\cdot\nabla)\bm{V}-\nabla p+\dfrac{1}{\mathrm{Re}}\nabla^2\bm{V}+\bm{f_B}\qquad\cdots\,(3)$

この速度場$\,\bm{V}^{\,*}\,$は連続の式((1)式)を満足していないので、速度と圧力を補正して時間発展を求める。

$\quad \bm{V}^{\,n+1}=\bm{V}^{\,*}+\bm{V}^{\,\prime}\quad,\qquad p^{\,n+1}=p^{\,n}+p^{\,\prime}\qquad\cdots\,(4)$

この(4)式を(1),(2)式に代入して、式(3)を考慮すると以下の式が得られる。

$\left\{\begin{array}{l}\nabla\cdot\bm{V}^{\,\prime}=-\nabla\cdot\bm{V}^{\,*}\hspace{18mm}\cdots\,(5)\\
\dfrac{\bm{V}^{\prime}}{\Delta t}=-\dfrac{1}{\mathrm{Re}}\nabla p^{\,\prime}\hspace{24mm}\cdots\,(6)\end{array}\right.$

(6)式を(5)式に代入すると圧力補正$\,p^{\,\prime}\,$に対するポアソン方程式が得られる。

$\quad\nabla^2 p^{\,\prime}=\dfrac{\mathrm{Re}}{\Delta t}\,\nabla\cdot\bm{V}^{\,*}\qquad\cdots\,(7)$

◆ 手順

１．式(3)を使って$\,\bm{V}^{\,*}\,$を求める
２．求めた$\,\bm{V}^{\,*}\,$から式(7)をSOR法で解いて、$p^{\,\prime}$を求める
３．式(6)を使って、$\bm{V}^{\,\prime}\,$を求める
４．式(4)を用いて、$p^{\,n+1}\,,\,\bm{V}^{\,n+1}\,$を求める

◆ 問題
図のような流路を考える。

計算領域は高さ$\,L_y$、横幅$\,L_x\,$の横長の長方形であり、左境界の上半分の部分($\,x=0\,,\,L_y\,/\,2\ge y\ge L_y\,$)から流体が流入し、右境界から流出する。
メッシュはスタガード格子を用いて解く。

◆ 境界条件
・流入境界
入口の外側($\,i=0\,$のセル)では、$u_{0,j}=U_{in}\,\,,\,v_{0,j}=2V_{in}-v_{1,j}\,\,,\,p^{\,\prime}_{0,j}=p^{\,\prime}_{1,j}\,$となる。
圧力補正値$\,p^{\,\prime}\,$の境界条件は法線方向の勾配$\,0\,$のノイマン条件となる。
・流出境界
出口の外側($\,i=N_i+1\,$のセル)では、各速度成分および物理量が移流によってスムーズに流出すると考え、境界法線方向の勾配を$\,0\,$とすると、$\,u_{N_i+1,j}=u_{N_i,j}\,\,,\,v_{N_i+1,j}=v_{N_i,j}\,$となる。
また、流出境界で断面内の圧力が一様で一定値と見なせるとし、その出口圧力$\,P_{out}\,$を与える。条件は
$\quad p^{\,\prime}_{N_i+1,j}=p^{\,\prime}_{N_i,j}\,\,,\,p_{N_i+1,j}=2P_{out}-p_{N_i,j}\,$となる。

◆ 計算結果
$\quad L_y=1\,\,,\,L_x=5\,\,,\,U_{in}=2\,\,,\,V_{in}=0\,\,,\,P_{out}=0\,$で、格子を$\,65\,\times\,13\,$で分割した結果を示す。
一番目は流れのベクトルプロットである。入口部分に剥離領域が生じているのが解る。次が圧力と流線である。３番目に横方向の$\,4/5\,$の位置の断面の流速分布である。流入流速で正規化してあるので、平均流速は$\,1/2\,$となっている。

◆ 参考文献

数値流体解析の基礎：肖鋒、長崎孝夫(2020)、コロナ社
流体力学第２版：杉山弘、遠藤剛、新井隆景(2014)、森北出版

In [2]:

import numpy as np
from matplotlib import pyplot as plt
#from mpl_toolkits.mplot3d import axes3d
import japanize_matplotlib
# 事前に  > pip install japanize_matplotlib が必要

import warnings
# 警告非表示
#warnings.simplefilter('ignore')

# 条件設定
Lx = 5.
Ly = 1.
nx = 65      # 分割数 X
ny = 13      # 分割数 Y
dx = Lx / float(nx)  # 格子サイズ X
dy = Ly / float(ny)  # 格子サイズ Y
re = 150     # レイノルズ数
vis = 1./re  # 動粘性係数
idx = 1./dx
idy = 1./dy
dx2 = dx * dx
dy2 = dy * dy
idx2 = 1./ dx2
idy2 = 1./ dy2
uin = 2.     # 流入速度
cn = 0.4     # 時間刻み係数
eps= 1.e-6   # SORの収束判定基準
alp = 1.8    # SORの加速係数
nx1 = nx+1
ny1 = ny+1
nx2 = nx+2
ny2 = ny+2
ja = ny // 2    # 入口のサイズ
dt = cn/(uin/dx+2.*vis*(idx2+idy2))   # 時間刻み

# 配列
u  = np.zeros((ny2,nx2))   # 流れ x方向成分 u
v  = np.zeros((ny2,nx2))   # 流れ y方向成分 v
p  = np.zeros((ny2,nx2))   # 圧力
pc = np.zeros((ny2,nx2))   # 圧力補正
un = np.zeros((ny2,nx2))   # 流れ更新分 u
vn = np.zeros((ny2,nx2))   # 流れ更新分 v
fe = np.zeros((ny2,nx2))   # 検査体積の右界面の流速
fn = np.zeros((ny2,nx2))   # 検査体積の上界面の流速

# グラフ条件設定
nxh = nx // 2
nxx = nx // 5 * 2
uc = np.zeros((ny,nxh))     # ベクトル表示用 u
vc = np.zeros((ny,nxh))     # ベクトル表示用 v
x, y = np.meshgrid(np.linspace(0, Lx, nxh), np.linspace(0, Ly, ny))
pd = np.zeros((ny,nx))     # 圧力表示用
xp, yp = np.meshgrid(np.linspace(0, Lx, nx), np.linspace(0, Ly, ny))
uu = np.zeros(ny)

# 圧力ポアソン方程式の係数
a_e = np.zeros((ny2,nx2))   # 方程式の係数 east
a_w = np.zeros((ny2,nx2))   # 方程式の係数 west
a_n = np.zeros((ny2,nx2))   # 方程式の係数 north
a_s = np.zeros((ny2,nx2))   # 方程式の係数 south
a_p = np.zeros((ny2,nx2))   # 方程式の係数 圧力
b  = np.zeros((ny2,nx2))   # 方程式の係数 div V

for j in range(1,ny1):
    for i in range(1,nx1):
        a_w[j][i]=idx2
        a_e[j][i]=idx2
        a_s[j][i]=idy2
        a_n[j][i]=idy2
        a_p[j][i]=2.*(idx2+idy2)
        if j==1:
            a_p[j][i]=a_p[j][i]-a_s[j][i]
            a_s[j][i]=0
        if j==ny:
            a_p[j][i]-a_p[j][i]-a_n[j][i]
            a_n[j][i]=0
        if i==1:
            a_p[j][i]=a_p[j][i]-a_w[j][i]
            a_w[j][i]=0
        if i==nx:
            a_p[j][i]-a_p[j][i]-a_e[j][i]
            a_e[j][i]=0

# 流入口の流れ設定
for j in range(ja+1,ny1):
    u[j][0]=uin     # 流入
    un[j][0]=uin    # 流れ更新分も設定

#
# 計算ループ
#
nt = 120     # 計算ループ回数
time = 0.

for n in range(1,nt+1):

#　境界条件設定

    for i in range(nx1):  # 境界条件 上下
        u[0][i] = -u[1][i]        # 下 u=0
        u[ny1][i] = -u[ny][i]     # 上 u=0
                  
    for j in range(ny1):  # 境界条件 左右
        v[j][0] = -v[j][1]        # 左 v=0
        u[j][nx1] = u[j][nx]      # 右 du/dx=0
        v[j][nx1] = v[j][nx]      #    dv/dx=0
        p[j][nx1] = -p[j][nx]     #    p=0

#　仮の 流速 V* の計算

    for j in range(ny1):   # 流束 u
        for i in range(nx1):
            u1 = (u[j][i]+u[j][i+1])/2.
            v1 = (v[j][i]+v[j][i+1])/2.
            fe[j][i] = u1*(u[j][i+1]+u[j][i])/2.-vis*idx*(u[j][i+1]-u[j][i])
            fn[j][i] = v1*(u[j+1][i]+u[j][i])/2.-vis*idy*(u[j+1][i]-u[j][i])

    for j in range(1,ny1): # u の時間発展
        for i in range(1,nx1):
            dp = (p[j][i]-p[j][i+1])*idx
            dudt=(fe[j][i-1]-fe[j][i])*idx+(fn[j-1][i]-fn[j][i])*idy+dp
            un[j][i]=u[j][i]+dudt*dt

    for j in range(ny):    # 流束 v
        for i in range(nx1):
            u1 = (u[j][i]+u[j+1][i])/2.
            v1 = (v[j][i]+v[j+1][i])/2.
            fe[j][i] = u1*(v[j][i+1]+v[j][i])/2.-vis*idx*(v[j][i+1]-v[j][i])
            fn[j][i] = v1*(v[j+1][i]+v[j][i])/2.-vis*idy*(v[j+1][i]-v[j][i])

    for j in range(1,ny):   # v の時間発展
        for i in range(1,nx1):
            dp = (p[j][i]-p[j+1][i])*idy
            dvdt=(fe[j][i-1]-fe[j][i])*idx+(fn[j-1][i]-fn[j][i])*idy+dp
            vn[j][i]=v[j][i]+dvdt*dt

    for j in range(1,ny1):   # div V の計算
        for i in range(1,nx1):
            b[j][i]=-((un[j][i]-un[j][i-1])*idx+(vn[j][i]-vn[j-1][i])*idy)/dt
            pc[j][i]=0

#　圧力のポアソン方程式の計算 ( SOR法 )

    for itx in range(1000):
        err=0
        for i in range(1,nx1):
            for j in range(1,ny1):
                a1=a_w[j][i]*pc[j][i-1]+a_e[j][i]*pc[j][i+1]
                a2=a_s[j][i]*pc[j-1][i]+a_n[j][i]*pc[j+1][i]
                pcnew=(a1+a2+b[j][i])/a_p[j][i]
                err1=np.abs(pcnew-pc[j][i])
                pc[j][i]=pc[j][i]+alp*(pcnew-pc[j][i])
                if err1 > err:
                    err=err1
                continue
        if err < eps:
               break
        continue

# V と 圧力の更新
    for j in range(1,ny1):
        for i in range(1,nx1):
            u[j][i]=un[j][i]+(pc[j][i]-pc[j][i+1])*idx*dt
            if j != ny:
                v[j][i]=vn[j][i]+(pc[j][i]-pc[j+1][i])*idy*dt
            p[j][i]=p[j][i]+pc[j][i]

    time=time + dt

# 結果の表示
print("n=",n,"   time=",'{:.2f}'.format(time),"   itx=",itx,"  err=",'{:.2e}'.format(err))

for j in range(ny):   # V 速度ベクトル
    for i in range(nxh):  # x 方向の表示は半分に
        m = i*2
        uc[j][i]=(u[j+1][m+1]+u[j+1][m])/2.  # セルの中心
        vc[j][i]=(v[j+1][m+1]+v[j][m+1])/2.

for j in range(ny):   # 圧力
    for i in range(nx):
        pd[j][i]=p[j+1][i+1]

for j in range(ny):
    uu[j]=uc[j][nxx]/uin

# グラフ表示
fig=plt.figure(figsize=(12,4))
ax1 = fig.add_subplot(1, 1, 1)
ax1.set(xlabel='x',ylabel='y')
ax1.set_title('流れ',fontsize=20)
ax1.quiver(x,y,uc,vc,color='blue',scale=80,headwidth=4,headaxislength=4,width=0.0015)
plt.show()

fig=plt.figure(figsize=(12,4))
ax2 = fig.add_subplot(1, 1, 1)
ax2.set(xlabel='x',ylabel='y')
ax2.set_title('圧力',fontsize=20)
ax2.contourf(xp, yp, pd, cmap='cool')
ax2.streamplot(x,y,uc,vc,color='k',density=0.5,linewidth=1)
plt.show()

plt.figure(figsize=(5,5))
plt.title('流れの速度 at Lx = 4',fontsize=20)
plt.ylabel('速度(入口比)',fontsize=16)
plt.xlabel('縦方向',fontsize=16)
plt.plot(uu)
plt.show()

n= 120    time= 1.57    itx= 149   err= 9.44e-07

In [ ]:

円柱まわりの低レイノルズ数流れ

円柱まわりの低レイノルズ数流れ(渦度法)

$\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}}$

円柱まわりの低レイノルズ数流れ　(渦度法)

広い領域に置かれた円柱(無限の高さ)の軸に垂直に一様流が当たっている。
この場合の流れを、非定常の方程式を定常になるまで解くことによって解く。
円柱から遠ざかるほど、流れは一様流に近づき変化は少なくなるので、格子は粗くてよい。

したがって、座標軸の動径方向に$ \,r=e^{\,\xi}\,$という変換を施した。

メッシュは左図のようになる。

返還式は次のようになる。

$\quad\left\{\begin{array}{l}x=e^{\,\xi}\cos\theta\\
y=e^{\,\xi}\sin\theta\qquad\cdots\,(1)\end{array}\right.$

極座標の$\,\theta=0\,$の方向から一様流が当たるとするので、図の右からとなる。

式(1)から流れ関数-渦度法の基礎方程式は

$\quad\left\{\begin{array}{l}\ppdr{\phi}{\xi^2}+\ppdr{\phi}{\theta^2}=-\omega\,e^{\,2\xi}\hspace{75mm}\cdots\,(2)\\
\pdr{\omega}{t}+e^{-2\xi}\left(\pdr{\phi}{\theta}\,\pdr{\omega}{\xi}-\pdr{\phi}{\xi}\,\pdr{\omega}{\theta}\right)=\dfrac{e^{-2\xi}}{\mathrm{Re}}\left(\ppdr{\omega}{\xi^2}+\ppdr{\omega}{\theta^2}\right)\qquad\cdots\,(3)\end{array}\right.$
ここで、$\omega\,,\,\phi\,$は以下のように定義される。

$\quad\omega=\pdr{v}{x}-\pdr{u}{y}\quad\cdots\,(4)\qquad u=\pdr{\phi}{y}\quad\cdots\,(5)\qquad v=-\pdr{\phi}{x}\quad\cdots\,(6)$

この関係から、流れ関数と極座標における速度成分関係式が導かれる。

$\quad\left\{\begin{array}{l}v_r=\dfrac{1}{r}\,\pdr{\phi}{\theta}=e^{-\,\xi}\pdr{\phi}{\theta}\hspace{20mm}\cdots\,(7)\\
v_{\theta}=-\pdr{\phi}{r}=-e^{-\,\xi}\pdr{\phi}{\xi}\hspace{18mm}\cdots\,(8)\end{array}\right.$

◆ 計算
計算では、極座標の $0\,\le\theta\,\lt\,\pi\,$の範囲で計算し、反対側は$\,x\,$軸で折り返す形でコピーする事とする。

・流れ関数の境界条件

$\quad\phi=0\hspace{30mm}$(対称線、円柱上)

$\quad\phi=y=e^{\,\xi}\sin\theta\hspace{25mm}$(遠方)

・渦度の境界条件

$\quad\omega=0\hspace{30mm}$(遠方、対称線)

$\quad\omega=-e^{-2\,\xi}\ppdr{\phi}{\xi^2}\quad\to\quad \omega_0=-2\phi_P\,/\,(\Delta\xi)^2\qquad$(円柱上)

・初期条件

$\quad\left\{\begin{array}{l}\phi=(r-(1/r))\sin\theta=(e^{\,\xi}-e^{-\,\xi})\sin\theta\\
\omega=0\end{array}\right.\hspace{40mm}$(全域)

◆ 参考文献

流体解析の基礎：河村哲也(2014)、朝倉書店

In [12]:

# 円柱まわりの低レイノルズ数流れ(渦度法)
#
import numpy as np
from matplotlib import pyplot as plt
import japanize_matplotlib
import warnings
# 警告非表示
warnings.simplefilter('ignore')

# -------  条件設定
na = 40
ny = 42
nx = na * 2
nxc = na-1
nyc = ny-1
re = 40.        # レイノルズ数
r1 = 1./re
dy = 0.05       # 半径の刻み幅
dt = 0.01       # 時間刻み
td = 1./dt
lmax = 2000     # 繰り返し計算回数
km = 40         # SOR法の最大反復回数
const = 1.0     # SOR法の加速係数
eps = 0.0001    # 繰り返し演算打ち切り誤差

dx = np.pi/float(na)
dxi = 1./dx
dyi = 1./dy
dx2 = dxi*dxi
dy2 = dyi*dyi
fct = .5/(dx2+dy2)

# ---------  配列初期化
psi=np.zeros((ny,na))  # ψ
omg=np.zeros((ny,na))  # ω
tmp=np.zeros((ny,na))  # ωバックアップ
u=np.zeros((ny,nx))    # 流速 u
v=np.zeros((ny,nx))    # 流速 v
omgd=np.zeros((ny,nx)) # ω表示用

# -------- 極座標
th = np.linspace(0, 2*np.pi,nx)
r = np.zeros(ny)
for j in range(nyc):
    r[j+1]=np.exp(j*dy)
TH,R=np.meshgrid(th,r)

#print("th",th)
#print("r=",r)

# --- 初期条件
for i in range(na):  # r 方向の最内側は円筒部分
    psi[0][i]=0
    omg[0][i]=0
    u[0][i]=0
    v[0][i]=0

for j in range(1,ny):
    for i in range(na):
        psi[j][i]=r[j]*np.sin(th[i])   # 遠方  e^ξ sin θ

# ----- 【　メインループ　】------
#

for l in range(1,lmax+1):
    fff=(l-1)/30.
    if fff>= 1.:
        fff=1.

    # STEP 1   境界条件
    for i in range(na):
        omg[1][i]=-2.*psi[2][i]*dy2*fff  # 円筒上
        psi[1][i]=0.                      # φ=0
    
    for i in range(na):
        psi[nyc][i]=r[nyc]*np.sin(th[i]) # 遠方
        omg[nyc][i]=0.
    
    for j in range(ny):                  # 対称線上
        psi[j][0]=0.
        omg[j][0]=0.
        psi[j][nxc]=0.
        omg[j][nxc]=0.
    
    # STEP2 ポアソン方程式を解いて ψ を求める
    
    for k in range(1,km):
        err=0.
        for j in range(2,nyc):
            for i in range(1,nxc):
                psx2=(psi[j][i+1]+psi[j][i-1])*dx2
                psy2=(psi[j+1][i]+psi[j-1][i])*dy2
                rhs=(psx2+psy2+omg[j][i]*np.exp(2.*j*dy))*fct
                err=err+(rhs-psi[j][i])*(rhs-psi[j][i])
                psi[j][i]=psi[j][i]*(1.-const)+rhs*const
        if err < 0.00001:
                break;

    # STEP3 ω の時間発展を求める
    for j in range(2,nyc):
        for i in range(1,nxc):
            tmp[j][i]=omg[j][i]
            
            omx2=(omg[j][i+1]-2.*omg[j][i]+omg[j][i-1])*dx2
            omy2=(omg[j+1][i]-2.*omg[j][i]+omg[j-1][i])*dy2
            psx=(psi[j][i+1]-psi[j][i-1])*.5*dxi
            psy=(psi[j+1][i]-psi[j-1][i])*.5*dyi
            omx=(omg[j][i+1]-omg[j][i-1])*.5*dxi
            omy=(omg[j+1][i]-omg[j-1][i])*.5*dyi
            rhs=((omx2+omy2)*r1+psx*omy-psy*omx)*np.exp(-2.*j*dy)
            omg[j][i]=omg[j][i]+dt*rhs
    
    # STEP4 誤差のチェック
    err1=0.
    for j in range(2,nyc):
        for i in range(1,nxc):
            bb = np.abs(omg[j][i]-tmp[j][i])
            if bb>=err1:
                err1=bb

    if l>10 and err1<=eps:
        break;

# 結果の表示
print( "loop=",l,"(time=",l*dt,")    err1=",'{:.2e}'.format(err1))
for j in range(2,nyc):
    for i in range(1,nxc):
#        vth=-(psi[j][i+1]-psi[j][i-1])*.5*dxi/r[j]
        vth=-(psi[j][i+1]-psi[j][i-1])*.5*dxi/r[j]+np.pi
        vrs=(psi[j+1][i]-psi[j-1][i])*.5*dyi/r[j]
        vth2=2.*np.pi-vth
        u[j][i]=vrs*np.cos(vth)
        u[j][nx-i]=vrs*np.cos(vth2)
        v[j][i]=vrs*np.sin(vth)
        v[j][nx-i]=vrs*np.sin(vth2)
        omgd[j][i]=omg[j][i]
        omgd[j][nx-i]=omg[j][i]

fig=plt.figure(figsize=(10,10))
ax1 = plt.subplot(111, polar=True) 
ax1.set_title('流れ',fontsize=25)
ax1.quiver(TH,R,u,v,color='blue') 
ax1.axis('off')
plt.show() 
fig=plt.figure(figsize=(10,10))
ax2 = plt.subplot(111, polar=True) 
ax2.set_title('ω(渦度)',fontsize=26)
#ax2.contourf(TH,R,omgd,cmap='autumn_r')
ax2.contourf(TH,R,omgd,cmap='PuBu_r')
ax2.axis('off')
plt.show() 

loop= 2000 (time= 20.0 )    err1= 2.67e-04

In [ ]:

1.クーロンの法則

2.ガウスの法則

(1) 積分形

(2) 微分形

3.静電ポテンシャル

(1) ストークスの定理

(2) 静電ポテンシャル

4.電荷のテスト

参考資料

グラジエントのテスト

前進差分

中心進差分

後退差分

テスト

1. 面の式

2. 面の最大傾斜

2-1. 最大傾斜

2-2. 傾斜方向

3. 面同士の角度

4. 参考資料

制限3体

1.制限３体問題

1.1.回転座標系

1.2.計算モデル

1.3.運動方程式

1.4.無次元化

1.5.方程式の数値積分

2.方程式の解法

2.1 Euler 法

2.2 Runge-Kutta法

2.3 Runge-Kutta-Gill法

2.4 Runge-Kutta-Gill法の実際の処理

3.ラグランジュ点

3.1 ラグランジュ点とは

3.2 第1ラグランジュ点 L\(_1\)

3.3 L\(_1\,\)の計算

4.運動軌跡

1.準備

1.1.極座標における速度、加速度

1.2.運動方程式

2.ケプラーの第2法則（面積速度一定の法則）

2.1 導出

2.2 角運動量保存

2.3 面積速度

3.ケプラーの第1法則（楕円軌道の法則）

3.1 導出

3.1.1 step-1

3.1.2 step-2

3.ケプラーの第3法則

1.楕円とは

2.楕円の方程式

3.極座標における楕円

4.楕円の性質

5.楕円の面積

1．音波

1. 地球電離圏

2. 酸素の電離

3. プラズマ振動

1.ガス粒子の衝突

2.宇宙流体での衝突