投稿者「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}}\)



傾きの計算





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



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

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

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

年月日
  • 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
  • マスク越し 爽やか匂う秋の風

    面の角度


    \(\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. 面の式


     点 \(A, B, C\) がある。この3点で構成される面の式を求める。

    ベクトル \(\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体

    \(\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.無次元化

    1. 長さの単位として、両星間の距離\(\,\bm{a}\,\)を採用する。
    2. 質量の単位として、両星の質量の和に万有引力定数を掛けた\(\,\bm{G(M_1+M_2)}\,\)を採用する。
    3. 時間の単位として、以下の値を採用する
      ケプラーの第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法の実際の処理

    手順は以下の通り

    1. 出発点の\(\,(x,y)\,\)を\(\,(x_0,y_0)\,\)とする\(\qquad\)初期値として\(\,q_0=0\,\)とする
    2. 定数として\(\,crk0=1-1/\sqrt{2}\quad,\quad crk1=1+1/\sqrt{2}\,\)とする
    3. \(y(x_0)=y_0\,\)をもとに\(\,k_1=h\cdot f(x_0,y_0)\,\)を計算\(\quad k_1=h\cdot f(x_0,y_0)\)
    4. \(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\)
    5. \(k_1\,,\,q_0\,\)から\(\,q_1\,\)を計算する\(\quad q_1=q_0+3(k_1-2q_0)/2-k_1/2\)
    6. \(x_0+h/2\,\)と\(\,y1\,\)から\(\,k_2\,\)を計算する\(\quad k_2=h\cdot f(x_0+h/2,y_1)\)
    7. \(y_1\,,\,k_2\,,\,q_1\,\)から\(\,y_2\,\)を計算する\(\quad y_2=y_1+\,crk0\,(k_2-q_1)\)
    8. \(q_1\,,\,k_2\,\)から\(\,q_2\,\)を計算する\(\quad q_2=q_1+3\,crk0\,(k_2-q_1)-\,crk0\,k_2\)
    9. \(x_0+h/2\,\)と\(\,y2\,\)から\(\,k_3\,\)を計算する\(\quad k_3=h\cdot f(x_0+h/2,y_2)\)
    10. \(y_2\,,\,k_3\,,\,q_2\,\)から\(\,y_3\,\)を計算する\(\quad y_3=y_2+\,crk1\,(k_3-q_2)\)
    11. \(q_2\,,\,k_3\,q_3\,(\)前回の\(),\,k2\,\)から\(\,q_3\,\)を計算する\(\quad q_3=q_2+3\,crk1\,(k_3-q_3)-crk1\,k_2\)
    12. \(x_0+h\,\)と\(\,y3\,\)から\(\,k_4\,\)を計算する\(\quad k_4=h\cdot f(x_0+h,y_3)\)
    13. \(y_3\,,\,k_4\,,\,q_3\,\)から\(\,y_4\,\)を計算する\(\quad y_4=y_3+(k_4-2q_3)/6\)
    14. \(q_3\,,\,k_4\,\)から\(\,q_4\,\)を計算する\(\quad q_4=q_3+(k_4-2q_3)/2-k_4/2\quad\)次回の\(\,q_0\,\)になる
    15. (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.楕円とは

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

    2.楕円の方程式


    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\)となる。



    音速

    \(\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

    電離層

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

     宇宙での流体を考える際、ガスの力学的な時間に対して電子・分子の衝突時間が十分に短い場合、流体近似が成り立つ。
     これは速度空間が十分に緩和して、連続体として見なせる状態である。気体を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\,\,\)どうもプラズマは、デバイ長やそれ以外の要素が大きいようで、単純な計算ではできないみたい。この検討はここまでとする。


    参考文献








    ロッシュポテンシャル



    ロッシュポテンシャル


    \(\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}}\)



    角柱周りの流れ


    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法


    拡大流れ_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}}\)

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



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



    \(\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}}\)

    室内気流の自然対流(渦度法)



    室内気流の自然対流(渦度法)