投稿者「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-12-20
  • 山眠る最終バスの発った跡
  • 2025-12-20
  • 山眠るランプの宿の闇暗し
  • 2025-11-15
  • 猿橋の柱状節理紅葉渓
  • 2025-11-15
  • 廃線に彩り添える紅葉渓
  • 2025-10-18
  • ずっしりと腰を労り南京哉
  • 2025-09-20
  • 天高し雑草取る手で腰叩く
  • 2025-08-23
  • 処暑の朝きな粉多めのヨーグルト
  • 2025-06-21
  • 紫陽花や瑞々しさの妖精か
  • 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}}\)

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



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


    \(\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}}\)
    ◆ はじめに
    非圧縮性の流れ方程式は、連続の式とナビエストークス方程式(無次元化した)を連立して解くことになる。
    \(\left\{\begin{array}{l}\nabla\cdot\bm{V}=0\hspace{53mm}\cdots\,(1)\\
    \pdr{\bm{V}}{t}+(\bm{V}\cdot\nabla)\bm{V}=-\nabla p+\dfrac{1}{\mathrm{Re}}\Delta\bm{V}\qquad\cdots\,(2)\end{array}\right.\)

    この方程式は非線形で、かつ速度\(\,\bm{V}\,\)については時間発展形になっているが、圧力\(\,p\,\)についてはそうなっていない。そこで、速度\(\,\bm{V}\,\)を時間発展的に求める場合、各時間ステップでの連続の式(式(1))を満足するように圧力\(\,p\,\)を求める必要がある。
    非圧縮性ナビエストークス方程式を数値的に解く方法はいくつかあるが、ここでは圧力を消去する渦度法をここで説明する。

    ◆ 導出
    式(1)(2)を2次元で表すと
    \(\left\{\begin{array}{l}\pdr{u}{x}+\pdr{v}{y}=0\hspace{40mm}\cdots\,(3)\\
    \pdr{u}{t}+u\pdr{u}{x}+v\pdr{u}{y}=-\pdr{p}{x}+\dfrac{1}{\mathrm{Re}}\left(\ppdr{u}{x^2}+\ppdr{u}{y^2}\right)\qquad\cdots\,(4)\\
    \pdr{v}{t}+u\pdr{v}{x}+v\pdr{v}{y}=-\pdr{p}{y}+\dfrac{1}{\mathrm{Re}}\left(\ppdr{v}{x^2}+\ppdr{v}{y^2}\right)\qquad\,\cdots\,(5)\end{array}\right.\)

    まず、式(4)を\(\,y\,\)で偏微分する。左辺は
    \(\quad\pdr{}{t}\left(\pdr{u}{y}\right)+\pdr{u^2}{x\,\partial y}+u\,\ppdr{u}{x\,\partial y}+\pdr{v}{y}\pdr{u}{y}+v\,\ppdr{u}{y^2}\qquad\)となる。変形すると
    \(\quad=\pdr{}{t}\left(\pdr{u}{y}\right)+u\,\ppdr{u}{x\,\partial y}+\pdr{u}{y}\left(\pdr{u}{x}+\pdr{v}{y}\right)+v\,\ppdr{u}{y^2}\qquad\)式(3)を代入すると
    \(\quad=\pdr{}{t}\left(\pdr{u}{y}\right)+u\,\ppdr{u}{x\,\partial y}+v\,\ppdr{u}{y^2}\qquad\cdots\,(6)\quad\)となる。
    右辺は
    \(\quad -\ppdr{p}{x\,\partial y}+\dfrac{1}{\mathrm{Re}}\,\pdr{}{y}\left(\ppdr{u}{x^2}+\ppdr{u}{y^2}\right)\qquad\cdots\,(7)\quad\)となる。

    続けて、式(5)を\(\,x\,\)で偏微分する。左辺は式(4)と同様に変形すると
    \(\quad\pdr{}{t}\left(\pdr{v}{x}\right)+u\,\ppdr{v}{x^2}+v\,\ppdr{v}{x\,\partial y}\qquad\cdots\,(8)\quad\)となる。
    右辺は
    \(\quad -\ppdr{p}{x\,\partial y}+\dfrac{1}{\mathrm{Re}}\,\pdr{}{x}\left(\ppdr{v}{x^2}+\ppdr{v}{y^2}\right)\qquad\cdots\,(9)\quad\)となる。

    ここで、式(8)から式(6)を引いて偏微分の順序を入れ替えると
    \(\quad\pdr{}{t}\left(\pdr{v}{x}-\pdr{u}{y}\right)+u\,\pdr{}{x}\left(\pdr{v}{x}-\pdr{u}{y}\right)+v\,\pdr{}{y}\left(\pdr{v}{x}-\pdr{u}{y}\right)\qquad\cdots\,(10)\quad\)となる。

    式(9)から式(7)を引くと
    \(\quad\dfrac{1}{\mathrm{Re}}\left(\pdr{}{x}\ppdr{v}{x^2}+\pdr{}{x}\ppdr{v}{y^2}-\pdr{}{y}\ppdr{u}{x^2}-\pdr{}{y}\ppdr{u}{y^2}\right)\)
    \(\quad=\dfrac{1}{\mathrm{Re}}\left(\ppdr{}{x^2}\pdr{v}{x}+\ppdr{}{y^2}\pdr{v}{x}-\ppdr{}{x^2}\pdr{u}{y}-\ppdr{}{y^2}\pdr{u}{y}\right)\)
    \(\quad=\dfrac{1}{\mathrm{Re}}\left[\ppdr{}{x^2}\left(\pdr{v}{x}-\pdr{u}{y}\right)+\ppdr{}{y^2}\left(\pdr{v}{x}-\pdr{u}{y}\right)\right]\qquad\cdots\,(11)\quad\)となる。

    ここで渦度\(\,\omega\,\)を流速ベクトル\(\,\bm{V}\,\)の回転とする。
    \(\quad\omega=\nabla\times\bm{V}=\pdr{v}{x}-\pdr{u}{y}\hspace{40mm}\cdots\,(12)\)
    この式(12)を式(10)と式(11)に代入して等号で結ぶと
    \(\quad\pdr{\omega}{t}+u\,\pdr{\omega}{x}+v\,\pdr{\omega}{y}=\dfrac{1}{\mathrm{Re}}\left(\ppdr{\omega}{x^2}+\ppdr{\omega}{y^2}\right)\qquad\cdots\,(13)\)

    流れ関数\(\,\phi(x,y)\,\)を考える。流速ベクトルを流れ関数の回転と定義する。
    \(\quad\bm{V}=\nabla\times\,(0,0,\phi)=\left(\pdr{\phi}{y}\,,\,-\pdr{\phi}{x}\,,\,0\right)\qquad\cdots\,(14)\)
    2次元の場合、\(\bm{V}=(\,u\,,\,v\,,\,0)\,\)の渦度の\(\,z\,\)成分\(\,\omega\,\)は
    \(\quad\omega=\pdr{v}{x}-\pdr{u}{y}=-\left(\ppdr{\phi}{x^2}+\ppdr{\phi}{y^2}\right)\qquad\cdots\,(15)\)

    よって、式(1)(2)は流れ関数\(\,\phi\,\)と渦度\(\,\omega\,\)の連立方程式として、以下の形になる。
    \(\left\{\begin{array}{l}\ppdr{\phi}{x^2}+\ppdr{\phi}{y^2}=-\omega\hspace{60mm}\cdots\,(16)\\
    \pdr{\omega}{t}+\pdr{\phi}{y}\,\pdr{\omega}{x}-\pdr{\phi}{x}\,\pdr{\omega}{y}=\dfrac{1}{\mathrm{Re}}\left(\ppdr{\omega}{x^2}+\ppdr{\omega}{y^2}\right)\qquad\cdots\,(17)\end{array}\right.\)

    ◆ 解き方
    式(16)、(17)を差分近似する。差分については「差分について」を参照のこと。

    \(\quad a_2\phi_{i-1,j}+b_2\phi_{i,j}+c_2\phi_{i+1,j}+a_4\phi_{i,j-1}+b_4\phi{i,j}+c_4\phi_{i,j+1}=-\omega_{i,j}\qquad\cdots\,(18)\)
    \(\quad\dfrac{\omega_{i,j}^{n+1}-\omega_{i,j}}{\Delta t}+(a_3\phi_{i,j-1}+b_3\phi_{i,j}+c_3\phi_{i,j+1})(a_1\omega_{i-1,j}+b_1\omega_{i,j}+c_1\omega_{i+1,j})\)
    \(\hspace{28mm}-(a_1\phi_{i-1,j}+b_1\phi{i,j}+c_1\phi_{i+1,j})(a_3\omega_{i,j-1}+b_3\omega_{i,j}+c_3\omega_{i,j+1})\)
    \(\hspace{15mm}=\dfrac{1}{\mathrm{Re}}(a_2\omega_{i-1,j}+b_2\omega_{i,j}+c_2\omega_{i+1,j}+a_4\omega_{i,j-1}+b_4\omega_{i,j}+c_4\omega_{i,j+1})\qquad\cdots\,(19)\)

    まず、式(18)をSOR法(参照)で解くために、式(18)を\(\,\phi_{i,j}\,\)について解いて

    \(\quad\phi_{i,j}^*=-\dfrac{1}{b_2+b_4}(a_2\phi_{i-1,j}^{(\nu+1)}+c_2\phi_{i+1,j}^{(\nu)}+a_4\phi_{i,j-1}^{(\nu+1)}+c_4\phi_{i,j+1}^{(\nu)}+\omega_{i,j})\qquad\cdots\,(20)\)

    \(\quad\phi_{i,j}^{(\nu+1)}=(1-\alpha)\phi_{i,j}^{(\nu)}+\alpha\phi_{i,j}^*\qquad\cdots\,\,(21)\)

    として、反復解法(SOR法)で求める。ただし、\(\,\nu\,\)は反復回数、\(\,\alpha\,\)は加速係数である。


    参考文献

    • 流体解析の基礎:河村 哲也(2014)、朝倉書店
    • 流体力学第2版:杉山 弘、遠藤 剛、新井 隆景(2014)、森北出版

    流れ関係特性値


    \(\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}}\)
    ◆ レイノルズ数 ( Re )
    レイノルズ数(Reynolds number)は流体力学において慣性力と粘性力との比で定義される無次元量である。流れの中でのこれら2つの力の相対的な重要性を定量している。

    \(\quad\mathrm{Re}=\dfrac{\rho\,U\,L}{\mu}\)

  • \(\rho\):流体の密度 kg/m\(^3\quad\) (例:水:\(1000\) kg/m\(^3\) 空気:\(1.2\) kg/m\(^3\))
  • \(U\):流体の速度 m/s
  • \(L\):代表長さ m
  • \(\mu\):流体の粘度 Pa\(\cdot\)s , kg/m\(\cdot\)s\(\quad\) (例:水:\(1.01\) mPa\(\cdot\)s 空気:\(0.0181\) mPa\(\cdot\)s at \(20\)℃)
  • \(\star\)
    【例】内径 \(2.5\) cm の配管に水が流れた場合、流速 \(U=0.04\) m/s の時は
    \(\quad \mathrm{Re}=\dfrac{1000\times 0.04\times 0.025}{1.01\times 10^{-3}}=1000\quad\)となる。
    \(\qquad\)空気の場合、Re=\(1000\) とするには、流速は \(U=0.6\) m/s 必要となる。
    \(\quad \mathrm{Re}=\dfrac{1.2\times 0.6\times 0.025}{0.018\times 10^{-3}}=1000\quad\)となる。

    ◆ プラントル数 ( Pr )
    プラントル数( Prandtl number)は熱伝導に関する無次元の物性値であり、熱エネルギー拡散の度合いと運動エネルギー拡散の度合いの比です。プラントル数が大きいほど熱エネルギーの拡散に対して運動エネルギーの拡散の度合いが強いことを意味します。プラントル数は流体に固有の物性値です。流動現象において見られる壁面近傍の境界層を例にすると、プラントル数が \(1.0\) の場合、温度境界層と速度境界層の厚みは等しくなります。プラントル数が \(1.0\) より大きくなると、温度境界層に対して速度境界層の方が厚くなり、プラントル数が \(1.0\) より小さい場合はその逆となります。

    \(\quad\mathrm{Pr}=\dfrac{\nu}{\alpha}=\dfrac{\mu\,C_P}{k}\)

  • \(\nu=\mu/\rho\):流体の動粘度 m\(^2\)/s
  • \(\alpha\):熱拡散率(温度伝導率) m\(^2\)/s
  • \(\mu\):流体の粘度 Pa\(\cdot\)s , kg/m\(\cdot\)s\(\quad\) (例:水:\(1.01\) mPa\(\cdot\)s 空気:\(0.0181\) mPa\(\cdot\)s at \(20\)℃)
  • \(C_P\):流体の比熱(定圧比熱) J/kg\(\cdot\)K\(\quad\)(例:水:\(4182\) J/kg\(\cdot\)K 空気:\(1006\) J/kg\(\cdot\)K at \(20\)℃)
  • \(k\):熱伝導率 W/m\(\cdot\)K , J/s\(\cdot\)m\(\cdot\)K\(\quad\)(例:水:\(0.602\) W/m\(\cdot\)K 空気:\(0.0257\) W/m\(\cdot\)K at \(20\)℃)
  • \(\star\)
    水(\(300\)K \(1\)atm):\(5.466\) 空気(\(300\)K \(1\)atm):\(0.717\) 水銀(\(300\)K \(1\)atm):\(0.025\) 潤滑油(\(300\)K \(1\)atm):\(2040\)

    ◆ グランホフ数 ( Gr )
    グラスホフ数(Grashof Number)は、伝熱現象、物質移動現象に関して、流れ場における粘性力に対する浮力の相対的な影響を示す無次元量である。粘性力と浮力の大きさの比を表した無次元数です。自然対流において流れが層流から乱流に遷移するかどうかを特徴付けます。レイノルズ数が強制対流の流れを特徴付ける無次元数ならば、グラスホフ数は自然対流の流れを特徴付ける無次元数です。グラスホフ数が大きい、つまり浮力の影響が強い場合は、流体はより自由に流れようとするため流動は乱流場となります。一方、グラスホフ数が小さい場合は、流体の粘度による流れの抑制効果が高いため層流場となります。層流と乱流の境界となるグラスホフ数はアプリケーションによって異なりますが、概ね \(10^8\)~\(10^9\) のオーダーです。

    \(\quad\mathrm{Gr}=\dfrac{g\beta\Delta\theta L^3}{\nu^2}=\dfrac{g\rho^2\beta\Delta\theta L^3}{\mu^2}\)

  • \(g\):重力加速度 m/s\(^2\) \(9.807\) m/s\(^2\)
  • \(\beta\):体膨張係数 \(1\)/K (例:水:\(0.2\times 10^{-3}/\)K 空気:1/K)
  • \(\Delta\theta\):温度差 K
  • \(\star\)
    空気で代表温度差:\(100\) K 代表長さ:\(1\) m では \(\mathrm{Gr}=1.466\times 10^{10}\) となる。

    差分について


    \(\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}}\)
    ◆ はじめに
    流れ解析において、微分方程式の数値解を求める手法として、まず『差分法』を挙げることが出来る。 関数 \(f(x)\,\)の\(\,x=x_i\,\)における微分は

    \(\quad \dfrac{df}{dx}=\di\lim_{h\to 0}\dfrac{f(x_i+h)-f(x_i)}{h}\quad\cdots\,(1)\)

    で表される。ここで、極限をとらずに\(\,h\,\)が十分小さいとして、”\(\,\sim\,\)” の記号を使って、以下のように表す。

    \(\quad \dfrac{df}{dx}\,\sim\,\dfrac{f(x_i+h)-f(x_i)}{h}\quad\cdots\,(2)\)

    図で言うと、点 \(P\,\)の傾きを\(\,PA\,\)の傾きとするのが『前進差分』である。逆に\(\,BP\,\)とするのが『後退差分』(式(3))である。

    \(\quad \dfrac{df}{dx}\,\sim\,\dfrac{f(x_i)-f(x_i-h)}{h}\quad\cdots\,(3)\)

    点\(\,P\,\)での微分を\(\,BA\,\)の傾きとするのが『中心差分』である。図の差を\(\,h=h_1=h_2\,\)とする。

    \(\quad \dfrac{df}{dx}\,\sim\,\dfrac{f(x_i+h)-f(x_i-h)}{2h}\quad\cdots\,(4)\)

    ◆ 2次微分
    式(4)を\(\,x\,\)で微分すると

    \(\quad \left.\dfrac{d^2f}{dx^2}\right|_{x=x_i}=\dfrac{df^{\prime}}{dx}\,\sim\,\dfrac{f^{\prime}(x_i+h)-f^{\prime}(x_i-h)}{h}\quad\cdots\,(5)\)

    ここで、\(f^{\prime}(x_i+h)\,\)を\(\,f(x_i)\,\)の前進差分(式(1))として展開すると、

    \(\quad \left.\dfrac{df}{dx}\right|_{x=x_i}\,\sim\,\dfrac{f(x_i+h)-f(x_i)}{h}\quad\cdots\,(6)\)

    \(f^{\prime}(x_i-h)\,\)を\(\,f(x_i)\,\)の後退差分(式(2))として展開すると、

    \(\quad \left.\dfrac{df}{dx}\right|_{x=x_i}\,\sim\,\dfrac{f(x_i)-f(x_i-h)}{h}\quad\cdots\,(7)\)

    式(6),(7)を式(5)に代入すると

    \(\quad \left.\dfrac{d^2f}{dx^2}\right|_{x=x_i}\,\sim\,\dfrac{f(x_i+h)-2f(x_i)+f(x_i-h)}{h^2}\quad\cdots\,(8)\)

    ◆ 不均等格子での差分
     左図のような\(\,h_1\ne h_2\,\)の場合、不均等格子となる。三点の線形結合で微分を差分で表すことを考える。\(\,f(x_i)\,\)を\(\,f_i\,\)と表すこととする。

    \(\,a\,f_{i-1}+b\,f_i+c\,f_{i+1}\quad\cdots\,(10)\)

    テーラー展開は次のようになる。

    \(\quad g(e)=g(d)+g^{\prime}(d)(e-d)+\dfrac{g^{\prime\prime}(d)}{2!}(e-d)^2+\cdots\)

    \(\,a\,f_{i-1}\,\)をテーラー展開すると

    \(\quad af(x_i-h_1)=af(x_i)+af^{\prime}(x_i)(-h_1)+a\dfrac{f^{\prime\prime}(x_i)}{2!}(-h_1)^2+\cdots\quad (11)\)

    \(\,c\,f_{i+1}\,\)をテーラー展開すると

    \(\quad cf(x_i+-h_2)=cf(x_i)+cf^{\prime}(x_i)(h_2)+c\dfrac{f^{\prime\prime}(x_i)}{2!}(h_2)^2+\cdots\quad (12)\)

    式(11)と(12)を式(10)に代入すると

    \(\quad a\,f_{i-1}+b\,f_i+c\,f_{i+1}=(a+b+c)f_i+(ch_2-ah_1)f_i^{\prime}+(ah_1^2+ch_2^2)\dfrac{f_i^{\prime\prime}}{2}+\cdots\,(13)\)

    これより、1階微分\(\,f_i^{\prime}\,\)を\(\,f_{i-1}\,\,,\,\,f_i\,\,,\,\,f_{i+1}\,\)の線形結合で表すには

    \(\,a+b+c=0\ \ ,\quad ch_2-ah_1=1\ \ ,\quad 1/2(ah_1^2+ch_2^2)=0\ \)とおいて、\(a\ ,\ b\ ,\ c\ \)を求める。

    \(\,\,a=-\dfrac{h_2}{h_1(h_1+h_2)}\,\,,\quad b=\dfrac{h_2-h_1}{h_1h_2}\,\,,\quad c=\dfrac{h_1}{h_2(h_1+h_2)}\quad\)となる。具体的な近似式は

    \(\quad\left.\dfrac{df}{dx}\right|_{x=x_i}=\dfrac{-h_2^2\,f_{i-1}+(h_2-h_1)(h_2+h_1)\,f_i+h_1^2\,f_{i+1}}{h_1h_2(h_1+h_2)}\quad\cdots\,(14)\)

    同様に、2階微分\(\,f_i^{\prime\prime}\,\)を\(\,f_{i-1}\,\,,\,\,f_i\,\,,\,\,f_{i+1}\,\)の線形結合で表すには

    \(\,a+b+c=0\ \ ,\quad ch_2-ah_1=0\ \ ,\quad 1/2(ah_1^2+ch_2^2)=1\ \)とおいて、\(a\ ,\ b\ ,\ c\ \)を求める。

    \(\,\,a=\dfrac{2}{h_1(h_1+h_2)}\,\,,\quad b=-\dfrac{2}{h_1h_2}\,\,,\quad c=\dfrac{2}{h_2(h_1+h_2)}\quad\)となる。近似式は以下のようになる。

    \(\quad\left.\dfrac{d^2f}{dx^2}\right|_{x=x_i}=\dfrac{2\,h_2\,f_{i-1}-2\,(h_2+h_1)\,f_i+2\,h_1\,f_{i+1}}{h_1h_2(h_1+h_2)}\quad\cdots\,(15)\)

    流れ解析の練習-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}}\)

    流れ解析の練習-1 キャビティ流れ 3D


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

    SOR法による連立一次方程式の解法


    \(\def\bm{\boldsymbol}\)
    \(\def\di{\displaystyle}\)
    \(\def\ve{\varepsilon_0}\)
    \(\def\dd#1#2{\dfrac{\partial #1}{\partial #2}}\)
    \(\def\dda#1#2{\dfrac{\partial^2 #1}{\partial #2}}\)

    ラグランジュ微分


    \(\def\bm{\boldsymbol}\)\(\def\di{\displaystyle}\)\(\def\ve{\varepsilon_0}\)\(\def\dd#1#2{\dfrac{\partial #1}{\partial #2}}\)\(\def\dda#1#2{\dfrac{\partial^2 #1}{\partial #2}}\)

    1. 流体粒子と流体運動の扱い

    流体を考える際、2つの捉え方がある。
    [1] ラグランジュの方法( Lagrange )
    物体(固体)の運動と同様、各流体粒子に着目して、各流体粒子が時間の経過とともに、どのように動くかを追跡していく方法である。
    この方法は、一つの流体粒子の経路や加速度を知る上では便利ではあるが、流れ場が全体的にどのようになっているかを知るには適さない。
    [2] オイラーの方法( Euler )
    この方法は、特定の流体粒子を追跡するのではなく、流れ場全体の様子をそれぞれの時刻に一度に調べる方法である。すなわち、流速や圧力を、座標\(\,\,x\,,\,y\,,\,z\,,\,\)および時間\(\,t\,\)の関数として表す。この方法は、流体の運動を調べる流体力学において、一般に用いられる方法である。

    2. テーラー展開

    まずは、数学の準備をします。
    [1] 1 変数でのテーラー展開
    \(\quad\)関数\(\,f(x)\,\)において、点\(\,a\,\)の周りでテーラー展開すると
    \begin{equation}f(x)=f(a)+f^{\prime}(a)(x-a)+\dots+\dfrac{f^{(n)}(a)}{n!}(x-a)^n+\cdots
    \end{equation}\(\quad\)ここで、\(\,x\,\to\,x+\Delta x\,\quad a\,\to\,x\,\) と置き換えると
    \begin{equation}
    f(x+\Delta x)=f(x)+(\Delta x)\,f^{\prime}(x)+\dfrac{(\Delta x)^2}{2!}f^{\prime\prime}(x)+\cdots+\dfrac{(\Delta x)^n}{n!}f^{(n)}(x)+\cdots
    \end{equation}[2] 多変数のテーラー展開\begin{align}f(x+\Delta x,& y+\Delta y,z+\Delta z)=f(x,y,z)+\left(\Delta x \dd{}{x}+\Delta y\dd{}{y}+\Delta z\dd{}{z}\right)f(x,y,z)\\ &+\dfrac{1}{2!}\left(\Delta x \dd{}{x}+\Delta y\dd{}{y}+\Delta z\dd{}{z}\right)^2f(x,y,z)+\cdots \end{align}\(\quad\)ここで、\(\bm{x}=^t(\,x,\,y,\,z)\,\,\,\,\bm{h}=^t(\Delta x,\Delta y,\Delta z)\,\)とおくと、\begin{equation} f(\bm{x}+\bm{h})=f(\bm{x})+\sum_{n=0}^{\infty}\dfrac{1}{n!}(\nabla\cdot\bm{h})^nf(\bm{x})\qquad\text{と表すことができる}\end{equation}

    3. ラグランジュ微分

    ラグランジュ微分とは、オイラーの方法とラグランジュの方法の両方の流儀を混ぜ合わせたような概念である。時々刻々と移動していく流体の「ある一部分」を追いかけながら、その「一部分」が持つ物理量\(\,A\,\)の時間的な変化を考える。
    時刻\(\,t\,\)に、位置\(\,\bm{x}=(x,y,z)\,\)にある速度\(\,\bm{u}=(u_x,u_y,u_z)\,\)の流体の持つ物理量を\(\,A(x,y,z,t)\,\)と表す。\(\,\Delta t\,\)秒後にはこの流体はおよそ\(\,(\,x+u_x\Delta t,\,y+u_y\Delta t,\,z+u_z\Delta t\,)\,\)付近に移動していることであろう。つまり、\(\,\Delta t\,\)秒後がの移動後の地点での物理量は\(\,A(x+u_x\Delta t,y+u_y\Delta t,z+u_z\Delta t,t+\Delta t)\,\)となる。
    これをテーラー展開の一次で近似すると\begin{align}A(x+u_x\Delta t,& y+u_y\Delta t,z+u_z\Delta t,t+\Delta t)\\&\doteq A(x,y,z,t)+\dd{A}{x}u_x\Delta t+\dd{A}{y}u_y\Delta t+\dd{A}{z}\Delta t+\dd{A}{t}\Delta t\end{align}変化分を計算すると\begin{align}\Delta A&\doteq A(x+u_x\Delta t,y+u_y\Delta t,z+u_z\Delta t,t+\Delta t)-A(x,y,z,t)\\&\doteq\dd{A}{x}u_x\Delta t+\dd{A}{y}u_y\Delta t+\dd{A}{z}u_z\Delta t+\dd{A}{t}\Delta t\\&=\left(\dd{A}{x}u_x+\dd{A}{y}u_y+\dd{A}{z}u_z+\dd{A}{t}\right)\Delta t\qquad\cdots\quad\text{(1)}\end{align}これを時間経過による変化割合で表すと\begin{equation}\dfrac{\Delta A}{\Delta t}\doteq\dd{A}{x}u_x+\dd{A}{y}u_y+\dd{A}{z}u_z+\dd{A}{t}\qquad\cdots\quad\text{(2)}\end{equation}(1)式の近似には\(\,\Delta t\,\)の2乗以上に比例する程度の誤差があり、それを\(\,\Delta t\,\)で割った(2)式には\(\,\Delta t\,\)あるいは\(\,\Delta t\,\)の2乗以上に比例するような誤差が含まれているが、ここで\(\,\Delta t\,\to\,0\,\)の極限を考えれば無視できるようになるので、それを等式を使って次のように表す。\begin{equation}\dfrac{DA}{Dt}=\dd{A}{x}u_x+\dd{A}{y}u_y+\dd{A}{z}u_z+\dd{A}{t}\qquad\cdots\quad\text{(3)}\end{equation}これが、「ラグランジュ微分」「物質微分」などと呼ばれているものである。これは、流体と一緒に流れている人から見た、その地点での物理量\(\,A\,\)の時間的な変化率である。
    物理量\(\,A\,\)を省いて、順序を変えると以下のようになる。
    \begin{equation}\dfrac{D}{Dt}\equiv u_x\dd{}{x}+u_y\dd{}{y}+u_z\dd{}{z}+\dd{}{t}=\bm{u}\cdot\nabla+\dd{}{t}\quad\cdots\quad\text{(4)}\end{equation}流体の加速度をその点での微分とラグランジュ微分で比べたのが、下の図である。

    ベクトル解析



    \(\def\bm{\boldsymbol}\)
    \(\def\di{\displaystyle}\)
    \(\def\ve{\varepsilon_0}\)
    \(\def\dd#1#2{\dfrac{\partial #1}{\partial #2}}\)
    \(\def\dda#1#2{\dfrac{\partial^2 #1}{\partial #2}}\)
    光の速度で、ベクトル解析の公式を使ったので、まとめておきます。

    1.はじめに

    (1) ベクトルの内積
    \(\quad\bm{a}\cdot\bm{b}=a_xb_x+a_yb_y+a_zb_z=\di\sum^3_ia_ib_i\)
    (2) ベクトルの外積
    \(\quad\bm{a}\times\bm{b}=(a_yb_z-a_zb_y,a_zb_x-a_xb_z,a_xb_y-a_yb_x)\)
    \(\quad(\bm{a}\times\bm{b})_i=\di\sum_{j,k=1}^3\epsilon_{i,j,k}\,a_jb_k\)
    \begin{align}\hspace{-20mm}\epsilon_{i,j,k}=1&\,:(i,j,k)=(1,2,3),(2,3,1),(3,1,2)\\
    -1&\,:(i,j,k)=(1,3,2),(2,1,3),(3,2,1)\\
    0&\,: other\end{align}
    3. スカラー三重積
    \(\quad\bm{A}\cdot(\bm{B}\times\bm{C})=\bm{B}\cdot(\bm{C}\times\bm{A})=\bm{C}\cdot(\bm{A}\times\bm{B})\)
    \begin{align}
    \bm{A}\cdot(\bm{B}\times\bm{C})&=A_x(\bm{B}\times\bm{C})_x+A_y(\bm{B}\times\bm{C})_y+A_z(\bm{B}\times\bm{C})_z\\[2pt]
    &=A_x(B_yC_z-B_zC_y)+A_y(B_zC_x-B_xC_z)+A_z(B_xC_y-B_yC_x)\\[2pt]
    &=A_xB_yC_z-A_xB_zC_y+A_yB_zC_x-A_yB_xC_z+A_zB_xC_y-A_zB_yC_x\\[2pt]
    &=B_x(C_yA_z-C_zA_y)+B_y(C_zA_x-C_xA_z)+B_z(C_xA_y-C_yA_x)\\[2pt]
    &=B_x(\bm{C}\times\bm{A})_x+B_y(\bm{C}\times\bm{A})_y+B_z(\bm{C}\times\bm{A})_z\\
    &=\bm{B}\cdot(\bm{C}\times\bm{A})\\[-40pt]
    \end{align}
    \(\qquad\)また、\(\quad\bm{A}\cdot(\bm{B}\times\bm{C})=\left|\begin{array}{ccc}A_x & A_y & A_z\\ B_x & B_y & B_z\\C_x & C_y & C_z\end{array}\right|\qquad\)である。

    4. ベクトル三重積
    \(\quad\bm{A}\times(\bm{B}\times\bm{C})=(\bm{A}\cdot\bm{C})\bm{B}-(\bm{A}\cdot\bm{B})\bm{C}\)
    \begin{align}
    (\bm{A}\times(\bm{B}\times\bm{C}))_x&=A_y(\bm{B}\times\bm{C})_z-A_z(\bm{B}\times\bm{C})_y\\
    &=A_y(B_xC_y-B_yC_x)-A_z(B_zC_x-B_xC_z)\\
    &=(A_yC_y+A_zB_z)B_x-(A_yB_y+A_zB_z)C_x\\
    &=(A_xC_x+A_yC_y+A_zB_z)B_x-(A_xB_x+A_yB_y+A_zB_z)C_x\\
    &=(\bm{A}\cdot\bm{C})B_x-(\bm{A}\cdot\bm{B})C_x\\[6pt]
    &\hspace{-12mm}\text{同様にして}\\[2pt]
    (\bm{A}\times(\bm{B}\times\bm{C}))_y&=(\bm{A}\cdot\bm{C})B_y-(\bm{A}\cdot\bm{B})C_y\\[2pt]
    (\bm{A}\times(\bm{B}\times\bm{C}))_z&=(\bm{A}\cdot\bm{C})B_z-(\bm{A}\cdot\bm{B})C_z\\[8pt]
    &\hspace{-18mm}\therefore\quad \bm{A}\times(\bm{B}\times\bm{C})=(\bm{A}\cdot\bm{C})\bm{B}-(\bm{A}\cdot\bm{B})\bm{C}
    \end{align}
    5. 微分演算子
    \(\quad\nabla=\left(\,\dd{}{x}\,,\,\dd{}{y}\,,\,\dd{}{z}\,\right)\)

    2.勾配、発散、回転

    勾配( grad ):grad\(\,\phi(\bm{r})=\nabla\,\phi(\bm{r})\equiv\left(\,\dd{}{x}f\,,\,\dd{}{y}f\,,\,\dd{}{z}f\,\right)\)

    発散 ( div ):div\(\,\bm{A}(\bm{r})=\nabla\cdot\bm{A}(\bm{r})\equiv\,\dd{}{x}A_x+\dd{}{y}A_y+\dd{}{z}A_z\)

    回転 ( rot ):rot\(\,\bm{A}(\bm{r})=\left(\dd{}{y}A_z-\dd{}{z}A_y\,,\,\dd{}{z}A_x-\dd{}{x}A_z\,,\,\dd{}{x}A_y-\dd{}{y}A_x\,\right)\)
    \(\hspace{8mm}\)または \(\quad(\nabla\times\bm{A})_i=\di\sum_{j,k=1}^3\epsilon_{ijk}(\nabla)_j(\bm{A})_k=\di\sum_{j,k=1}^3\epsilon_{ijk}\,\partial_jA_k\)

    3.計算例-1

    \(\quad\bm{r}=(x,y,z)\quad r=|\bm{r}|\ne 0\quad\)とする。

    (1)\(\,\nabla\,r\)
    \begin{align}
    \nabla\,r&=\nabla(\sqrt{x^2+y^2+z^2})\\[2pt]
    &=\left(\dfrac{x}{\sqrt{x^2+y^2+z^2}}\,,\,\dfrac{y}{\sqrt{x^2+y^2+z^2}}\,,\,\dfrac{z}{\sqrt{x^2+y^2+z^2}}\,\right)\\[2pt]
    &=\dfrac{\bm{r}}{r}
    \end{align}
    (2)\(\,\nabla\left(\dfrac{1}{r}\right)\)
    \begin{align}
    \nabla\left(\dfrac{1}{r}\right)&=\left(\dd{}{x}\dfrac{1}{\sqrt{x^2+y^2+z^2}}\,,\,\dd{}{y}\dfrac{1}{\sqrt{x^2+y^2+z^2}}\,,\,\dd{}{z}\dfrac{1}{\sqrt{x^2+y^2+z^2}}\,\right)\\[3pt]
    &=\left(\,\dfrac{-x}{(x^2+y^2+z^2)^{3/2}}\,,\,\dfrac{-y}{(x^2+y^2+z^2)^{3/2}}\,,\,\dfrac{-z}{(x^2+y^2+z^2)^{3/2}}\,\right)\\[3pt]
    &=-\dfrac{\bm{r}}{r^3}
    \end{align}
    (3)\(\,\nabla\cdot\bm{r}\)
    \begin{equation}
    \nabla\cdot\bm{r}=\dd{}{x}x+\dd{}{y}y+\dd{}{z}z=3
    \end{equation}
    (4)\(\,(\nabla\cdot\nabla)\,r\)
    \begin{align}
    (\nabla\cdot\nabla)r&=\left(\dda{}{x^2}+\dda{}{y^2}+\dda{}{z^2}\right)\sqrt{x^2+y^2+z^2}\\[2pt]
    &=\dd{}{x}\dfrac{x}{\sqrt{x^2+y^2+z^2}}+\dd{}{y}\dfrac{y}{\sqrt{x^2+y^2+z^2}}+\dd{}{z}\dfrac{z}{\sqrt{x^2+y^2+z^2}}\\[2pt]
    &=\dfrac{1}{\sqrt{x^2+y^2+z^2}}-\dfrac{x^2}{(x^2+y^2+z^2)^{3/2}}\\
    & +\dfrac{1}{\sqrt{x^2+y^2+z^2}}-\dfrac{y^2}{(x^2+y^2+z^2)^{3/2}}\\
    & +\dfrac{1}{\sqrt{x^2+y^2+z^2}}-\dfrac{z^2}{(x^2+y^2+z^2)^{3/2}}\\[2pt]
    &=\dfrac{3}{r}-\dfrac{(x^2+y^2+z^2)}{r^3}=\dfrac{2}{r}
    \end{align}
    (5)\(\,\nabla\times\bm{r}\)
    \begin{equation}
    \nabla\times\bm{r}=\left(\dd{}{z}y-\dd{}{y}z\,,\,\dd{}{x}z-\dd{}{z}x\,,\,\dd{}{y}x-\dd{}{x}y\right)=\bm{0}
    \end{equation}
    (6)\(\,\nabla\times(\nabla\phi)\hspace{10mm}\mathrm{rot}\,\,\mathrm{grad}\,\,\phi\)
    \begin{align}
    \nabla\times(\nabla\phi)&=\left(\,\dd{}{y}(\dd{\phi}{z})-\dd{}{z}(\dd{\phi}{y})\,\,,\right.\\
    &\hspace{20mm}\dd{}{z}(\dd{\phi}{x})-\dd{}{x}(\dd{\phi}{z})\,\,,\\
    &\hspace{30mm}\left.\dd{}{x}(\dd{\phi}{y})-\dd{}{y}(\dd{\phi}{x})\,\right)\\
    &=0\hspace{40mm}\because\,\,\dd{}{y}(\dd{\phi}{z})=\dd{}{z}(\dd{\phi}{y})
    \end{align}
    (7)\(\,\nabla\cdot(\nabla\times\bm{A})\hspace{10mm}\mathrm{div}\,\,\mathrm{rot}\,\,\bm{A}\)
    \begin{align}
    \nabla\cdot(\nabla\times\bm{A})&=\dd{}{x}(\dd{A_z}{y}-\dd{A_y}{z})+\dd{}{y}(\dd{A_x}{z}-\dd{A_z}{x})+\dd{}{z}(\dd{A_y}{x}-\dd{A_x}{y})\\[3pt]
    &=\dfrac{\partial^2A_z}{\partial x\partial y}-\dfrac{\partial^2A_y}{\partial x\partial z}+\dfrac{\partial^2A_x}{\partial y\partial z}-\dfrac{\partial^2A_z}{\partial y\partial x}+\dfrac{\partial^2A_y}{\partial z\partial x}-\dfrac{\partial^2A_x}{\partial z\partial y}\\[3pt]
    &=0\hspace{40mm}\because\,\,\,\dfrac{\partial^2A_z}{\partial x\partial y}=\dfrac{\partial^2A_z}{\partial y\partial x}
    \end{align}
    (8)\(\,\nabla\times(\nabla\times\bm{A})\hspace{10mm}\mathrm{rot}\,\,\mathrm{rot}\,\,\bm{A}\)
    \begin{align}
    (\nabla\times(\nabla\times\bm{A}))_x&=\dd{}{y}\left(\dd{A_y}{x}-\dd{A_x}{y}\right)-\dd{}{z}\left(\dd{A_x}{z}-\dd{A_z}{x}\right)\\[3pt]
    &=\dfrac{\partial^2}{\partial y\partial x}A_y-\dda{}{y^2}A_x-\dda{}{z^2}A_x+\dfrac{\partial^2}{\partial z\partial x}A_z\\[3pt]
    &=\dda{}{x^2}A_x+\dfrac{\partial^2}{\partial y\partial x}A_y+\dfrac{\partial^2}{\partial z\partial x}A_z-\left(\dda{}{x^2}+\dda{}{y^2}+\dda{}{z^2}\right)A_x\\[3pt]
    &=\dd{}{x}\left(\dd{A_x}{x}+\dd{A_y}{y}+\dd{A_z}{z}\right)-\Delta A_x\\[3pt]
    &=\dd{}{x}(\nabla\cdot\bm{A})-\Delta A_x\\[3pt]
    &\hspace{-20mm}\text{同様に}\\
    (\nabla\times(\nabla\times\bm{A}))_y&=\dd{}{y}(\nabla\cdot\bm{A})-\Delta A_y\\[3pt]
    (\nabla\times(\nabla\times\bm{A}))_z&=\dd{}{z}(\nabla\cdot\bm{A})-\Delta A_z\\[3pt]
    &\hspace{-20mm}\text{よって}\\[2pt]
    \nabla\times(\nabla\times\bm{A})&=\nabla(\nabla\cdot\bm{A})-\Delta\bm{A}
    \end{align}

    光の速度


    \(\def\bm{\boldsymbol}\)
    \(\def\di{\displaystyle}\)
    \(\def\ve{\varepsilon_0}\)
    \(\def\dd#1#2{\dfrac{\partial #1}{\partial #2}}\)
    \(\def\dda#1#2{\dfrac{\partial^2 #1}{\partial #2}}\)
    オーソドックスな形で、光(電磁波)の真空中の速度を求めてみる。

    [1]マックスウェル方程式

    光も電磁波なので、マックスウェルの方程式に従う。マックスウェル方程式は4つに分けられる。
    (1)電場の発散
    \[\mbox{div }\bm{E}=\nabla\cdot\bm{E}=\dfrac{\rho}{\ve}\tag{1.1}\]
    電場の源が電荷であり、電場が電荷から放射状であることを表す。
    ここで \(\,\rho\,\)は電荷密度で、\(\ve\,\)は真空の誘電率である。
    (2)電場の回転
    \[\mbox{rot }\bm{E}=\nabla\times\bm{E}=-\dd{\bm{B}}{t}\tag{1.2}\]
    磁束密度の変化により電場が生じる Faradayの法則=電磁誘導の法則
    (3)磁場の発散
    \[\mbox{div }\bm{B}=\nabla\cdot\bm{B}=0\quad\tag{1.3}\]
    磁場には源がない
    (4)磁場の回転
    \[\mbox{rot }\bm{B}=\nabla\times\bm{B}=\mu_0(\bm{j}+\ve\dd{\bm{E}}{t})\tag{1.4}\]
    電流及び電場の変化が磁場を生むことを表す。Ampére-Maxwellの法則
    ここで\(\,\mu_0\,\)は真空の透磁率である。

    [2]マックスウェル方程式の展開

    真空中を伝わる電磁波について考えたいので, 電荷密度はいたるところで\(\,0\,\)であるとする. よって電流密度も\(\,0\,\)であるので、次の4つの式が得られる。
    \begin{align}
    \mbox{div }\bm{E}&=0 \tag{2.1}\\
    \mbox{rot }\bm{E}&=-\dd{\bm{B}}{t} \tag{2.2}\\
    \mbox{div }\bm{B}&=0 \tag{2.3}\\
    \mbox{rot }\bm{B}&=\mu_0\ve\dd{\bm{E}}{t} \tag{2.4}
    \end{align}
    場の回転を改めて場とみなして、回転を調べる。式\(\,(2.2)\,\)を回転すると
    \begin{align}
    \nabla\times(\nabla\times\bm{E})&=-\dd{}{t}(\nabla\times\bm{B})\\
    &=-\mu_0\dd{}{t}\left(\bm{j}+\ve\dd{\bm{E}}{t}\right)\\
    &=-\mu_0\dd{\bm{j}}{t}-\mu_0\ve\dda{\bm{E}}{t^2}\tag{2.5}
    \end{align}
    となり、電流密度の他には電場のみの閉じた形となる。ここでは、ベクトル解析のベクトルの回転の回転からラプラシアンを導く、以下の公式を用いる。
    \[\nabla\times(\nabla\times\bm{A})=\nabla(\nabla\cdot\bm{A})-\Delta\bm{A}\qquad\because\,\,\Delta\equiv\nabla\cdot\nabla\tag{2.6}\]
    この公式を\(\,(2.5)\,\)式の右辺に適用する。
    \begin{align}
    \nabla\times(\nabla\times\bm{E})&=\nabla(\nabla\cdot\bm{E})-\Delta\bm{E}\\[3pt]
    &=0-\Delta\bm{E}\qquad\,\because\,\,\nabla\cdot\bm{E}=0\quad(2.1)\,\,\text{より}\\[3pt]
    &=-\Delta\bm{E}\tag{2.7}
    \end{align}
    \(\,(2.5)\,\)式と\(\,(2.6)\,\)をあわせると
    \[\Delta\bm{E}-\mu_0\ve\dda{\bm{E}}{t^2}-\mu_0\dd{\bm{j}}{t}=0\tag{2.8}\]
    ここの電流は\(\,0\,\)なので、第3項は消せる。ここで、\(x\,\)成分で考えると。
    \[\dda{\bm{E}}{x^2}-\mu_0\ve\dda{\bm{E}}{t^2}=0\tag{2.10}\]
    の形の波動方程式が得られる。

    [3]波動方程式の一般解

    この方程式の解は次のような形で表される。例えば\(\,\bm{E}\,\)の振幅を\(\,z\,\)軸にとると、解は
    \[E_z=f(x-vt)+g(x+vt)\tag{3.1}\]
    の形で表せる。ここで、\(f\,\)と\(\,g\,\)は任意の関数であり、\(f(x − vt)\,\,\)は\(\,x\,\)軸の正の向きに速さ\(\,v\,\)で
    進む波動を表し、\(g(x + vt)\,\,\)は\(\,x\,\)軸の負の向きに速さ\(\,v\,\)で進む波動を表す。
    ここで、式\(\,(2.2)\,\)を成分に分解する。
    \begin{align}
    \dd{E_z}{y}-\dd{E_y}{z}&=-\dd{B_x}{t}\\
    \dd{E_x}{z}-\dd{E_z}{x}&=-\dd{B_y}{t}\\
    \dd{E_y}{x}-\dd{E_x}{y}&=-\dd{B_z}{t}
    \end{align}
    ここで、\(E_x=0\,\,,\,\,E_y=0\,\)なので、次のような結果となる
    \begin{align}
    \dd{B_x}{t}&=0\\
    \dd{B_y}{t}&=f^{\prime}(x-vt)+g^{\prime}(x+vt) \tag{3.2}\\
    \dd{B_z}{t}&=0
    \end{align}
    式\(\,(3.2)\,\)を\(\,t\,\)で積分すると
    \[B_y=-\dfrac{1}{v}f(x-vt)+\dfrac{1}{v}g(x+vt)\tag{3.3}\]
    となる。つまり、\(E_z\,\)と\(\,B_y\,\)は、同じ関数\(\,f\,\)と\(\,g\,\)で表され、両者が互いに組み合って、離れることなく、同じ速さ\(\,v\,\)を持つ波動となって\(\,x\,\)軸を伝搬する。

    電場と磁場の進行波

    [4]電磁波の伝搬する速さ

    波動方程式の一般解と波動方程式を再掲する。
    \begin{align}
    &E_z=f(x-vt)+g(x+vt)\tag{3.1}\\
    &\dda{\bm{E}}{x^2}-\mu_0\ve\dda{\bm{E}}{t^2}=0\tag{2.10}
    \end{align}
    この\(\,(3.1)\,\)式を\(\,(2.10)\,\)に代入する。このとき、\(x-vt=p\,\,,\,\,x+vt=q\,\,\)とする。
    \[\dd{E_z}{x}=\dd{f}{x}+\dd{g}{x}=\dd{f}{p}\dd{p}{x}+\dd{g}{q}\dd{q}{x}=\dd{f}{p}+\dd{g}{q}\quad\text{なので}\]
    \[\dda{E_z}{x^2}=\dda{f}{p^2}+\dda{g}{q^2}\qquad\text{となる。次に時間で微分する。}\]
    \[\dd{E_z}{t}=\dd{f}{t}+\dd{g}{t}=\dd{f}{p}\dd{p}{t}+\dd{g}{q}\dd{q}{t}=-v\dd{f}{p}+v\dd{g}{q}\quad\text{なので}\]
    \[\dda{E_z}{t^2}=v^2\dda{f}{p^2}+v^2\dda{g}{q^2}=v^2\dda{E_z}{x^2}\qquad\text{となるので}\quad (2.10)\,\text{より}\]
    \[v^2=\dfrac{1}{\mu_0\ve}\qquad\text{よって}\qquad v=\dfrac{1}{\sqrt{\mu_0\ve}}\tag{4.1}\]
    真空の誘電率\(\,\ve\,\)と透磁率\(\,\mu_0\,\)の値

    \(\qquad\ve=\dfrac{10^7}{4\pi c^2}\quad[\,\mathrm{C}\cdot\mathrm{N}^{-1}\cdot\mathrm{m}^{-2}\,\,]\hspace{20mm}\mu_0=4\pi\times 10^{-7}\quad[\,\,\mathrm{N}\cdot\mathrm{A}^{-2}\,\,]\quad\)を代入すると

    \(\qquad\quad v=\dfrac{1}{\sqrt{\mu_0\ve}}=c\qquad\)となり、電磁波の速度が光速と一致する。

    まあ、誘電率の定義からして、光速を使っているので当然の結果である。







    二重振り子

    1.経緯

    放送大学の「物理と化学のための数学演習」という授業のレポートの課題のひとつに「二重振り子」があったので、忘れないようにここにメモっておきます。

    2.レポート課題

    天井から長さ \(\ell\) の紐で吊るされた質量 \(M\) の錘に、もう一つの同じ長さ \(\ell\) の紐で同じ質量の錘が垂れ下がった、二重振り子の微小振動を考える。この問題は自由度2の連成振動となっている。以下の問いに答えよ。
    \(\textbf{[1]}\) それぞれの紐が鉛直線となす角度を \(\theta_1\,,\,\theta_2\,\)とし、近似式、\(\sin\theta\simeq\theta\,\)を用いて、2つの錘の運動方程式を導き、この連成振動の基準座標と基準振動数を求めよ。紐の質量は小さいとして無視する。
    \(\textbf{[2]}\) 2つの紐を変形しない同じ質量 \(m\) の棒で置き換えるとどうなるか。棒の重心はその真ん中にあるとする。

    ※ 2. の問題は、私のやり方が拙かったせいか、解析解が得られませんでした。そこで、1. だけをここにメモります。

    3.方程式の導出

    最初のオモリの座標を \((x_1,y_1)\,\) , そこから垂れ下がったオモリの座標を \((x_2,y_2)\,\) とする。

    \(x_1=\ell\sin\theta_1\qquad y_1=-\ell\cos\theta_1\)
    \(x_2=x_1+\ell\sin\theta_2=\ell(\sin\theta_1+\sin\theta_2)\)
    \(y_2=y_1-\ell\cos\theta_2=-\ell(\cos\theta_1+\cos\theta_2)\)
    \(\dot{x}_1=\ell\cos\theta_1\dot{\theta}_1\simeq\ell\dot{\theta}_1\)
    \(|\dot{x}_1|^2=\ell^2\dot{\theta}_1^2\)
    \(\dfrac{1}{2}|\dot{x}_1|^2=\dfrac{\ell^2}{2}\dot{\theta}_1^2\)
    \(\dot{x}_2=\ell(\cos\theta_1\dot{\theta}_1+\cos\theta_2\dot{\theta}_2)\simeq\ell(\dot{\theta}_1+\dot{\theta}_2)\)
    \(\dfrac{1}{2}|\dot{x}_2|^2=\dfrac{\ell^2}{2}(\dot{\theta}_1+\dot{\theta}_2)^2\)

    \(\theta_1\,,\,\theta_2\,\)が小さいとするので、運動エネルギーは \(\quad T=\dfrac{M\ell^2}{2}\dot{\theta}_1^2+\dfrac{M\ell^2}{2}(\dot{\theta}_1+\dot{\theta}_2)^2\)

    ポテンシャル・エネルギーは\(\quad U=-Mg\ell(\cos\theta_1+\cos\theta_1+\cos\theta_2)\)

    よって、ラグラジアン\(\,L\,\)は

    \(L=T-U=\dfrac{M\ell^2}{2}(2\dot{\theta}_1^2+2\dot{\theta}_1\dot{\theta}_2+\dot{\theta}^2)+Mg\ell(2\cos\theta_1+\cos\theta_2)\)

    \(\dfrac{\partial L}{\partial\theta_1}\,=-2Mg\ell\sin\theta_1\simeq -2Mg\ell\theta_1\hspace{12mm}\dfrac{\partial L}{\partial\dot{\theta}_1}=M\ell^2(2\dot{\theta}_1+\dot{\theta}_2)\)

    \(\dfrac{\partial L}{\partial\theta_2}\,=-Mg\ell\sin\theta_2\simeq -Mg\ell\theta_2\hspace{16mm}\dfrac{\partial L}{\partial\dot{\theta}_2}=M\ell^2(\dot{\theta}_1+\dot{\theta}_2)\)

    オイラー・ラグランジュ方程式は\(\quad\dfrac{d}{dt}\left(\dfrac{\partial L}{\partial\dot{q}_1}\right)=\dfrac{\partial L}{\partial q}\,\,\)なので

    \(2\ddot{\theta}_1+\ddot{\theta}=-2\dfrac{g}{\ell}\theta_1\qquad \ddot{\theta}_1+\ddot{\theta}_2=-\dfrac{g}{\ell}\theta_2\quad\)の関係が得られる。
    これより
    \(\left(\begin{array}{@{\,}c@{\,}}\ddot{\theta}_1\\[2pt] \ddot{\theta}_2\end{array}\right)=\dfrac{g}{\ell}\left(\begin{array}{@{\,}rr@{\,}}-2&1\\[2pt]2&-2\end{array}\right)\left(\begin{array}{@{\,}c@{\,}}\theta_1\\[2pt] \theta_2\end{array}\right)\qquad\)の連立微分方程式が導かれる。
    \(\hspace{28mm}\)行列\(\,\,K\,\)

    4.方程式の解法

    この行列 \(K\) を対角化することで、方程式を解く。まず \(K\) の固有値と固有ベクトルを求める。
    固有値 \(\lambda\) は、固有ベクトル \(\textbf{x}\) と \(K\textbf{x}=\lambda\textbf{x}\) の関係にあるので、\(\det{(K-\lambda I)}=0\) より
    \(\lambda^2+4\lambda+2=0\qquad \lambda=-2\pm\sqrt{2}\qquad\)が得られる。

    \(\textbf{(4-1)} \,\,\lambda_1=-2+\sqrt{2}\) のとき

    \(\quad(K-\lambda_1 I)=\left(\begin{array}{@{\,}rr@{\,}}-\sqrt{2}&1\\[2pt]2&-\sqrt{2}\end{array}\right)\left(\begin{array}{@{\,}c@{\,}}x_1\\[2pt] x_2\end{array}\right)=0\qquad\) より

    \(\quad x_1=\dfrac{1}{\sqrt{3}}\qquad x_2=\dfrac{\sqrt{2}}{\sqrt{3}}\qquad\quad \textbf{x}_1=\left(\begin{array}{@{\,}r@{\,}}\frac{1}{\sqrt{3}}\\[2pt]\frac{\sqrt{2}}{\sqrt{3}}\end{array}\right)\)

    \(\textbf{(4-2)} \,\,\lambda_2=-2-\sqrt{2}\) のとき

    \(\quad(K-\lambda_2 I)=\left(\begin{array}{@{\,}rr@{\,}}\sqrt{2}&1\\[2pt]2&\sqrt{2}\end{array}\right)\left(\begin{array}{@{\,}c@{\,}}x_1\\[2pt] x_2\end{array}\right)=0\qquad\) より

    \(\quad x_1=\dfrac{1}{\sqrt{3}}\qquad x_2=-\dfrac{\sqrt{2}}{\sqrt{3}}\qquad\quad \textbf{x}_2=\left(\begin{array}{@{\,}r@{\,}}\frac{1}{\sqrt{3}}\\[2pt]-\frac{\sqrt{2}}{\sqrt{3}}\end{array}\right)\)

    5.正則行列の逆行列

    固有ベクトルを縦に並べた行列 \(P\) を用いて \(P^{-1}KP\) として \(K\) を対角化する。

    \(\quad P=\left(\begin{array}{@{\,}rr@{\,}}\frac{1}{\sqrt{3}}&\frac{1}{\sqrt{3}}\\[3pt]\frac{\sqrt{2}}{\sqrt{3}}&-\frac{\sqrt{2}}{\sqrt{3}}\end{array}\right)\qquad\quad P^{-1}=\left(\begin{array}{@{\,}rr@{\,}}\frac{\sqrt{3}}{2}&\frac{\sqrt{3}}{2\sqrt{2}}\\[3pt]\frac{\sqrt{3}}{2}&-\frac{\sqrt{3}}{2\sqrt{2}}\end{array}\right)\)

    この \(P^{-1}\) を用いて、新しい変数 \(\Theta_1\,\,,\,\,\Theta_2\,\,\)を定義する。

    \(\quad\left(\begin{array}{@{\,}c@{\,}}\Theta_1\\[2pt] \Theta_2\end{array}\right)=P^{-1}\left(\begin{array}{@{\,}c@{\,}}\theta_1\\[2pt]\theta_2\end{array}\right)=\left(\begin{array}{@{\,}rr@{\,}}\frac{\sqrt{3}}{2}&\frac{\sqrt{3}}{2\sqrt{2}}\\[3pt]\frac{\sqrt{3}}{2}&-\frac{\sqrt{3}}{2\sqrt{2}}\end{array}\right)\left(\begin{array}{@{\,}c@{\,}}\theta_1\\[2pt]\theta_2\end{array}\right)\)

    元の連立微分方程式は、変数 \(\Theta_1\,\,,\,\,\Theta_2$\,\,\)に対する以下の微分方程式と等価である。

    \(\quad\dfrac{d}{dt}\left(\begin{array}{@{\,}c@{\,}}\Theta_1\\[2pt] \Theta_2\end{array}\right)=\dfrac{g}{\ell}\left(\begin{array}{@{\,}rr@{\,}}-2+\sqrt{2}&0\\[3pt]0&-2-\sqrt{2}\end{array}\right)\left(\begin{array}{@{\,}c@{\,}}\Theta_1\\[2pt] \Theta_2\end{array}\right)\)

    最終的に、運動は2つの固有振動の足し合わせで表され、それぞれの固有振動数 \(\omega_1\,\,,\,\,\omega_2\,\,\)は、

    \(\quad\omega_1=\sqrt{2+\sqrt{2}}\sqrt{\dfrac{g}{\ell}}\quad,\quad \omega_1=\sqrt{2-\sqrt{2}}\sqrt{\dfrac{g}{\ell}}\quad\) となる。







    弾性衝突


    放送大学のオンライン授業で「力学演習」を履修している。そこでの練習問題で一次元の弾性衝突問題があり、すっかり忘れているので再確認のため、計算してみました。

    図のような2球で考える。向きは\(\,x\,\)軸方向に正として、初速度をそれぞれ\(\,v_1\,\,,\,v_2\,\)とする。衝突後の速度をそれぞれ\(\,v_3\,\,,\,v_4\,\)とする。弾性衝突なので、運動量とエネルギーが保存する。

    ・運動量保存
    \(\qquad m_A\,v_1+m_B\,v_2=m_A\,v_3+m_B\,v_4\quad\cdots\,\)(1)
    ・エネルギー保存
    \(\qquad \dfrac{1}{2}m_A\,v_1^2+\dfrac{1}{2}m_B\,v_2^2=\dfrac{1}{2}m_A\,v_3^2+\dfrac{1}{2}m_B\,v_4^2\quad\cdots\,\)(2)

    式(1)から、\(\quad v_4=\dfrac{m_A(v_1-v_3)+m_Bv_2}{m_B}\quad\cdots\,\)(3)
    式(2)から、\(\quad v_4^2=\dfrac{m_A(v_1^2-v_3^2)+m_Bv_2^2}{m_B}\quad\cdots\,\)(4)
    式(3)を2乗すると、
    \(\quad v_4^2=\dfrac{m_A^2v_1^2-2m_A^2v_1v_3+m_A^2v_3^2+m_B^2v_2^2+2m_Am_Bv_1v_2-2m_Am_Bv_2v_3}{m_B^2}\quad\cdots\,\)(5)
    式 (4)\(\,=\,\)(5) なので、
    \(\quad m_Am_Bv_1^2-m_Am_Bv_3^2+m_B^2v_2^2\)
    \(\qquad\qquad =m_A^2v_1^2-2m_A^2v_1v_3+m_A^2v_3^2+m_B^2v_2^2+2m_Am_Bv_1v_2-2m_Am_Bv_2v_3\quad\cdots\,\)(6)

    \(\quad m_A^2v_1^2-2m_Av_1v_3+m_A^2v_3^2+2m_Am_Bv_1v_2-2m_Am_Bv_2v_3-m_Am_Bv_1^2+m_Am_Bv_3^2=0\,\cdots\,\)(7)

    \(\quad(m_A+m_B)v_3^2-2(m_Av_1+m_Bv_2)v_3+m_Av_1^2+2m_Bv_1v_2-m_Bv_1^2=0\quad\cdots\,\)(8)

    この式(8) を\(\,v_3\,\)の2次方程式として解くと
    \(\quad v_3=\dfrac{m_Av_1+m_Bv_2\pm\sqrt{m_B^2(v_1-v_2)^2}}{m_A+m_B}=\dfrac{m_Av1+m_Bv_2\pm m_B|v_1-v_2|}{m_A+m_B}\quad\cdots\,\)(9)

    ここで、球 A の速度が B の速度より大きいケースで考えてみる。 \(\,v_1\ge v_2\)
    \(\,\pm\,\)の\(\,+\,\)では、\(\,v_3=v_1\,\)となり、すり抜ける形なので\(\,-\,\)を採用する。

    \(\quad v_3=\dfrac{m_Av_1-m_Bv_1+2m_Bv_2}{m_A+m_B}\quad\cdots\,\)(10)

    ここで、両方の球の質量が同じ場合、\(\,v_3=v_2\,\)となり、A と B の運動量が入れ替わることになる。

    ◇ 大小球の正面衝突

    大小の球を正面衝突させた場合を考える。
    Aは\(\,m_A=3m\,\)とし、Bは\(\,m_B=m\,\)とする。
    速度はともに\(\,v_1=v_2=v_0\,\)(符号は正負となる)として、正面からの衝突とする。

    式(10)と式(1)により

    \(\quad v_3=\dfrac{3mv_0-mv_0-2mv_0}{3m+m}=0\quad v_4=\dfrac{3mv_0-mv_0}{m}=2v_0\,\)となる。

    結局、Aは停止しBは速度\(\,2v_0\,\)で反対側に跳ね返ることになる。

    ◇ 3球での衝突問題

    『中川のビジュアル物理学教室』のページの中の「解説:一次元衝突」に次のような3球の衝突が紹介されていました。(画像は中川氏のページのものを使わせてもらいました)

    静止している球B,Cに,その2倍の質量をもつ球Aを\(\,v_0\,\)で弾性衝突したらどうなるだろうか?

    この問題は,過去に青学大の入試問題として出題されたことがあるという事でした。

    まず、Aが静止しているBに衝突する時点を考える。式(10) に\(\,m_A=2m\,,\,m_B=m\,,\,v_1=v_0\,,\,v_2=0\,\)を代入すると、
    \(\quad v_A=v_3=\dfrac{2mv_0-mv_0}{3m}=\dfrac{1}{3}v_0\quad\cdots\,\)(11)
    静止していたBの速さ\(\,v_B\,\)は、式(11)の値を式(1)に代入すると\(\,v_B=v_4=\dfrac{4}{3}v_0\,\)となる。
    Bがこの速度\(\,v_B\,\)で、静止しているCに衝突すると
    BとCは同質量なので、式(11)より、Bは停止してCは速度\(\,v_C=v_B\,\)で動き出すことになる。
    停止した(1回動いて停止した形になる)BにAが\(\,v_A=\dfrac{1}{3}v_0\,\)で衝突することになる。
    同じく、式(11)に代入すると、再度のAの速度\(\,v_A’\,\)は
    \(\quad v_A’=\dfrac{2m\frac{1}{3}v_0-m\frac{1}{3}v_0}{3m}=\dfrac{1}{9}v_0\quad\cdots\,\)(12)
    これを式(1)に代入すると\(\quad2m\dfrac{1}{3}=2m\dfrac{1}{9}v_0+mv_B’\,\)より
    Aは\(\,\dfrac{1}{9}v_0\,\)で動き出し、Bは\(\,\dfrac{4}{9}v_0\,\)で動き出す。
    Cは先程、Bの衝突で\(\dfrac{4}{3}v_0\,\)で動き出しているので、3球はバラバラになることになる。







    振り子の実験

    1.経緯

    私が通学(通信制なので、通ってはいないが)している放送大学のゼミで、自宅での実験を動画撮影して皆で共有するという企画があったので、振り子の実験をやってみました。
    振り子は小学\(5\)年生用の実験セットをAmazonから仕入れました。
    (「ふりこのはたらき」 \(730\)円)

    2.実験

    ひもの長さは\(30\)cmとして、振れ角測定用に角度を書いた紙を付けました。大きく振った位置から重りを降ろした。
    (\(2020.09.05\) 実施)
    撮影は iPhone をそのままの設定で(\(30\)f)三脚無しで行いました。時間は約\(2\)分間でした。
    動画は DropBox 経由で、Windows マシンに落とした。再生で簡単にコマ落としが出来、時間の表示が \(1/100\) 秒単位で出てくるので MovieMaker を使用した。(最近は MovieMaker はサポートされていないようですね)

    撮影した動画の抜粋


    編集は mac の iMovie を使用した。縦のプロジェクトが作れないので、ビデオクリップで横にして、編集して、出来た動画を QuickTime のplayer で開いて、縦に変換した。.mov 型の動画になるので、vlc で .mp4 に変換した。
    (このページでは動画のサイズが\(8\)Mb なので抜粋版をつくりました)

    3.実験結果

    動画をMovieMaker でコマ落とし再生して、最大振幅の角度と時間を記録したのが下図です。

    動画が\(30\)fなので、時間刻みが\(1/30\)sec単位となり、グラフは階段状になるが、振れ角が大きいと周期が長いという結果となった。

    4.検討1:振幅の小さい場合

    糸の張力\(\,T\,\)、重力\(\,W\,\)の合力\(\,f^{\prime\prime}\)が重りに掛かる。
    \(\quad W=mg\)
    \(\quad f^{\prime\prime}=mg\sin\theta\)
    円周方向の重りの動く長さを\(\,s\,\)とすると
    \(\quad m\dfrac{d^2s}{dt^2}=-mg\sin\theta\)
    \(\quad s=l\theta\,\)なので\(\quad\dfrac{d^2\theta}{dt^2}=-\dfrac{g}{l}\sin\theta\)
    振幅が小さいときは
    \(\quad\theta\approx\sin\theta\quad\)とすることが出来るので
    \(\quad\dfrac{d^2\theta}{dt^2}=-\dfrac{g}{l}\theta\quad\)とすることが出来、次のような一般解が得られる。

    \(\quad\theta=A\sin\sqrt{\dfrac{g}{l}}t+B\cos\sqrt{\dfrac{g}{l}}t\)
    ここで初期条件として\(\quad t=0\quad d\theta/dt=0\quad\theta=\theta_0\quad\)とすると
    \(\quad\theta=\theta_0\cos\sqrt{\dfrac{g}{l}}t\quad\)が得られる。
    周波数\(\,f\,\)は\(\quad f=\dfrac{1}{2\pi}\sqrt{\dfrac{g}{l}}\quad\)となる。
    周期\(\,T\,\)は\(\quad T=\dfrac{1}{f}=2\pi\sqrt{\dfrac{l}{g}}\quad\)である。

    今回は紐の長さを\(\,30\,\)cm としたので、

    \(\quad T=2\times 3.141593\times\sqrt{\dfrac{0.3}{9.80655}}=1.099..\,\)(s) となった。結果のグラフを見ると、振れ角\(20\,\)度以下ではほぼこの周期となっているようである。

    5.検討2:振幅が大きい場合

    エネルギー保存則を用いて、振幅の大きい振り子の運動を計算する。振り子の最下点を原点とし、振りだし高さを\(\,h_0\,\)、\(\Delta h\,\)だけ降りた点の円周方向の速度を\(\,v\,\)とすると
    \(\quad\dfrac{1}{2}mv^2=mgl(\cos\theta-\cos\theta_0)\quad\)なので
    \(\quad\therefore v=\sqrt{2gl(\cos\theta-\cos\theta_0)}\)
    ここで、\(\,v=\dfrac{ds}{dt}=l\dfrac{d\theta}{dt}\quad\)なので
    \(\quad l\dfrac{d\theta}{dt}=\sqrt{2gl(\cos\theta-\cos\theta_0)}\quad\)より
    \(\quad dt=\sqrt{\dfrac{l}{2g}}\dfrac{d\theta}{\sqrt{\cos\theta-\cos\theta_0}}\)

    これを積分すると周期を求めることが出来る。
    \(\,t\,\)の積分範囲を\(\,4\,\)分の\(\,1\,\)周期とすると、\(\theta=0\,\to\,\theta_0\)
    \(\displaystyle\int^{\frac{T}{4}}_0dt=\dfrac{T}{4}=\sqrt{\dfrac{l}{2g}}\displaystyle\int_0^{\theta_0}\dfrac{d\theta}{\sqrt{\cos\theta-\cos\theta_0}}\)

    ここからの積分は振り子の計算その2を参考にしてください。

    プロセスとしては\(\quad\sin\phi=\sin(\theta/2)/\sin(\theta_0/2)\quad\)とおいて、いくつかの計算を経て以下の積分となる。
    \(\quad T=4\sqrt{\dfrac{l}{g}}\displaystyle\int_0^{\frac{\pi}{2}}\dfrac{d\phi}{\sqrt{1-\sin^2(\theta_0/2)\sin^2\phi}}\)
    この積分は楕円積分と呼ばれる。

    5.1 楕円積分の計算

    楕円積分を参照のこと。

    第1種の楕円積分は以下のような式となる。

    \(\quad I=\displaystyle\int_0^{\frac{\pi}{2}}\dfrac{d\phi}{\sqrt{1-a^2\sin^2\phi}}=\dfrac{\pi}{2}\displaystyle\sum_{n=0}^{\infty}\left\{\dfrac{2n-1)!!}{2n!!}\right\}^2a^{2n}\)
    ここで\(\quad n!!\quad\)は二重階乗である。

    これより、\(\,T\quad\)は以下のようになる。
    \(\quad T=4\sqrt{\dfrac{l}{g}}\displaystyle\int_0^{\frac{\pi}{2}}\frac{d\phi}{\sqrt{1-\sin^2(\theta_0/2)\sin^2\phi}}\)

    \(\qquad=2\pi\sqrt{\dfrac{l}{g}}\displaystyle\sum_{n=0}^{\infty}\left\{\dfrac{(2n-1)!!}{2n!!}\right\}^2\left(\sin\dfrac{\theta_0}{2}\right)^{2n}\)

    具体的な計算は\(\quad b=\sin\dfrac{\theta_0}{2}\quad\)とおいて、振り子の長さ\(\,l=0.3\,\)m を使うと

    \(\quad T=2\pi\sqrt{\dfrac{l}{g}}\left\{1+\left(\dfrac{1}{2}\right)^2b^2+\left(\dfrac{1}{2}\dfrac{3}{4}\right)^2b^4+\left(\dfrac{1}{2}\dfrac{3}{4}\dfrac{5}{6}\right)^2b^6\right.\)

    \(\qquad\qquad\qquad \left.+\left(\dfrac{1}{2}\dfrac{3}{4}\dfrac{5}{6}\dfrac{7}{8}\right)^2b^8+\left(\dfrac{1}{2}\dfrac{3}{4}\dfrac{5}{6}\dfrac{7}{8}\dfrac{9}{10}\right)^2b^{10}+\cdots\,\right\}\)

    \(\qquad=1.099\times\left\{1+\dfrac{1}{4}b^2+\dfrac{9}{64}b^4+\dfrac{25}{256}b^6+\dfrac{1225}{16384}b^8\right.\)

    \(\qquad\qquad\qquad\left.+\dfrac{3969}{65536}b^{10}+\dfrac{53361}{1048576}b^{12}+\cdots\,\right\}\)

    この値と実験結果を併せて表示する。







    測定データが \(\quad 1/30\,\)sec 単位であるが、まあ妥当な結果となったようである。

    積分その1

    積分その1

    【1】\(\,\,\displaystyle\int_0^{\infty}xe^{-ax}dx\quad\)の計算

    \(\quad\displaystyle\int_0^{\infty}xe^{-ax}dx\quad\cdots\,\)(1)

    まず、\(\,\displaystyle\int xe^{-ax}dx\,\,\)の不定積分から計算する。

    部分積分で考えると

    \(\quad\displaystyle\int xe^{-ax}dx=x(-ae^{-ax})+\dfrac{1}{a}\displaystyle\int e^{-ax}dx\)
    \(\qquad=-axe^{-ax}+\dfrac{1}{a}\displaystyle\int e^{-ax}dx\quad\cdots\,\)(2)

    (1)式の定積分を考える時、まず次の極限を計算する。

    \(\quad\displaystyle\lim_{x\to\infty}\dfrac{x}{e^x}\quad\cdots\,\)(3)

    \(\quad x\ge 0\,\,\)のとき\(\,\,e^x\ge\dfrac{x^2}{2}\,\,\)を考える。

    \(\quad f(x)=e^x-\dfrac{x^2}{2}\quad\cdots\,\)(4)

    とすると、\(x\ge 0\quad\)のとき、\(f'(x)\gt 0\,\,,\quad f^{\prime\prime}\gt 0\quad\)なので、\(f(x)\gt 0\quad\)である。よって

    \(\,\,e^x\ge\dfrac{x^2}{2}\,\,\)が成り立つ。そこで

    \(\quad 0\ge \dfrac{x}{e^x}\ge \dfrac{x}{\frac{x^2}{2}}\quad\)つまり\(\quad 0\ge \dfrac{x}{e^x}\ge \dfrac{2}{x}\quad\)となり、

    \(\quad\displaystyle\lim_{x\to\infty}\dfrac{2}{x}=0\quad\)なので\(\quad\displaystyle\lim_{x\to\infty}\dfrac{x}{e^x}=0\quad\cdots\,\)(5)

    となる。

    この結果 (2)式を\(\,0\,\)から\(\,\infty\,\)まで積分すると

    \(\quad\displaystyle\int_0^{\infty}xe^{-ax}dx=\left[-axe^{-ax}\right]_0^{\infty}+\dfrac{1}{a}\displaystyle\int_0^{\infty}e^{-ax}dx\)

    \(\qquad =0-0+\dfrac{1}{a}\left[-\dfrac{1}{a}e^{-ax}\right]_0^{\infty}=\dfrac{1}{a^2}\quad\cdots\,\)(6)

    グラフは\(\,y=xe^{-x}\,\)である。\(x\,\to\,\infty\,\)まで積分すると斜線部分の面積は\(\,1\,\)となります。







    ビリアル定理

    1-1. クラウジウスのビリアル定理

    ビリアル定理は 19世紀に Clausius によって考案された。
    「系の平均活力は、その(平均)ビリアル(の大きさ)に等しい。」
    ここで登場する「活力」(”vis viva”)は、今日の運動エネルギーに相当します。 そもそも「活力」は古典力学の草創期にライプニッツが導入した量で、今日の運動エネルギーの2倍に相当する量 \(\,mv^2\,\)ですが、クラウジウスは \(\,mv^2/2\,\)とちょうど運動エネルギーと同じ量として用いています。ビリアル virial はラテン語のvis(「力」)からクラウジウスが作った造語です。ここでは、(力)×(位置ベクトル)をビリアルと呼ぶことにします。

    距離の逆2乗則に従う重力クーロン力の中心力で相互作用しあっている多体系では、長時間平均した運動エネルギー\(\,\langle\,T\,\rangle\,\)と平均の全ポテンシャルエネルギー\(\,\langle\,V\,\rangle\,\)は次の関係式を満たす。

    \(\quad 2\langle\,T\,\rangle+\langle\,V\,\rangle=0\quad\cdots\,\)(1)

    1-2. ビリアル定理の古典力学的証明

    ここで位置ベクトル\(\,\vec{r}\,\)と運動量\(\,\vec{p}\,\)の内積の総和を以下のように考える。

    \(\quad G=\displaystyle\sum_i\vec{p}_i\,\cdot\,\vec{r}_i\quad\cdot\,\)(2)

    (2)式を時間\(\,t\,\)で微分する。

    \(\quad\dfrac{dG}{dt}=\displaystyle\sum_i\vec{p}_i\cdot\dfrac{d\vec{r}_i}{dt}+\displaystyle\sum_i\dfrac{d\vec{p}_i}{dt}\cdot\vec{r}_i\)
    \(\qquad\quad=\displaystyle\sum_i\,m_i(\vec{v}_i)^2+\displaystyle\sum_i\vec{F}_i\cdot\vec{r}_i\)
    \(\qquad\quad=\displaystyle\sum_i2\times\dfrac{1}{2}m_i(\vec{v}_i)^2+\displaystyle\sum_i\vec{F}_i\cdot\vec{r}_i\)
    \(\qquad\quad=2T+\displaystyle\sum_i\vec{F}_i\cdot\vec{r}_i\quad\cdots\,\)(3)

    ここでは、次の関係式を使用している。

    \(\quad \vec{p}_i=m_i\vec{v}_i\,,\qquad\vec{v}_i=\dfrac{d\vec{r}_i}{dt}\,,\qquad\vec{F}_i=\dfrac{d\vec{p}_i}{dt}\quad\cdots\,\)(4)

    (3)式の両辺を\(\,0\,\)から時間\(\,t\,\)の範囲で積分して\(\,t\,\)で割り、\(t\to\infty\,\)の極限をとって長時間平均する。すると粒子が動き得る範囲は有限なので\(\,G\,\)も有限だから、左辺は\(\,0\,\)に収束する。

    \(\quad\displaystyle\lim_{t\to\infty}\dfrac{1}{t}\displaystyle\int_0^t\dfrac{dG}{dt}=\displaystyle\lim_{t\to\infty}\dfrac{G(t)-G(0)}{t}=0\quad\cdots\,\)(5)

    したがって、

    \(\quad 0=2\langle\,T\,\rangle+\left\langle\displaystyle\sum_i\vec{F}_i\cdot\vec{r}_i\right\rangle\quad\cdots\,\)(6)

    つまり、ビリアル定理を得る。

    \(\quad\langle\,T\,\rangle=-\dfrac{1}{2}\left\langle\displaystyle\sum_i\vec{F}_i\cdot\vec{r}_i\right\rangle\quad\cdots\,\)(7)

    次にポテンシャルエネルギー\(\,V\,\)が中心力ポテンシャルで、粒子間の距離に反比例する形で、系のポテンシャル\(\,V\,\)が各粒子対の相互作用の和によって書き表される場合、以下のように表される。

    \(\quad V(r_1,\cdots,r_N)=\displaystyle\sum_{i\lt j}\dfrac{a_{ij}}{|\vec{r}_i-\vec{r}_j|}\quad\cdots\,\,\)(8)

    粒子\(\,i\,\)に作用する力の合計\(\,\vec{F}_i\,\)は、ポテンシャルを位置座標で微分することで以下のように表すことが出来る。

    \(\quad\vec{F}_i=-\nabla_{r_i}V=\displaystyle\sum_{j\ne i}\dfrac{a_{ij}(\vec{r}_i-\vec{r}_j)}{|\vec{r}_i-\vec{r}_j|^3}=\displaystyle\sum_{j\ne i}\vec{F}_{ij}\quad\cdots\,\)(9)

    ここで、\(\,\vec{F}_{ij}\,\)は、粒子\(\,j\,\)から粒子\(\,i\,\)に働く力である。

    \(\quad\vec{F}_{ij}=a_{ij}\dfrac{\vec{r}_i-\vec{r}_j}{|\vec{r}_i-\vec{r}_j|^3}\quad\cdots\,\)(10)

    この\(\,\vec{F}_{ij}\,\)を用いると、ビリアル(力と位置ベクトルの内積)は次のように表せる。

    \(\quad\displaystyle\sum_i\vec{F}_i\cdot\vec{r}_i=\displaystyle\sum_{i,j(i\lt j)}\vec{F}_{ij}\cdot\vec{r}_i+\displaystyle\sum_{i,j(j\lt i)}\vec{F}_{ij}\cdot\vec{r}_i\)
    \(\qquad\qquad\quad=\displaystyle\sum_{i,j(i\lt j)}\vec{F}_{ij}\cdot\vec{r}_i+\displaystyle\sum_{i,j(i\lt j)}\vec{F}_{ji}\cdot\vec{r}_j\qquad\because\,i\leftrightarrow j\)
    \(\qquad\qquad\quad=\displaystyle\sum_{i,j(i\lt j)}\vec{F}_{ij}\cdot(\vec{r}_i-\vec{r}_j)\quad\because\,\vec{F}_{ji}\!=\!-\vec{F}_{ij}\,\,\cdots\,\)(11)

    (10)式を(11)式に代入すると

    \(\quad\displaystyle\sum_i\vec{F}_i\cdot\vec{r}_i=\!\displaystyle\sum_{i,j(i\lt j)}\dfrac{a_{ij}}{|\vec{r}_i-\vec{r}_j|}\)
    \(\qquad\qquad\qquad=\displaystyle\sum_{i,j(i\lt j)}\!V_{ij}=\!\displaystyle\!\sum_{i,j(i\lt j)}\!V_{ji}\quad\cdots\,\)(12)

    よって(6)式は以下のように表すことが出来る。

    \(\quad 0=2\langle\,T\,\rangle+\left\langle\displaystyle\sum_i\vec{F}_i\cdot\vec{r}_i\right\rangle=2\langle\,T\,\rangle+\left\langle\displaystyle\sum_{i,j(i\lt j)}V_{ij}\right\rangle\)
    \(\qquad\qquad=2\langle\,T\,\rangle+\langle\,V\,\rangle\quad\cdots\,\)(13)







    ドメイン変更

    tamami.tech というドメインを取得したけど、更新料が 3,400 円と高額なので、tamami8.net という 更新料 1,480 円のドメインに変更しました。
    まあ変更と行っても、”*.tech” のドメインは 2021年の6月まで有効ですが、HTML のソースを順次書き換え中です。

    最近、ロリポップではサーバーレンタルとドメイン取得をまとめてやると、更新料無料らしい。まあ、しょうがないと諦めています。

    懸垂線


    ひもの両端を持って垂らしたときにできる曲線を懸垂曲線( カテナリー catenary )といいますが、この曲線の方程式を導いてみる。

    ひもの底を原点\(O\,\)とし、水平方向に\(x\,\)軸、鉛直方向に\(y\,\)軸をとる。ひもの曲線を\(\,y=f(x)\,\)とおいて関数\(\,f(x)\,\)を求める。

    ひもの\(\,x>0\,\)の部分に点\(\,(x,f(x))\,\)をとり、これを\(P\,\)とします。また、ひもの\(O\,\)から\(P\,\)の部分(両端を含む)を\(C\,\)とします。

    【\(\,f(x)\,\)に対する微分方程式を立てる】

    ひもの一部\(\,C\,\)に加わる力の釣り合いの式から\(\,f(x)\,\)についての微分方程式を立てる。\(C\,\)には次の3つの力が加わっている。

    • 点\(\,O\,\)に水平方向のひもの張力\(\,T_0\)
    • 点\(\,P\,\)に接線方向上向きにひもの張力\(\,T\)
    • \(\,C\,\)の重心に鉛直下向きに重力\(\,G\)

    ここで、原点にはたらく張力\(\,T_0\,\)は\(\,x\,\)に関係なく一定なので定数とする。

    \(T\,\)について、\(P\,\)における\(\,f(x)\,\)の接線と\(\,x\,\)軸の正の部分とのなす角を\(\,\theta\,\)とすると、\(T\,\)の\(\,x\,\)成分、\(y\,\)成分はそれぞれ\(\,T\cos\theta\,,\,T\sin\theta\,\)となる。

    この角度\(\,\theta\,\)と\(\,f(x)\,\)は、\(x\,\)における微分が接線なので
    \(\quad f'(x)=\tan\theta\quad\)となる。

    \(C\,\)の長さを\(\,l\,\)とすると曲線の長さの式から
    \(\quad l=\displaystyle\int_0^x\sqrt{1+\{f'(x)\}^2}dx\quad\)で求められる。

    線密度を\(\,\rho\,\)、重力加速度を\(\,g\,\)とすると
    \(\quad G=\rho lg=\rho g\quad\)となる。

    これらより釣り合いの式は以下のようになる。

    • 鉛直方向:\(\,T\sin\theta=\rho lg\)
    • 水平方向:\(\,T\cos\theta=T_0\)

    \(T\,\)を消去して、定数をまとめて \(\,k=\dfrac{\rho g}{T_0}\,\,\)とおくと

    \(\quad f'(x)=k\displaystyle\int_0^x\sqrt{1+\{f'(x)\}^2}dx\quad\)となる。

    両辺を\(\,x\,\)で微分すると

    \(\quad f^{\prime\prime}(x)=k\sqrt{1+\{f'(x)\}^2}\quad\)となる。

    境界条件は\(\quad f(0)=0\,\,,\quad f'(0)=0\quad\)となる。

    【微分方程式を解く】

    この微分方程式を解く。\(\quad z=f'(x)\quad\)とおくと

    \(\quad\dfrac{dz}{dx}=k\sqrt{1+z^2}\)
    \(\quad\dfrac{dz}{\sqrt{1+z^2}}=k\,dx\quad\)となる。

    先程の条件で\(\quad f'(x)=\tan\theta=z\quad\)なので
    \(\quad dz=\dfrac{1}{\cos^2\theta}d\theta\quad\)となる。よって、

    \(\quad \dfrac{dz}{\sqrt{1+z^2}}=\dfrac{1}{\sqrt{1+\tan^2\theta}}\dfrac{d\theta}{\cos^2\theta}=\dfrac{d\theta}{\sqrt{\dfrac{\cos^2\theta+\sin^2\theta}{\cos^2\theta}}\cos^2\theta}\)
    \(\qquad =\dfrac{d\theta}{\cos\theta}\quad\)となる。

    このままでは積分出来ないので、\(\quad\dfrac{d}{d\theta}\sin\theta=\cos\theta\quad\)を利用して

    \(\quad\dfrac{d\theta}{\cos\theta}=\dfrac{d(\sin\theta)}{\cos^2\theta}=\dfrac{d(\sin\theta)}{(1-\sin\theta)(1+\sin\theta)}\)

    \(\quad =\dfrac{1}{2}\left(\dfrac{1}{1-\sin\theta}+\dfrac{1}{1+\sin\theta}\right)\,d(\sin\theta)\quad\)となる。

    この形なら積分可能なので\(\,\,\sin\theta\,\,\)で積分する。

    \(\quad\displaystyle\int\dfrac{dz}{\sqrt{1+z^2}}=\dfrac{1}{2}\displaystyle\int\left(\dfrac{1}{1-\sin\theta}+\dfrac{1}{1+\sin\theta}\right)\,d(\sin\theta)\)

    \(\quad=\dfrac{1}{2}\left\{-\log(1-\sin\theta)+\log(1+\sin\theta)\right\}\)

    \(\quad=\dfrac{1}{2}\log\left(\dfrac{1+\sin\theta}{1-\sin\theta}\right)\)

    よって\(\quad z=\tan\theta\quad\)として

    \(\quad f'(x)=\displaystyle\int\dfrac{dz}{\sqrt{1+z^2}}=\dfrac{1}{2}\log\left(\dfrac{1+\sin\theta}{1-\sin\theta}\right)=kx+C_1\quad\) となる。

    境界条件\(\quad f'(x)=0\quad x=0\,\,,\,\theta=0\quad\)より\(\quad C_1=0\,\,\)。

    \(\quad\dfrac{1+\sin\theta}{1-\sin\theta}=e^{2kx}\quad\)より、まず\(\quad\sin\theta\quad\)について解く。

    \(\quad 1+\sin\theta=e^{2kx}(1-\sin\theta)\qquad (e^{2kx}+1)\sin\theta=e^{2kx}-1\)

    \(\quad \sin\theta=\dfrac{e^{2kx}-1}{2^{2kx}+1}\quad\)となる。ここから以下の公式を使う。

    \(\quad 1+\dfrac{1}{\tan^2\theta}=\dfrac{\tan^2\theta+1}{\tan^2\theta}=\dfrac{\sin^2\theta+\cos^2\theta}{\sin^2\theta}=\dfrac{1}{\sin^2\theta}\)

    \(\quad\dfrac{1}{\tan^2\theta}=\left(\dfrac{e^{2kx}+1}{2^{2kx}-1}\right)^2-1\)

    \(\qquad=\dfrac{(e^{2kx}+1)^2-(e^{2kx}-1)^2}{(e^{2kx}-1)^2}\)

    \(\qquad=\dfrac{4e^{2kx}}{(e^{2kx}-1)^2}\quad\)より

    \(\quad\tan\theta=\dfrac{e^{2kx}-1}{2e^{kx}}=\dfrac{e^{kx}-e^{-kx}}{2}=\sinh kx\)

    \(\quad f'(x)=z=\tan\theta=\sinh kx\quad\)なので

    \(\quad f(x)=\displaystyle\int\dfrac{e^{kx}-e^{-kx}}{2}dx=\dfrac{e^{kx}+e^{-kx}}{2k}+C_2\)

    境界条件が\(\quad f(0)=0\quad\)なので\(\quad C_2=-k\quad\)より

    \(\quad f(x)=\dfrac{e^{kx}+e^{-kx}}{2k}-k=\dfrac{\cosh kx-1}{k}\)







    CO回転スペクトル

    1. ミクロの世界の量子状態

    1-1. 箱の中の粒子

    一次元の粒子の波動方程式を考える。\(\,x\,\)軸上に、\(\,x=0\,\)から\(\,x=L\,\)までの長さが\(\,L\,\)の区間(\(\,0\ge\,x\ge\,L\,\))を考え、この区間の内部でだけ粒子が運動できるモノとする。このような状況の粒子を、1次元の箱の中の粒子と呼ぶ。
    箱の外に出られないように、箱の外の位置エネルギーを\(\,\infty\,\)とする。また、箱の中では粒子に力が働かず自由に運動できるように、位置エネルギーを\(\,0\,\)とする。
    解くべき方程式は、定常状態の波動方程式 \(\,\hat{H}\phi(x)=E\phi(x)\,\) であり、粒子の質量を\(\,m\,\)とするとハミルトニアン\(\,\hat{H}\,\)は以下のようになる。

    \(\quad \hat{H}=-\dfrac{\hbar^2}{2m}\Delta+U(x)\quad\cdots\,\)(1.1)

    変数は\(\,x\,\)のみの1次式なので、ラプラシアン\(\,\Delta\,\)は\(\,x\,\)の2次微分 \(\,d^2/dx^2\,\)に等しい。箱の外では位置エネルギーが無限大だから、\(\,\hat{H}\phi=E\phi\,\)が成り立つためには、\(\phi=0\,\)にならなければならない。
    すなわち、箱の外の確率密度は\(\,0\,\)になり、粒子は箱の外には存在しない。次に箱の中では、位置エネルギーを\(\,U(x)=0\,\)としたので、波動方程式は次のような2次微分を含む方程式になる。

    \(\quad -\dfrac{\hbar^2}{2m}\ \dfrac{d^2\phi}{dx^2}=E\phi\quad\cdots\,\)(1.2)

    この式では\(\,\phi\,\)と\(\,x\,\)以外はすべて定数であり、次の形の基本的な微分方程式となっている

    \(\quad \dfrac{d^2\phi}{dx^2}=-k^2\phi\quad\cdots\,\)(1.3)

    ここで右辺の\(\,k\,\)は、求めるべきエネルギー\(\,E\,\)を含み、次式で与えられる定数である

    \(\quad k=\left(\dfrac{2mE}{\hbar^2}\right)^{1/2}\quad\cdots\,\)(1.4)

    式 (1.3)の一般解は指数関数または三角関数を含む次式の形に表される

    \(\quad \phi(x)=ae^{ikx}+be^{-ikx}=A\sin kx+B\cos kx\quad\cdots\,\)(1.5)

    波動関数の重要な性質として、連続性の条件があり、境界条件を考慮しなければならない。
    箱の外では\(\,\phi=0\,\)であるから、箱の中(\(\,0\ge x\ge L\,)\)の\(\,\phi\,\)も箱の両端で値が\(\,0\,\)にならなくてはいけない。
    そこで境界条件として \(\,\phi(0)=\phi(L)=0\,\) を使うと、\(\,x=0\,\)で\(\,\cos\,\)関数の部分は\(\,0\,\)にならないので、\(\sin\,\)の関数部分だけが残る。
    また、\(\,x=L\,\)で\(\,\sin kx\,\)の値が\(\,0\,\)になるには \(\,kL=n\pi\,\)となり、その結果、波動関数とエネルギーがそれぞれ次のように決まる。

    \(\quad\phi(x)=A\sin(n\pi x/L)\quad\cdots\,\)(1.6)

    \(\quad E=\dfrac{(n\pi\hbar/L)^2}{2m}\quad\cdots\,\)(1.7)

    ここで整数\(\,n\,\)は\(\,1,2,3\,\)などの自然数であり、これは量子数と呼ばれる。

    (1.7) 式から、箱の中の粒子に許されるエネルギーは、図のように飛び飛びのエネルギー準位になる。

    この結果は有限な空間に閉じ込められた粒子のエネルギーが飛び飛びになって量子化されることを示している。
    ここで注目するべきことは、エネルギーが最低の基底状態、すなわち\(\,n=1\,\)の状態でも、エネルギーが\(\,0\,\)にならないことである。このエネルギーをゼロ点エネルギーという。

    (1.7) 式で与えられる波動関数は図の中央に示したように、エネルギー準位を区別する量子数\(\,n\,\)が増えるにつれ上下の振動が激しくなる。波動関数の値が\(\,0\,\)になる位置を節(ふし)という。

    節になる位置は、\(\phi\,\)でも\(\,\phi^2\,\)も同じであり、図の右に示した\(\,\phi^2\,\)の図から明らかなように、節の位置に粒子が観測される確率は\(\,0\,\)である。

    1-2. 調和振動子

    波動方程式の簡単な例として、左図のようなバネ(バネ定数を\(\,k\,\))につるされた質量\(\,m\,\)の錘(おもり)の振動を考える。位置エネルギーは、\(\,U(x)=\frac{1}{2}kx^2\,\) であり、釣り合った位置から\(\,x\,\)だけ錘がずれると、図の右に示した放物線に沿って上がって行き、常に中心に引き戻される力が働いて振動する。このような振動を調和振動子という。

    この場合のハミルトニアンの運動エネルギーは箱の中の粒子の場合と同じになり、位置エネルギーは\(\,U(x)=\frac{1}{2}kx^2\,\) であるから、この調和振動子の波動方程式は次のように表される。

    \(\left(-\dfrac{\hbar^2}{2m}\dfrac{d^2}{dx^2}+\dfrac{1}{2}kx^2\right)\Psi=E\Psi\quad\cdots\,\)(1.8)

    この微分方程式は、次のように波動関数を仮定すると解くことができる。

    \(\quad\Psi(x)=e^{-ax^2}\,\displaystyle\sum_{i=0}b_ix^i\quad\cdots\,\)(1.9)

    解いた結果、次式に従うエネルギー準位が得られる。

    \(\quad E_n=\left(n+\dfrac{1}{2}\right)h\nu\qquad(n=0,1,2,\cdots)\quad\cdots\,\)(1.10)

    ここで、\(\,n=0,1,2,\cdots\,\) は振動の量子数(振動量子数)であり、\(h\,\)はプランク定数で、\(\nu\,\)はこの振動の固有振動数である。


    調和振動子に許されるエネルギー準位は図の左に示すように等間隔になる。その間隔は、振動のエネルギー量子\(\,h\nu\,\)であり,固有振動数\(\,\nu\,\)はバネの力の定数が\(\,k\,\)で質量が\(\,m\,\)の古典力学の振動子と全く同じで、次式によって与えられる。

    \(\quad\nu=\dfrac{1}{2\pi}\sqrt{\dfrac{k}{m}}\quad\cdots\,\)(1.11)

    量子論の振動子には、古典力学と明らかに違った注目すべき特徴がある。エネルギーが最低の状態(基底状態)でも、そのエネルギーは\(\,0\,\)にならず、\(\,\dfrac{1}{2}h\nu\,\)というエネルギーをもって振動する。このエネルギーをゼロ点エネルギーといい、その時の振動をゼロ点振動という。

    図の左に示したように、調和振動子の波動関数は中心の周りで振動子、遠くに行くと減衰する。値が\(\,0\,\)になる節の数は量子数\(\,n\,\)に等しく、エネルギーが高いほど多くなる。

    1-3. 2粒子系の運動方程式

    2個の粒子を含む場合すなわち2粒子系の波動方程式は、以下のように取り扱うと1個の粒子の問題になり簡単になる。図の左に示したように、質量\(\,m_1\,\)と\(\,m_2\,\)の2個の球が距離$r$で結ばれているとして、その運動を考える。
    ここで重心を座標原点に固定すると、\(\,r_1\,\)と\(\,r_2\,\)の和が\(\,r\,\)になり、2個の粒子の相対運動は長さが\(\,r\,\)の棒の先端につけた1つの球の運動の問題に単純化される。そのとき先端の粒子の質量\(\,\mu\,\)は、2個の粒子の質量の逆数の和の逆数になる。

    \(\quad\)換算質量\(\quad\mu=\dfrac{1}{\frac{1}{m_1}+\frac{1}{m_2}}\)

    こうして2個の粒子の相対運動は、換算質量をもつ1つの粒子の運動として取り扱うことが出来る。したがって、その波動方程式は次のように、1つの粒子の方程式と同じになる。

    \(\quad\left(-\dfrac{\hbar^2}{2\mu}\Delta+U\right)\Psi=E\Psi\quad\cdots\,\)(1.13)

    このような問題では、原点からの距離\(\,r\,\)が特別に重要なので、直交座標系\(\,(x,y,z)\,\)の代わりに、図の右に示す3次元の極座標\(\,(r,\theta,\varphi\,)\,\)を変数にとる。すると体積要素\(\,dxdydz\,\)とラプラシアン\(\,\Delta\,\)は、以下のように表される。

    \(\quad\)体積要素\(\quad dxdydz=r^2\sin\theta dr d\theta d\varphi\quad\cdots\,\)(1.14)

    \(\quad\)ラプラシアン\(\quad \Delta=\dfrac{1}{r^2}\dfrac{\partial}{\partial r}\left(r^2\dfrac{\partial}{\partial r}\right)+\dfrac{1}{r^2}\Lambda\quad\cdots\,\)(1.15)

    \(\quad\)ルジャンドリアン\(\quad \Lambda=\dfrac{1}{\sin\theta}\dfrac{\partial}{\partial\theta}\left(\sin\theta\dfrac{\partial}{\partial\theta}\right)+\dfrac{1}{\sin^2\theta}\dfrac{\partial^2}{\partial\varphi^2}\quad\)(1.16)

    ここで、(1.15)式のラプラシアンに含まれる\(\,\Lambda\,\)は角度に依存する演算子であり、ルジャンドリアンと呼ばれる。

    1-4. 剛体回転子

    2量子系の波動方程式は、換算質量\(\,\mu\,\)を使うと原点から距離\(\,r\,\)のところに繋がれた1個の粒子の回転運動になる。距離\(\,r\,\)が一定なので\(\,r\,\)に関する微分は消えてしまい、また位置エネルギーも\(\,0\,\)とおいてよいので、その波動方程式は次のように単純な式になる。

    \(\quad-\dfrac{\hbar^2}{2I}\Delta\Psi=E\Psi\quad\cdots\,\)(1.17)

    ここで、次式で表される\(\,I\,\)は、重心を軸とする回転運動の慣性モーメントである。

    \(\quad I=\mu r^2\quad\cdots\,\)(1.18)

    (1.17)式の波動方程式は三角関数を組み合わせた式を仮定して解くことが出来、その結果、剛体回転子に許されるエネルギーは次式のように、整数値をとる回転量子数\(\,J\,\)によって飛び飛びの準位になる。

    \(\quad E=\dfrac{J(J+1)\hbar^2}{2I}\qquad (J=0,1,2,\cdots)\quad\cdots\,\)(1.19)

    剛体回転子の波動方程式は、数学でよく知られている球面調和関数\(\,Y_{J,m}(\theta,\varphi)\,\)というものになり、それは\(\,J\,\)のほかにもう一つの量子数\(\,m\,\)にも依存している。\(\,J\,\)が決まれば、それに従って\(\,m\,\)のとりうる範囲が決まる。\(\,m\,\)は\(\,J\,\)から\(\,J-1,J-2,\cdots,0,\cdots,-J\,\)まで、合計\(\,2J+1\,\)通りの値が可能である。

    \(\quad Y_{00}(\theta,\varphi)=\left(\dfrac{1}{4\pi}\right)^{1/2}\quad\cdots\,\)(1.20)
    \(\quad Y_{10}(\theta,\varphi)=\left(\dfrac{3}{4\pi}\right)^{1/2}\cos\theta\quad\cdots\,\)(1.21)
    \(\quad Y_{1\pm1}(\theta,\varphi)=\mp\left(\dfrac{3}{8\pi}\right)^{1/2}\sin\theta\ e^{\pm i\varphi}\quad\cdots\,\)(1.22)
    \(\quad Y_{20}(\theta,\varphi)=\left(\dfrac{5}{16\pi}\right)^{1/2}(3\cos^2\theta-1)(\quad\cdots\,\)(1.23)
    \(\quad Y_{2\pm1}(\theta,\varphi)=\mp\left(\dfrac{15}{8\pi}\right)^{1/2}\sin\theta\ \cos\theta\ e^{\pm i\varphi}\quad\cdots\,\)(1.24)
    \(\quad Y_{2\pm2}(\theta,\varphi)=\mp\left(\dfrac{15}{32\pi}\right)^{1/2}\sin^2\theta\ e^{\pm i2\varphi}\quad\)(1.25)

    2.分子の回転

    2-1. 回転運動のエネルギー

    質量\(\,m\,\)、速度\(\,v\,\)の直線運動の大きさが運動量\(\,p=mv\,\)で、運動に伴うエネルギーは

    \(\quad E=\dfrac{1}{2}mv^2=\dfrac{1}{2}\dfrac{p^2}{m}\quad\cdots\,\)(2.1)

    である。一方、回転運動の大きさを表す量が角運動量であり、回転の軸から距離\(\,r\,\)にあって、運動量\(\,p\,\)で運動している物体の角運動量\(\,p_{\theta}\,\)は、

    \(\quad p_{\theta}=r\times p\quad\cdots\,\)(2.2)

    と表される。\(\,p_{\theta}\,\)は運動量を含むベクトルである。この角運動量を用いて、回転運動のエネルギーは

    \(\quad E=\dfrac{1}{2}\dfrac{p_{\theta}^2}{I}\quad\cdots\,\)(2.3)

    のように表される。式(2.1)と(2.3)を比較すると対応していることが分かる。

    直線運動を行う物体で運動の持続性(慣性)に影響する量である質量に対応するのが、回転運動における\(\,I\,\)であり、慣性モーメントと呼ばれる。慣性モーメントは

    \(\quad I=\sum_im_ir_i^2=\sum_im_i(x_i^2+y_i^2+z_i^2)\quad\)(2.4)

    であるが、2原子分子では換算質量\(\,\mu\,\)と原子間の距離\(\,r\,\)を使うと、以下の簡単な式で表すことが出来る

    \(\quad \mu=\dfrac{m_1\cdot m_2}{m_1+m_2}\quad\cdots\,\)(2.5)

    \(\quad I=\mu r^2\quad\cdots\,\)(2.6)

    式(2.3)には運動量\(\,p_{\theta}\,\)が含まれており、この式を量子の世界にもちこむには、その\(\,p_{\theta}\,\)を演算子に置き換えればよい。その結果

    \(\quad p_{\theta}=\sqrt{J(J+1)}\hbar\quad\cdots\,\)(2.7)

    と置き換えることが出来る。従って、式(2.3)の2原子分子の場合の量子論的な表現は式(1.19)のように

    \(\quad E=\dfrac{J(J+1)\hbar^2}{2I}=BJ(J+1)\qquad (J=0,1,2,\cdots)\quad\cdots\,\)(2.8)

    となる。ここで\(\,J\,\)は回転量子数であり、\(\,0\,\)を含む整数のみが許される。すなわち回転エネルギーは量子化されており、飛び飛びの準位のみが許されることになる。\(\,B\,\)は回転定数と呼ばれるもので次式で表される。

    \(\quad B=\dfrac{h}{8\pi^2I}\quad\)周波数単位}\(\qquad \dfrac{h}{8\pi^2c_0I}\quad\)波数単位\(\quad\cdots\,\)(2.9)

    2-2. 観測される回転スペクトル

    回転運動だけに起因するスペクトルはマイクロ波領域から遠赤外線領域に現れる。図は期待状態のCO分子の吸収スペクトルである。規則的な間隔の吸収線からなる。式(2.8)より2原子分子の場合の回転エネルギー準位は下図のようになることが理論的に予測できる。

    回転量子数の\(\,J=0,1,2,3,\cdots\,\)に対して、エネルギー準位は\(\, 0,2B,6B,12B,\cdots\,\)となる。従って隣接するエネルギー準位間の大きさは \(\,2B,4B,6B,\cdots\,\)である。これを横軸にエネルギーをとってプロットすると図の下に示したようにスペクトルの位置が\(\,2B\,\)の間隔で並ぶ。これは実測のスペクトルとぴったり一致する。

    このようにして、光子を1つ吸収することによる(2原子分子の)回転遷移は、\(\,J\,\)の値が1つだけ変化する状態の間に起こることがわかる。吸収はエネルギーが増える方向への遷移であるが、エネルギーが減る方向の瀬には発光スペクトルとして観測され、実際同様なスペクトルを示すことがわかっている。すなわち回転遷移の選択則は\(\,\Delta J=\pm 1\,\)である。ここでは、この選択則を実測との対応から導いたが、理論的にも同じ結果が導かれる。

    以上のようにして、実測スペクトルの間隔から回転定数\(\,B\,\)が決定できる。この\(\,B\,\)は(2.9)式に定義された量であるから、\(\,B\,\)がわかるということは慣性モーメント\(\,I\,\)が実験的に求められたことを意味する。これより\(\,B\to I\to r\,\)のようにして、分子の構造を与える情報\(\,r\,\)が得られることになる。

    2原子分子であるNO分子についても同様な解析が可能であり、実測スペクトルから分子構造に関する情報を導くことができる。それらの結果を以下の表にまとめる。

    \(\mu\) (kg) 回転定数\(\,B\) (Hz) \(B\,\)cm\(^{-1}\) 結合距離 \(r\,\mathring{\mathrm{A}}\)
    \(^{12}\)CO \(1.139\!\times\!10^{-26}\) \(5.793\times 10^{10}\) \(1.932\) \(1.128\)
    \(^{13}\)CO \(1.191\!\times\!10^{-26}\) \(5.538\times 10^{10}\) \(1.847\) \(1.128\)
    NO \(5.012\!\times\!10^{10}\) \(1.672\) \(1.151\)

    COの場合、\(J\!=\!1\!-\!0\,\)の回転遷移では、\(E=BJ(J+1)\,\) より、\(115.9\)GHz となり、\(^{13}\)COでは、\(110.8\)GHzとなる。\(J\!=\!2\!-\!1\,\)の場合は \(\,6B-2B=4B\,\)なので、\(^{12}\)CO は\(231.7\)GHzで、\(^{13}\)COは\(\,221.5\,\)GHzになる。

    2-3. 観測されたスペクトル


    爆発的星形成銀河(starburst銀河)M82の中心数100pc領域における、赤外線から電波に至る波長領域のスペクトル。大別して、(1)周波数に対して、連続的に放射強度が変化している「連続波放射」と、(2)ある特定の周波数でのみ、強い放射を示す「スペクトル線放射」とに分けられる。
    連続波放射は、主に波長の長い電波域におけるシンクロトロン放射(磁場の中を高速で運動している電子からの放射)と、サブミリ波~赤外線域にかけて、大きな山をつくっている成分、すなわち、星間ダスト(星間空間中に存在する、シリケイトなどの固体微粒子)からの熱放射の重ね合わせとしてよく説明でできる。
    この他、ミリ波付近には自由・自由遷移放射(電離領域中を運動する自由電子からの放射)も見られる。一方、この波長帯のスペクトル線放射としては、原子からの再結合線(桃色)・微細構造線(赤)・超微細構造線、分子からの回転線(青)などがみられる。







    ドメイン取得とSSL化

    Google Chrome を使っていると、http でのアクセスの場合、保護されていないページと表示されるので、SSL化をしようとトライしてみました。
    ドメインはデフォルトのままでは、無料SSL導入が出来ない(ロリポップでは)らしいので、しょうがないからドメインを導入しました。”tamami“という名前は、1970年に桜花賞をぶっちぎりで勝った「タマミ」からとっているので、これをそのまま使うことにしました。*.tokyo だと安かったんですけど、何となく “.tech“にしてしまいました。よく値段を確認しなかったんですけど、更新費用が年額 3,980円でサーバーのレンタル料(3,000円)より高い....
    結局、GMO としてこういう形でコストを回収しているんだっと言うことを、やっと勉強した。😭

    これで、ようやくSSL化が出来ました。内部のアドレス変換で RegEx というツールで変換したのが、失敗で \(\TeX\,\)の制御コードのバックラッシュ “\”が消えてしまいました。何か設定があるんだと思うけど、やってしまったことはしょうがないので、手動で復旧した。これで、 Google Chrome から文句は言われないかな。

    ネットの入試問題など-2

    問4

    \(\quad x^5=i\quad\)を解け。

    【解法】
    \(\quad i^5=i\quad\)なので、与式は\(\,\,x^5-i^5=0\quad\)と変形できる。

    \(\quad x^5-i^5=(x-i)(x^4+x^3i+x^2i^2+xi^3+i^4)\)
    \(\qquad =(x-i)(x^4+x^3i-x^2-xi+1)\quad\)ここで、\(x\ne0\quad\)なら

    \(\qquad =(x-i)\left(x^2+xi-1-\dfrac{1}{x}+\dfrac{1}{x^2}\right)=0\)

    ここで、\(\,x-\dfrac{1}{x}=t\,\,\)とすると、
    \(\quad t^2+2+it-1=t^2+it+1=0\quad\)より

    \(\quad\left(t+\dfrac{1}{2}i\right)^2+\dfrac{5}{4}=0\quad\)なので

    \(\quad t+\dfrac{1}{2}i=\pm\dfrac{\sqrt{5}}{2}i\quad\)よって\(\quad t=\dfrac{-1\pm\sqrt{5}i}{2}\)

    ここで\(\,\,x=\cos\theta+i\sin\theta\quad\)とおくと

    \(\quad t=\cos\theta+i\sin\theta-(\cos\theta+i\sin\theta)^{-1}\)

    \(\qquad =\cos\theta+i\sin\theta-\cos\theta+i\sin\theta=(2\sin\theta)i\)

    \(\quad \sin\theta=\dfrac{-1\pm\sqrt{5}}{4}\quad\)となる

    \(\quad \cos^2\theta=1-\sin^2\theta=\dfrac{10\pm2\sqrt{5}}{16}\quad\)より

    \(\quad \cos\theta=\pm\dfrac{\sqrt{10\pm2\sqrt{5}}}{4}\quad\)なので、解は \(\,i\,\)を頂点とする5角形の点。

    \(\quad x=i\,\,\,,\,\,\pm\dfrac{\sqrt{10+2\sqrt{5}}}{4}+\dfrac{-1+\sqrt{5}}{4}i\,\,,\)

    \(\qquad\qquad\quad \pm\dfrac{\sqrt{10-2\sqrt{5}}}{4}+\dfrac{-1-\sqrt{5}}{4}i\)

    また、\(\,i=\cos 5\theta+i\sin 5\theta=\cos\dfrac{\pi}{2}+i\sin\dfrac{\pi}{2}\quad\)なので

    \(\quad \theta=\dfrac{\pi}{10}\quad\)となり、\(\sin\dfrac{\pi}{10}=\dfrac{-1+\sqrt{5}}{4}\quad\)である。

    問5

    \(\quad \displaystyle\int_0^{\pi}e^x\sin^2x\,dx>8\quad\)を示せ

    【解答】
    \(\quad I=\displaystyle\int_0^{\pi}e^x\sin^2x\,dx\quad\)とする。

    \(\quad\sin^2x=\dfrac{1-\cos2x}{2}\quad\)なので、\(\,I=\dfrac{1}{2}\displaystyle\int_0^{\pi}e^x(1-\cos2x)dx\)

    \(\quad I=\dfrac{1}{2}\left\{\displaystyle\int_0^{\pi}e^x\,dx-\displaystyle\int_0^{\pi}e^x\cos2x\,dx\right\}\)

    ここで、\(\displaystyle\int_0^{\pi}e^x\cos2x\,dx=I_2\quad\)とおく。さきに微分を計算する。

    \(\quad e^x\sin2x=a\,\,,\quad e^x\cos2x=b\quad\)とおく

    \(\quad(e^x\cos2x)^{\prime}=e^x\cos2x-2e^x\sin2x=b-2a\quad\cdots\,\)(1)

    \(\quad(e^x\sin2x)^{\prime}=e^x\sin2x+2e^x\cos2x=a+2b\quad\cdots\,\)(2)

    \(\quad b=\dfrac{(1)\!+\!(2)\!\times\!2}{5}=\dfrac{1}{5}(e^x\cos2x+2e^x\sin2x)^{\prime}=e^x\cos2x\)

    \(\quad\therefore\,\,I_2=\dfrac{1}{5}\bigl[e^x\cos2x+2e^x\sin2x\bigr]_0^{\pi}=\dfrac{1}{5}(e^{\pi}-1)\)

    \(\quad I=\dfrac{1}{2}\left(e^{\pi}-1-\dfrac{1}{5}(e^{\pi}-1)\right)=\dfrac{2}{5}(e^{\pi}-1)\quad\)となるので

    \(\quad e^{\pi}>21\quad\)を示せばいいことになる。

    ここで\(\,\,y=e^x\quad\)のグラフで\(\quad x=3\quad\)での接線を考える。

    接線の式は\(\quad y=e^3(x-3)+e^3=e^3x-2e^3\quad\)となる。

    \(\quad e^{\pi}\quad\)は、接線の\(\,\,x=\pi\quad\)の点より上にあるので

    \(\quad e^{\pi}>e^3\pi-2e^3=e^3(\pi-2)>2.7^3\times1.1\fallingdotseq21.65\)

    よって、標題は証明された。

    問6

    \(\,x>0\,\)のとき\(\quad\dfrac{x^4+x^2+1}{x^3+x}\quad\)の最小値を求めよ

    【解答】
    与式の分母、分子を\(\quad x^2\quad\)で割ると

    \(\quad =\dfrac{x^2+1+\frac{1}{x^2}}{x+\frac{1}{x}}\quad\cdots\,\)(3) となる。

    ここで\(\quad x+\dfrac{1}{x}=t\quad\cdots\,\)(4) とおくと、式(3)は

    \(\quad \dfrac{t^2-1}{t}=t-\dfrac{1}{t}\quad\cdots\,\)(5) となる。

    ここで\(\quad x>0\quad\dfrac{1}{x}>0\quad\)なので、

    \(\quad t=x+\dfrac{1}{x}\ge 2\sqrt{x\dfrac{1}{x}}=0\quad\)なので、\(\,t\ge 2\)

    ここで、式(5)を\(\quad f_1(t)=t\quad f_2(t)=-\dfrac{1}{t}\quad\)と分けると

    \(\quad t\ge 2\quad\)で、\(f_1(t)\,\,,\,f_2(t)\quad\)はともに単調増加

    なので、最小値は\(\quad t=2\,\to\,x=1\quad\)最小値は\(\quad\dfrac{3}{2}\)