流体さんのえっち〜〜〜〜！（H定理）

※この記事は航空宇宙工学科アドベントカレンダー用記事です．

adventar.org

こんにちは，皆さんいかがお過ごしでしょうか？僕はそろそろ修論提出なのになんも結果出てなくて焦っています． みなさん見事にアイキャッチ画像に釣られましたね？残念！今回は流体の中でもH定理と呼ばれる法則の紹介っていう，真面目な話でした！

はじめに

H定理は身近にある水や空気などの密な流体では問題になりません． しかし，特に希薄な気体では保存則と並んで重要な第四の法則として知られています． これらの理由は何なのでしょうか？ H定理の導出と絡めながら答えを探り，Boltzmann方程式含めた気体分子運動論の本質は何か，古典物理的に考えてみたいと思います．

※やたら長くなってしまいました，僕が話したかった話は「H定理について考えてみる 」に書かれているので，そこから読んでみてもいいかもしれません．

複雑すぎるよBoltzmann方程式

H定理はBoltzmann方程式の重要な性質として知られているので，まずはBoltzmann方程式とは何ぞ？から始めましょう．ここ以降出てくる数式は全て総和規約に従います．
（なお2年ほど前に流体最強の方程式 - 空飛ぶセロ弾きという記事でBoltzmann方程式の導出を行ったので，細かい事が気になる人は参照してください．）

Boltzmann方程式（外力なし）は次式で表されます．

Boltzmann方程式

\begin{align}\label{boltzmann} \frac{\partial f}{\partial t}+\xi_i\frac{\partial f}{\partial X_i} &=J(f,f)\\ J(f,f)&=\frac{1}{m} \int_{\boldsymbol{\xi}_{1}}\int_{\boldsymbol{\alpha}}(f'f'_1-ff_1)B(|\alpha_jV_j|, V) d\Omega(\boldsymbol{\alpha}) d^3\boldsymbol{\xi}_1 \end{align}

...いや複雑すぎぃ！なのでざっくりと意味を紹介します． イメージさえ掴んでもらえれば十分です．
なお，各変数は末尾に掲載されているので，そちらを参照してください．

速度分布関数とは？

通常の流体力学では現実世界における位置（ x,y,z）に対して密度や温度などのマクロな量を使います． なぜなら十分密な流体な場合，流体の分子速度はMaxwell分布に従っており，変数としてそれらのスカラー量だけ決めれば良いからです． 一方，流れが希薄になると分子間の衝突は十分行われず，分子の速度は必ずしもMaxwell分布に従いません． そこで，分子の速度を分布関数として表す必要が出てきます．これを速度分布関数と呼びます．

f:id:kane-please:20191203190329p:plain
速度分布関数の意味合い，微小体積内の速度分布関数はその位置，速度を持った分子数を表す

左辺はなに？

Boltzmann方程式の左辺は物体微分（的なもの）です． すなわち，時刻 tにおける速度空間内の体積素片 dX d\xiは時刻 t+dtで dX' d\xi'に変化するが その内の分子数が時刻 tから時刻 t+dtにかけてどのように変化するかを表しています．（図を参照）
気持ちとしては，速度分布関数の時間変化を見たいけれども，分子には速度があるので速度分布関数の位置はどんどん変化していく． なのでNS方程式と同様に分子の移流も考えればええやん！って感じです．

f:id:kane-please:20191204005040p:plain
左辺のイメージ

右辺はなに？

右辺がBoltzmann方程式のキモです．右辺は分子間衝突による変化量を表しています．なお分子間衝突がない場合，この項は0になります． 特にBoltzmann方程式では仮定として分子間衝突は２体衝突のみで，３体衝突以上は起こらないと仮定しています． 分子間衝突による変化は以下の二つに分ける事ができます．

  • J_G：衝突による発生項
    時刻 tには dX_i d\xi_iになかったが，微小時間 dtのうちに衝突して速度が変化し，時刻 t+dtで dX' d\xi'内に入ってきたもの

  • J_L：衝突による消滅項
    時刻 tには dX_i d\xi_iにあったが，微小時間 dtのうちに衝突して速度が変化し，時刻 t+dtで dX' d\xi'から出たもの

すなわち \begin{equation} J(f,f) = J_G - J_L \end{equation} 分子間の衝突があると微小時間 dtでも速度は劇的に変化するので，対象の体積素片 dX d\xiの周囲だけでなく，速度空間全体に対してこれらの項は考える必要があります．

f:id:kane-please:20191204010550p:plain
衝突により外部から分子が飛び込んできたり（ J_G），内部から分子が飛び出ていく（ J_L）
ここから，さらに J_G， J_Lに関して詳しく見ていきます．

消滅項 J_Lの導出

さて，ここからは具体的に衝突の様子を考える必要があります． 消滅項を考える際には \xi_iの速度を持った分子０に \xi_{i1}の速度をもつ分子1が短い時間 dtの間に衝突します．これを絵にすると下のような図になります．

f:id:kane-please:20191205144306p:plain
衝突過程の様子
さて，ここで分子１は dtの時間で相対速度分の V_idtしか進むことができません．なので，分子１は必ず図の灰色の領域に入っていますが，分子１のような粒子はいくつ あるのでしょう？
結論からいうと期待値的に「１よりずっと小さい」です．つまり，ほとんどの場合０で時々１ということです． 理由は， \boldsymbol{X}というただでさえ小さい領域の中のさらに小さい領域 V_ie_idtdSだからです． なので，この領域内の分子数を数えるには以下のような統計的な仮説をおく必要があります．

  1. d\boldsymbol{X}内で状態は一様
  2. V_ie_idtdS内の分布も分子０の \xi_iによらず d\boldsymbol{X}と同じ fで与えられる（分子無秩序の仮定）
    特に今回は衝突前に対して分子無秩序が適用されているとします

この仮定を用いて， d\boldsymbol{X}の中で時刻 tから微小時間 dt内に速度 \xi_iをもつ分子の衝突数を求めます． 特に分子無秩序の仮定が適用できる，衝突をする直前で衝突数を考えます． これは次のように求められる．速度 \xi_{i1}と \xi_iが衝突する時，衝突数は
（上の灰色領域内の \xi_{i1}分子数）×（ d \boldsymbol{X}内の \xi_i分子数） で表されます．つまり， \begin{equation} \left( \frac{1}{m}f(X_i,\xi_{i1},t)d\boldsymbol{\xi_1} \times V_je_jdSdt\right)\left(\frac{1}{m}f(X_i,\xi_{i},t)d \boldsymbol{\xi}d \boldsymbol{X}\right) \end{equation}

ここで， \xi_iに衝突する分子の数を全て求めるには，

  • 衝突する分子の速度 \xi_{i1}に対して全積分
  • 衝突する微小面積 dSを全積分

すれば良い．以上から， d \boldsymbol{\xi}d \boldsymbol{X}内での衝突数 m^{-1}J_L d \boldsymbol{\xi}d \boldsymbol{X}dtが得られます．

\begin{align} J_L &= \frac{1}{m}\int_{全\xi_{i1},全S}f(\xi_{i1})f(\xi_i)(V_je_j)dS d \boldsymbol{\xi}_1 \\ &= \frac{1}{m}\int_{全\xi_{i1},全S}f(\xi_{i1})f(\xi_i)B(|\alpha_jV_j|, V) d\Omega d \boldsymbol{\xi}_1 \end{align}

ここで， (V_je_j)dSは衝突係数 Bと立体角 \Omegaに取り替えました．（よく見るBoltzmann方程式の形にするため，深い意味はない）

発生項 J_Gの導出

先ほどと同様に考えます． 今度は時刻 tに任意の分子速度 \bar{\xi_i},\bar{\xi_{i1}}が衝突して，時刻 t+dtに \xi_i, \xi_{i1}になる状況を考えます． 先ほどと同様にして， \bar{\xi_i}と \bar{\xi_{i1}}の衝突数は次のように考えられます．

\begin{equation} \left( \frac{1}{m}f(\bar{\xi_{i1}})d\boldsymbol{\bar{\xi_1}} \times \bar{V_j}e_jdSdt\right)\left(\frac{1}{m}f(\bar{\xi_{i}})d \bar{\boldsymbol{\xi}} d \boldsymbol{X}\right) \end{equation}

ただし，今度は \xi{i1}に関して全積分，のように単純ではありません． \bar{\xi_i}と \bar{\xi_{i1}}の衝突の後，速度が {\xi_i}と {\xi_{i1}} になる組み合わせの集合 D(\bar{\xi_i}, \bar{\xi_{i1}})に対して足し合わせます．なので， d \boldsymbol{\xi}d \boldsymbol{X}内での衝突数 m^{-1}J_G d \boldsymbol{\xi}d \boldsymbol{X}dtは次の式になります．

\begin{align} J_G d\boldsymbol{\xi} = \frac{1}{m}\int_{D}f(\bar{\xi_{i1}})f(\bar{\xi_i})B(|\alpha_j\bar{V_j}|, \bar{V}) d\Omega d \bar{\boldsymbol{\xi}_1}d\bar{\boldsymbol{\xi}} \end{align}

ここで，二体衝突の式を適用することにより，積分範囲を変更します．具体的には，以下のように積分変数の変換を行います． \begin{align} \xi_i &= \bar{\xi_i} + \alpha_i \alpha_j \bar{V_j}, & \xi_{i1} = \bar{\xi_{i1}} + \alpha_i \alpha_j \bar{V_j} \\ \leftrightarrow\ \bar{\xi_i} &= {\xi_i} + \alpha_i \alpha_j {V_j}, & \bar{\xi_{i1}} = {\xi_{i1}} + \alpha_i \alpha_j {V_j} \end{align} （ここからの式変形はただ煩雑なだけなので，割愛します）
結局再び \xi_{i1}に関する全積分で求めることができるようになります．

\begin{align} J_G &= \frac{1}{m}\int_{全\xi_{i1},全S}f(\bar{\xi_{i1}})f(\bar{\xi_i})B(|\alpha_jV_j|, V) d\Omega d \boldsymbol{\xi}_1 \end{align}

ゆえに結局衝突項は次のように表すことができます．これは冒頭のBoltzmann方程式と同様な意味となっています． \begin{align} J(f,f) &= J_G - J_L \\ &= \frac{1}{m}\int_{全\xi_{i1},全S}\left(f(\bar{\xi_{i1}})f(\bar{\xi_i}) - f({\xi_{i1}})f({\xi_i})\right)B(|\alpha_jV_j|, V) d\Omega d \boldsymbol{\xi}_1 \end{align}

まとめ

ボルツマン方程式左辺は速度分布関数の時間変化を表し，右辺は分子間衝突による影響を表す． ただ，仮定として

  • 分子の衝突は２体衝突のみ

  • 微小体積 d\boldsymbol{X}内では状態は一様

  • 衝突前の分子の分布に分子無秩序の仮定を適用

特に重要なので，分子無秩序に関してさらに言及すると，分子の衝突前に分子無秩序が適用されることの意味とは， 「衝突前には２分子は（速度分布関数 fに従うので）相関があるが，衝突後には（必ずしも fに従わないので）もはや相関は見られない」 ということを言っています．

H定理を証明する

Boltzmann方程式がどういうものか何となくでも理解していただけたでしょうか？ ここからはH定理の証明をしていきます． 証明に必要な関係式の導出は必要になり次第求めます．

H定理とは次の関係式を指します．

H定理

外力の働かない系を考える．次のような関数 H, H_iを考え，特に Hは局所H関数と呼ばれる． \begin{align} H(X_i,t) &= \int f \ln c_0^{-1} f d\boldsymbol{\xi}\\ H_i(X_i,t) &= \int \xi_i\ f \ln c_0^{-1} f d\boldsymbol{\xi} \end{align} この Hに対して以下の不等式が成立し，等号は fがMaxwell分布の時に限る． これをH定理と呼ぶ． \begin{equation} \frac{\partial H}{\partial t} + \frac{\partial H_i}{\partial X_i} \le 0 \end{equation}

以下，証明していきます．

Boltzmann方程式(\ref{boltzmann})の両辺に (1 + \ln c^{-1}_0 f)をかけて， \xi_iの全空間に渡り積分します． ここで， \begin{align} (1 + \ln c^{-1}_0f)\frac{\partial f}{\partial t}&=\frac{\partial}{\partial t}(f\ln c^{-1}_0f) \\ (1 + \ln c^{-1}_0f)\xi_i\frac{\partial f}{\partial t}&=\frac{\partial}{\partial X_i}(\xi_i f\ln c^{-1}_0f) \end{align} が成立することから左辺は \begin{equation}\label{hlhs} LHS = \frac{\partial H}{\partial t} + \frac{\partial H_i}{\partial X_i} \end{equation} と表されます．
右辺について，

\begin{equation} RHS = \frac{1}{m}\int (1 + \ln c^{-1}_0f)(f' f'_1 - f f_1)Bd \Omega d \boldsymbol{\xi}_1 d\boldsymbol{\xi} \end{equation}

うーん，このままではどうしようもないですね．そこで次の対称関係式を示しておきます．

対称関係式

\xi_iのある関数を \psi(\xi_i)として，以下の J_\psi(f,f)を定義しておきます． \begin{equation}\label{sym_formula} J_\psi(f,f) = \frac{1}{m}\int\psi(\xi_i)(f'f'_1 - ff_1)Bd \Omega d \boldsymbol{\xi}_1 d\boldsymbol{\xi} \end{equation} この積分に対して，次の二つの操作を考えます．

  1. (\xi_i, \xi_{i1}) \rightarrow (\xi_{i1}, \xi_{i})　と変数を書き換える（変数変換ではないことに注意）
    • この時，ただ見た目を変えただけなので，ヤコビアンを考える必要はない．また， (f'f'_1 - ff_1)も対称な式なので，見た目は変わらない

  2. (\xi_i, \xi_{i1}, \alpha_i) \rightarrow (\xi'_i, \xi'_{i1}, \alpha_i)　と変数変換を行った後に (\xi’_i, \xi’_{i1}) \rightarrow (\xi_{i}, \xi_{i1})と書き換える
    • 前者の変数変換のヤコビアンは1である．後者に関して， 書き換えにより (f'f'_1 - ff_1)は符号が逆になる．

以上の操作を順次式(\ref{sym_formula})に適用していきます．
式(\ref{sym_formula})に操作 a. を適用すると

\begin{equation}\label{sym_a} J_\psi(f,f) = \frac{1}{m}\int \psi(\xi_{i1})(f'f'_1 - ff_1)Bd \Omega d \boldsymbol{\xi}_1 d\boldsymbol{\xi} \end{equation}

式(\ref{sym_formula})に操作 b. を適用すると \begin{equation} J_\psi(f,f) = \frac{1}{m}\int \psi(\xi'_{i})(ff_1 - f'f'_1)Bd \Omega d \boldsymbol{\xi}_1 d\boldsymbol{\xi} \end{equation}

式(\ref{sym_a})に操作 b. を適用すると \begin{equation} J_\psi(f,f) = \frac{1}{m}\int \psi(\xi'_{i1})(ff_1 - f'f'_1)Bd \Omega d \boldsymbol{\xi}_1 d\boldsymbol{\xi} \end{equation}

ここで，以上の４式を全て足し合わせると次の対称関係式を得ることができます． \begin{equation}\label{sym} J_\psi(f,f) = \frac{1}{m}\int \frac{1}{4}(\psi + \psi_1 - \psi' - \psi'_1) (f'f'_1 - ff_1)Bd \Omega d \boldsymbol{\xi}_1 d\boldsymbol{\xi} \end{equation}

H定理の証明に戻る

さて，この対称関係式(\ref{sym})に \psi = 1 - \ln c^{-1}_0fを代入します． \begin{align} RHS &= \int \frac{1}{4} \left( \ln c^{-1}_0f + \ln c^{-1}_0f_1 - \ln c^{-1}_0f' - \ln c^{-1}_0f'_1 \right)(f'f'_1 - ff_1)Bd \Omega d \boldsymbol{\xi}_1 d\boldsymbol{\xi}\\ &= \frac{1}{4} \int (f'f'_1 - ff_1)\ln\left(\frac{ff_1}{f'f'_1} \right)Bd \Omega d \boldsymbol{\xi}_1 d\boldsymbol{\xi} \\ & \le 0 \label{hrhs} \end{align} 最後の式では， (x-y)\ln(x/y) \ge 0を使用しました．

以上の結果から，式(\ref{hlhs})(\ref{hrhs})を用いて，

\begin{equation} \frac{\partial H}{\partial t} + \frac{\partial H_i}{\partial X_i} \le 0 \end{equation} が示されます．

はぁ〜〜，ようやくこれでH定理証明完了です．

H定理について考えてみる

ようやく本題です．今まで長々と話してしまってすみません...
ここからは，H定理に関して考察することで色々面白いことがわかるので，紹介します．

そもそもHって何なの？

H定理とかいうものを証明しましたが，これが何の意味を持っているのでしょう？ 実際に Hを計算してみることで，わかります．

fに次のようにMaxwell分布を代入します

\begin{equation} f = \rho (2\pi RT)^{\frac{-3}{2}}\exp\left(-\frac{(\xi_i - v_i)^2}{2RT} \right) \end{equation}

よって Hは次のように計算できます．

\begin{align} H &= \int f \ln c^{-1}_0 f d \boldsymbol{\xi} \\ &= \int \rho (2\pi RT)^{\frac{-3}{2}}\exp\left(-\frac{(\xi_i - v_i)^2}{2RT} \right) \left(\ln c^{-1}_0 + \left(-\frac{(\xi_i - v_i)^2}{2RT} \right) + \ln \left(\rho (2\pi RT)^{\frac{-3}{2}} \right) \right)d \boldsymbol{\xi}\\ &= \rho(\ln c_0^{-1}) + \rho\ln \left(\rho (2\pi RT)^{\frac{-3}{2}} ) \right) + \int \rho (2\pi RT)^{\frac{-3}{2}}\left(-\frac{(\xi_i - v_i)^2}{2RT} \right) \exp\left(-\frac{(\xi_i - v_i)^2}{2RT} \right) d \boldsymbol{\xi} \label{h_calc_1} \end{align}

なお，ガウス積分を計算すると，

\begin{equation} \int \rho (2\pi RT)^{\frac{-3}{2}}\exp\left(-\frac{(\xi_i - v_i)^2}{2RT} \right)d \boldsymbol{\xi} = \rho \end{equation}

ここで，式(\ref{h_calc_1})の残った積分を計算するために，次のガウス積分の公式を使用します． \begin{equation} \int^{\infty}_{-\infty}x^{2n}\exp({-ax^2})dx = \frac{(2n-1)!!}{2^n a^n}\sqrt{\frac{\pi}{a}} \end{equation}

これを使用するために，式(\ref{h_calc_1})で x = (\xi_i - v_i)/\sqrt{2RT}と変数変換します． ここで注意すべきは， i=0,1,2なので変数変換の際にかける \sqrt{2RT}は三回かけられるます． すると，結局以下のように表せます．

\begin{equation} \frac{H}{\rho} = \ln \left(\rho (2\pi RT)^{\frac{-3}{2}} ) \right) + const. \end{equation}

ここで，単原子分子気体のエントロピー sは以下のように表せます． \begin{equation} s = R \ln \left(\rho^{-1} (2\pi RT)^{\frac{3}{2}} \right) \end{equation}

よって，以下の関係式を得られます．

\begin{equation} \frac{H}{\rho} = -R^{-1}s + const. \end{equation}

これでやっとH関数の正体もといH定理の意味がわかりました．

H定理の意味

H関数はエントロピーに対応しており， H定理はエントロピー非減少を意味していた！

この結果は次のように言い換えることができます．エントロピー非減少というのは，時間的に戻らない，つまり時間的方向性を持っていることになるので，

H定理はBoltzmann方程式の解が時間的方向性を持っていることを示している．

あっれれ〜〜？おっかしいぞ〜〜？

さてさて，ここで話を終わっても良いのですが，もう少し考えたいと思います．
皆さんに思い出して頂きたいのですが，Boltzmann方程式の導出では二体衝突などの古典物理しか用いていませんでした． そうなると，先ほど得られたH定理の結果はなんだか気持ちが悪いような気がします． それに気が付いたのがヨハン・ロシュミットです．ロシュミットはH定理に対して次のような反論を行いました．

ロシュミットの逆行性批判

可逆的な古典物理の議論しか用いていないのに，Boltzmann方程式からH定理のような非可逆過程が導かれるのは，どこか間違っているのでは？

言われてみればその通りではないでしょうか？Boltzmann方程式の導出で用いた二体衝突は，当然ですが，逆に衝突後の速度で衝突すれば衝突前の速度が得られます． つまり全く逆のことを行えば以前の状態を得られます．このような可逆性は古典物理特有のものと言えます． 一方で，Boltzmann方程式から導かれるH定理はエントロピー増大の法則を示しており，一度エントロピーが増大すると元には戻らない，不可逆な過程になっていることを示しています．

可逆性を仮定しながら導いてきたのに，なぜか不可逆な結果が得られた，なんだか矛盾しているような気がします．

この答えは，Boltzmann方程式に課した仮定である分子無秩序にあります． ここで次の例を考えてみましょう．

  • 衝突後に分子無秩序の仮定を適用する
衝突後に分子無秩序の仮定を適用

詳細を述べるのは煩雑さを招くだけなので，割愛しながら結果を紹介します． 衝突の生成，消滅項 J_G, J_Lの計算においては，衝突前の状態に対して衝突数を計算していました． これは衝突前の状態に対して分子無秩序の仮定を適用したからです．

では，逆に衝突後の状態に対して分子無秩序の仮定を適用したらどうなるのでしょうか？ その時の衝突項を紛らわしさを防ぐために　 \hat{J}　と変えておきます． それを考えるために，今回は衝突後の衝突数を考えます． ゆえに，下図のような灰色領域を今回は考えることになります．

f:id:kane-please:20191206145042p:plain
衝突過程の様子

考えてみればわかりますが，今回は衝突後の分布に対してのみ fを適用できるので，計算方法が変わってきます．

  • 消滅項の場合
    衝突後の速度が \bar{\xi_i}, \bar{\xi_{i1}}であるとして，衝突前が \xi_i, \xi_{i1}となる組み合わせに対して衝突数を計算． 以前の生成項と同様な計算になり， \begin{align} \hat{J}_L &= \frac{1}{m}\int_{全\xi_{i1},全S}f(\bar{\xi_{i1}})f(\bar{\xi_i})B(|\alpha_jV_j|, V) d\Omega d \boldsymbol{\xi}_1 \end{align}

  • 生成項の場合
    衝突後の速度が \xi_i, \xi_{i1}として，以前の消滅項と同様に衝突数を計算すると， \begin{align} \hat{J}_G &= \frac{1}{m}\int_{全\xi_{i1},全S}f({\xi_{i1}})f({\xi_i})B(|\alpha_jV_j|, V) d\Omega d \boldsymbol{\xi}_1 \end{align}

よって，今回の衝突項 \hat{J}は次のように計算されます． \begin{align} \hat{J} &= \hat{J}_G - \hat{J}_L \\ &= -J \end{align}

あれ，逆になってしまいました． ということは，この仮定のもとではH定理も逆の結果となり，「エントロピー非増大の法則」が得られてしまいますw

鋭い人はもうお気づきかと思いますが，本来可逆に思われるBoltzmann方程式に非可逆性をもたらしているのは，この分子無秩序の仮定です． そして，さらに重要なのが，

  • 分子衝突前に分子無秩序の仮定をし

  • 分子衝突後に適用はしない

というところです．つまり，こう言い換えることもできます．

Boltzmann方程式の解の時間的方向性はそもそも分子無秩序の仮定の置き方に起因する

結論

大変長くなってしまい申し訳ありません． 最後に僕が伝えたかったことをまとめておきたいと思います．

まとめ

  • H定理はエントロピー非減少の法則を示し，Boltzmann方程式の解は不可逆，つまり時間的方向性があることがわかる．
  • 一見，可逆に見えるBoltzmann方程式に時間的方向性をもたらしたのは，元を辿ると分子無秩序の仮定である．

今回の話から得られる教訓として，一見「分子無秩序の仮定」のようにとても妥当に見える仮定があったとしても，実は方程式の重大な性質を 決めてしまうほど強力な威力をもつことがあるので注意！ということでしょうか．
あと，個人的にエントロピー非減少の法則の証明って世の中にあんまりないかなーって思って今回示すことができてちょっと嬉しいです．

自分本当に文章力ないので，めちゃわかりにくくなってしまって申し訳ないです．．． ただ，自分以外の航空宇宙民は毎回本当に面白い文章書いてるし，ハイスペなので，ぜひ読んでみてください！

参考文献

  • 曾根良夫・青木一生（1994）"分子気体力学"

  • Carlo Cercignani（1975）"The Boltzmann Equation and Its Applications"

変数名一覧

変数 変数名
f 速度分布関数
f ' 　衝突後のf
t 時間
\xi 分子速度
X （速度分布関数の）位置
B(|\alpha_j V_j|, V) 衝突係数
\Omega 衝突断面積
J(f,f) 衝突項
J_G 衝突による発生項
J_L 衝突による消滅項
V_i 分子０と分子１の相対速度

f:id:kane-please:20191206160903j:plain:w50

【Python】 np.arrayに関する笑った仕様

クイズ

test1()〜test5()はそれぞれ実行可能でしょうか？それとも実行不可能でしょうか？

メインの話はtest3() test4() test5()ですけど，ついでなので，リストも絡めて見ました． Python書けると豪語される方ならきっと全問正解ですよね！

import numpy as np

def test1():
    a = [[10.0*j+k for k in range(5)]for j in range(5)]
    a /= 10
    print(a)

def test2():
    a = [[10*j+k for k in range(5)]for j in range(5)]
    a /= 10
    print(a)

def test3():
    a = [[10.0*j+k for k in range(5)]for j in range(5)]
    a = np.array(a)
    a /= 10
    print(a)

def test4():
    a = [[10*j+k for k in range(5)]for j in range(5)]
    a = np.array(a)
    a /= 10
    print(a)

def test5():
    a = [[10*j+k for k in range(5)]for j in range(5)]
    a = np.array(a)
    a = a/10
    print(a)

def main():
    test1()
    #test2()
    #test3()
    #test4()
    #test5()
    
if __name__=='__main__':
    main()

【NLP】【Python】CaboChaのPythonバイディングで詰まった

CaboChaをPythonで使いたい

今のバイト先でNLPのためのコードを主に書いています． そこで確か前にMeCabとCaboChaの二つをバインディングした仮想環境を作った気がするんですが，最近みたらなぜかその仮想環境が壊れて使えなくなってました． （まぁ，LLVM入れたり，CRF++入れ直したりしてたので．．．） そこで，再び両方とも入った仮想環境を作ろうとしたら詰まったので，同じ苦しみを味わないためにメモ．

環境

エラーコード

$ python setup.py install

warning: no library file corresponding to '-L/usr/local/Cellar/mecab/0.996/lib' found (skipping)  
clang++ -bundle -undefined dynamic_lookup -isysroot /Applications/Xcode.app/Contents/Developer/Platforms/MacOSX.platform/Developer/SDKs/MacOSX10.14.sdk build/temp.macosx-10.14-x86_64-2.7/CaboCha_wrap.o -L/usr/local/Cellar/cabocha/0.69/lib -lcabocha -lcrfpp -lmecab -liconv -lmecab -lstdc++ -o build/lib.macosx-10.14-x86_64-2.7/_CaboCha.so  
ld: library not found for -lcrfpp  
clang: error: linker command failed with exit code 1 (use -v to see invocation)  
error: command 'clang++' failed with exit status 1

解決策

この記事に書いてあることをコピペしました． www.maytry.net どうやらエンコードの違い？みたいですね． 正直エラーコードからそんなこと1mmもわからなかったので，本当に詰んでました．．．

まずは公式サイトから自分の持っているCaboChaに対応するバージョンの.tar.bz2ファイルをダウンロードしてください． そのあとそれを解凍して，そのディレクトリに入ります．あとは 以下のコードを実行すればおkです．

$ cd cabocha-0.69
$ ./configure --with-mecab-config=`which mecab-config` --with-charset=utf8 # これでエンコードのエラーが出ないはず
$ make
$ sudo make install
$ cd python
$ sudo python setup.py install

ちなみに，自分の場合この前にこのサイトを参考にして以下のものを試しました．

$ ./configure --with-mecab-config=`which mecab-config` --with-charset=UTF8 --with-posset=UNIDIC
$ make
$ make install
$ cd python
$ sudo pyhton setup.py install

この方法は試しましたが全く効果がありませんでした． 自分のはPython2でこのサイトはPython3だったので，そのために上手くいかなかったのかも．

正規表現チートシート

正規表現の参照リスト

正規表現のチートシート的なのが欲しくてメモ． 以下のサイトから完全にパクりました．自分用のメモが欲しくてここにもメモさせていただきました． qiita.com

文字 説明 同様 例 マッチする マッチしない
\d 任意の数字 [0-9]
\D 任意の数字以外 [^0-9]
\s 任意の空白文字 [\t\n\r\f\v]
\S 任意の空白文字以外 [^\t\n\r\f\v]
\w 任意の英数字 [a-zA-Z0-9_]
\W 任意の英数字以外 [\a-zA-Z0-9_]
\A 文字列の先頭 ^
\Z 文字列の末尾 $
. 任意の一文字 - a.c abc, acc, aac abbc, accc
^ 文字列の先頭 - ^abc abcdef defabc
$ 文字列の末尾 - abc$ defabc abcdef
* ０回以上の繰り返し - ab* a, ab, abb, abbb aa, bb
+ １回以上の繰り返し - ab+ ab, abb, abbb a, aa, bb
? ０回または１回 - ab? a, ab abb
{m} m回の繰り返し - a{3} aaa a, aa, aaaa
{m,n} m〜n回の繰り返し - a{2, 4} aa, aaa, aaaa a, aaaaa
[] 集合 - [a-c] a, b, c d, e, f
縦線 和集合（または） - a縦線b a, b c, d
() グループ化 - (abc)+ abc, abcabc a, ab, abcd

【Octave】Octaveの環境構築（MacOS Mojave）

本当はMATLABが使いたかった

最近読んでる論文がソースコードを公開しているので、お試しで一度動かしてみたいなぁと思って開いてみるとまさかの.mファイル...
調べてみたらMATLABやんけ！！ってなって、学校のPCでもいいけど、いじるなら手元に環境欲しいので仕方なくOctaveを入れることを決意しました。

結論

Octaveのwikiがあったので、これを読んだら環境を作る事ができました。

Octave for macOS - Octave

これ読むのも面倒って人は脳死で次のコマンドを打ちまくりましょう。

$ brew update
$ brew install qt
$ brew install octave --with-qt

この後、Qtの代わりにgnuplotを起動させる設定とかがあるけど、グラフは普通にかけたので、今回は飛ばしました。

上記のwikiの最後のCreate a launcher app with AppleScriptと言うところでApple Scriptを書く事でGUIのOctaveを作る事ができそう？です。

失敗例

最初にこの記事を参考にしてしまいました。

MacにOctaveをインストール : 徒然なるままに

この通りにOpenBLASをtapして入れたりしたのですが、最後のインストールのところでどうしてもエラーが出てしまいました。 エラーの内容も他のファイルに書き込まれている内容が予想と違うみたいなものだったので、今までの何かと衝突してる可能性があります。

どちらにしろ、やっぱり頼るべきは公式ドキュメントですね！

【流体力学】【C++】圧縮性CFDの基礎

流体力学 プログラミング

※この記事は「東京大学航空宇宙工学科／専攻 Advent Calendar 2018」用記事です．

-- 2020/11/2 修正 --
FVSのコードミスが発覚したため、結果を更新致しました。@ur_kinsk様、ご指摘頂きありがとうございます！

はじめに

自分は流体系の人間なので、今回も何か流体、特に数値流体の話がしたいなぁと思って考えてました。 すると「あ、そういえば、空気力学第四でCFDやったけど、あれムズすぎてわからなかった気がする...　」 ってなったので、今回は空気力学第四のパクリで圧縮性の話をしたいと思います。

Shock Captureの技術

圧縮性流体を考えるCFDの世界ではどのように衝撃波を捉えるかが重要になります。 CFD界隈では銀本と呼ばれる本（流体力学の数値計算法）を見ると、衝撃波捕獲法として重要な技法として二つあげています。

  • 流束差分法（FDS）
  • 流束分離法（FVS）

今回はその両方を紹介して、現在のCFD技術をちょこっと理解した気になりましょう！

一次元システム方程式

流体を少しでも齧ったことのある方ならきっとオイラー方程式くらいは知ってるはず。
今回の問題設定で考えるのは、そのオイラー方程式が使えるような状況、かつ一次元の衝撃波菅問題に関して考えたいと思います。 流体の変数というのは基本的に以下の三つの保存則から定まります。

  • 質量保存
  • 運動量保存
  • エネルギー保存

当たり前といえば当たり前ですが、実はこの三つを完全に満たすようにCFDを作るのは基本的にはできないはずで、今回は触れないですが、精度とこれらの保存則は割とトレードオフの関係にあります。 さて、ここで今回の問題で以上の保存則を表す方程式は以下の三つになります。

\begin{align} &\frac{\partial \rho}{\partial t} + \frac{\partial \rho u}{\partial x} = 0 \\ &\frac{\partial \rho u}{\partial t} + \frac{\partial (p + \rho^2 )}{\partial x} = 0 \\ &\frac{\partial e}{\partial t} + \frac{\partial(e+p)}{\partial x} = 0 \end{align}

ここで

\begin{equation} e = \frac{1}{2}\boldsymbol{v}^2 + e_{heat} \end{equation}

であるとして、 eは流体のもつ全エネルギー、 e_{heat}は流体の熱エネルギーであることに注意しましょう。 この式、普段古典的な流体力学の授業で紹介されるような式の形をしていません。（例えば、参考 http://www.se.fukuoka-u.ac.jp/iwayama/teach/gfd/2007/chap3.pdf ）
もちろん、粘性がないからとか、 eで全エネルギーを括ったからとかいうのはありますが、何より重要なのは、全て「時間微分項」＋「空間微分項」の形で表記されている点です。 ゆえに、上の式は以下のように書き換えることが可能です。 \begin{equation} \label{qt} Q_t + E_x = 0 \end{equation} \begin{equation} \boldsymbol{Q} = \left( \begin{array}{c} \rho \\ \rho u \\ e \end{array} \right),\ \boldsymbol{E} = \left( \begin{array}{c} \rho u \\ p + \rho u^2 \\ (e+p)u \end{array} \right) \end{equation} つまり、これが「保存形」と呼ばれる流体のシステム方程式の形になります。 確かにこれ、保存則の形そのまんまですしね。この式の上で、 Qは「保存量」 Eは「流束」と呼ばれます。

ここで、式(\ref{qt})の形のままでは何もできません。すなわち、この式を局所的に線形化します（係数マトリックスを凍結するといいます）。

\begin{equation} \label{freeze} \frac{\partial Q}{\partial t} + A \frac{\partial Q}{\partial x}=0 \end{equation}

ただし、係数マトリックス A\equiv \partial E/ \partial Qは A_{ij} = \partial E_{ij} / \partial Q_{ij}の ように計算できて、

\begin{equation} A = \left( \begin{array}{ccc} 0 & 1 & 0 \\ -\cfrac{3-\gamma}{2} u^2 & (3-\gamma u) & \gamma - 1 \\ \left(\cfrac{\gamma - 1}{2}u^2 - H \right)u & H - (\gamma - 1)u^2 & \gamma u \end{array} \right) \end{equation}

この Aを固有値分解すると次の三つの固有値を得ることができます。

\begin{align} \lambda_1 = u-c,&& \lambda_2 = u, && \lambda_3 = u+c \end{align}

ところで、式(\ref{freeze})を眺めると、これって輸送方程式 u_t + cu_x = 0と同じ形であることがわかります。つまり、 Aが輸送速度に対応していると予測が立ちます。 一般に Aは非対角成分を持ちますが、対角化することで独立な三つの波動方程式に分解でき、その三つの波の特性速度がこの行列の固有値として現れているのだと考えることができます。

実際、対角化は以下のように行われます。（後で使います） \begin{equation} A=R\Lambda R^{-1} \end{equation}

\begin{align} \boldsymbol{R} = \left( \begin{array}{ccc} 1 & 1 & 1 \\ u-c & u & u+c \\ H-uc& \cfrac{1}{2}u^2 & H+uc \end{array} \right),&& \boldsymbol{R}^{-1} = \left( \begin{array}{ccc} \cfrac{1}{2}\left(b_1 + \cfrac{u}{c}\right) & -\cfrac{1}{2}\left(\cfrac{1}{c} + b_2 u\right) & \cfrac{1}{2}b_2 \\ 1-b_1 & b_2u & -b_2 \\ \cfrac{1}{2}\left(b_1 - \cfrac{u}{c}\right) & \cfrac{1}{2}\left(\cfrac{1}{c} + b_2 u\right) & \cfrac{1}{2}b_2 \end{array} \right) \end{align}

なお、 \begin{align} b_1 = \frac{u^2}{2}\frac{\gamma - 1}{c^2}, && b_2 = \cfrac{\gamma - 1}{c^2} \end{align}

近似リーマン解法

CFDの基本的な手法として風上差分法というものが存在するのは知ってるでしょうか？ 風上差分法というのは、 \begin{equation} \frac{\partial u}{\partial x} = \frac{u_{n+1}-u_{n}}{\Delta x} \end{equation} のように移流項を風上側で差分をとって評価するもので 、圧縮性の流体を考えるときには 特に頻繁に使われます。なぜ中心差分法でも、後退差分法でもないのか？今回は詳しく説明しないですが、 流れのもつ「情報」というものは 基本的に「超音速」を超えた場合に風上からしか来ないからです。 なぜなら、流れの情報というものは基本的に「音速」で伝達するものであり、流れが超音速を超えると、後ろからの情報は音速でしか伝わらないので、 上流に伝播しなくなるからです。

さて、このような背景から今回は風上差分法を使います。 ここまでをまとめるとこんな感じです。

  • オイラー方程式には三つの波が存在
  • 情報は風上側から流れてくるので、風上差分法を用いるのが適切

さて、今までの話から、システム方程式の解法では三つの波が独立に存在し、それぞれに対して高解像度差分法を適用することが可能になります。 CFDの内部では次のグラフのような形でデータが存在しているとイメージしてください。

f:id:kane-please:20181205143345p:plain
CFDのイメージ
この次の時間ステップを求めるためには各セル境界面での流束を評価する必要があります。 このように、二つの領域に渡って不連続な初期値をもつ双曲型の連立方程式を「リーマン問題」と呼びます。
f:id:kane-please:20181205151616p:plain
リーマン問題
このリーマン問題の解から、セル境界面における通過する流束を求めることができる。 この流束を \tilde{E}_{j+1/2}とすると、

\begin{equation} \label{int} Q_j^{n+1}=Q_j^n - \frac{\Delta t}{\Delta x}(\tilde{E}_{j+1/2}^n - \tilde{E}_{j-1/2}^n) \end{equation}

として、次の時刻での保存量 Q_j^{n+1}を求めることができる。この式から、各セルにおける物理量分布を求めることが可能になります。 ところで、そもそもCFDの内部ではセル内で物理量が一定分布の近似したデータについて考えているので、その時点ですでに精度は失われており、 わざわざ計算コストの高いリーマン問題を厳密に解く必要はありません。 そこで、近似的にリーマン問題を満足するように \tilde{E}_{j+1/2}^nを決めてしまえばいい。 その近似手法の中でもよく知られた（古典的な）ものは以下の二つになります。

  • 流束差分法（FDS）
  • 流束分離法（FVS）

以下ではそれぞれの手法の説明をした後に、結果を比較したいと思います。

FDS (Flux Difference Splitting)

いわゆる「Roeスキーム」と呼ばれる手法になります。 この手法の特徴としては、離散化されたデータからRoeの近似リーマン解を用いて出来るだけ物理的に正しい j+1/2での 数値流束を評価することで風上差分を行う、ということです。 具体的に見てもらった方がわかりやすいかもしれないので、早速説明していきます。

Roeは非線形スカラー方程式に対する１次の風上差分を非線形のシステム方程式に拡張しました。 このスキームでは１次精度の数値流束は以下の式で評価されます。

\begin{equation} \label{ffds} \tilde{E}_{j+1/2} = \frac{1}{2}(E_{j+1} + E_{j} -|A|_{j+1/2}(Q_{j+1} - Q_j)) \end{equation}

ただし、 \begin{equation} |A|_{j+1/2} = R_{j+1/2}|\Lambda|_{j+1/2}R_{j+1/2}^{-1} \end{equation}

ここで、添字の j+1/2で表される量は以下のRoeの平均を用いて評価されます。 また、ここでの絶対値 |\cdot|の意味は全ての行列要素で絶対値をとる、という意味です。 このRoeの平均というものはRoeが提案した線形分解によるリーマン解法 であり、平均の値 Q_{ave}が左右の状態 Q_L, Q_Rの非線形関数として定義されています。 \begin{align} \rho_{ave} &= \sqrt{\rho_L \rho_R}\\ u_{ave} &= \frac{\sqrt{\rho_L}u_L + \sqrt{\rho_R}u_R}{\sqrt{\rho_L} + \sqrt{\rho_R}}\\ H_{ave} &= \frac{\sqrt{\rho_L}H_L + \sqrt{\rho_R}H_R}{\sqrt{\rho_L} + \sqrt{\rho_R}}\\ c_{ave}^2 &= (\gamma - 1)\left(H_{ave} - \frac{1}{2}u_{ave}^2\right) \end{align}

FDSでの解析はまず式(\ref{ffds})により流束 \tilde{E}_{j+1/2}を求めた後に式(\ref{int})に基づいて時間積分を行えば良いです。

式(\ref{ffds})の最後の項は数値粘性を表しており、証明は省略しますが、この項は物理量の「跳び」（例えば衝撃波付近は大きい）に対して 数値粘性を付加しています。逆にいうと、「跳び」の少ない衝撃波から離れた領域に対しては数値粘性を付加しないので、正確になります。 この点がFDSの優れている点です。

FVS (Flux Vector Splitting)

こちらは実際にはリーマン問題の近似解を実際に求めている訳ではないのですが、「スカラー方程式への風上差分をシステム方程式にも拡張している」 という意味でなぜか近似リーマン解法に分類されています。 この手法の特徴としては、「物理的な流束を各特性波の方向に応じて差分方向を切り替える」という点です。

考え方はとても簡単です。まず 流束ベクトルを単純に E^{+}(Q),E^{-}(Q)に分解します。ここで、 E^+は左からくる正の固有値に対応する流束、 E^-は右からくる負の固有値に対応する流束です。正の固有値を持つということは波（情報）が左からくるということなので後退差分を利用します。 逆に負の固有値を持つならば波（情報）が右からくるということなので前進差分を利用します。それだけです。

それでは具体的に手法を見ましょう。 流速の分解は実際には以下の式で行います。 \begin{equation} \Lambda^\pm = \frac{\Lambda \pm |\Lambda|}{2} \end{equation} \Lambda^+は正の固有値を持つもののみから、 \Lambda^-は負の固有値を持つもののみからなる。 これより新しく以下のように A^\pmが定義できます。 \begin{equation} A^\pm = R\Lambda^\pm R^{-1} \end{equation} さてこれらを用いて数値流束 \tilde{E}_{j+1/2}を表示すると、 \begin{equation} \tilde{E}_{j+1/2} = E_{j+1/2}^+ (Q_{j+1/2}) + E_{j-1/2}^-(Q_{j+1/2}) \end{equation} 今回は一次精度の流束分離を行い、その時には次のように置くことができる。 \begin{align} Q_{j+1/2,L} = Q_j, && Q_{j+1/2,R} = Q_{j+1} \end{align} ここで、 L, Rはそれぞれ左からくる波、右からくる波に対応しています。つまり E^+については、 Q_Lを、 E^-については、 Q_Rをそれぞれ使用する必要がある。

よって数値流束 \tilde{E}_{j+1/2}は次のように評価されます。 \begin{align} \tilde{E}_{j+1/2} &= E_{j+1/2}^+ + E_{j-1/2}^-\\ \label{ffvs} &= A_j^{+}Q_j + A_{j+1}^-Q_{j+1}\\ &= \frac{A_j + |A|_j}{2}Q_j + \frac{A_{j+1} + |A|_{j+1}}{2}Q_{j+1}\\ \label{flux_fvs} &= \frac{1}{2}(E_{j+1} + E_{j} -(|A|_{j+1}Q_{j+1} - |A|_{j}Q_j)) \end{align} あとはFDSの時と同様にして、式(\ref{ffvs})により流束 \tilde{E}_{j+1/2}を求めた後に式(\ref{int})に基づいて時間積分を行うだけです。

さて、ここでFDSとFVSの違いに関して軽く述べると、式(\ref{ffds})と(\ref{flux_fvs})を見比べると、最後の項だけ違うことがわかります。
そもそも対流項 \partial E/ \partial xの計算をしようとした時に、 E_{j-1/2} - E_{j+1/2}をすることを考えると、最初の二項は 二次精度中心差分を形成することがわかります。そして、残った最後の項が拡散項として働き、このスキームに風上差分のような特徴を加えていると考えることができます。 （一般にCFDは中心二次精度差分は正確だが不安定、風上一時精度差分は安定だが鈍るのが特徴で、その両者の組み合わせをするのが差分スキームの基本です）
この拡散項がFDSでは必ず正になり、拡散を引き起こしますが、一方FVSではここが必ずしも正にはならないので、そこが本質的に異なっていると考えることができます。

実装

ここに示すと醜くなるので、Githubに置いておきました。見てみたい方はみてみるといいかもしれません。
ただ、適当にパパッと作ったので、クソコードなのは許してください。 あと、Eigen使ってるので、そこだけ注意お願いします。

結果

さて、それではFDSとFVSが実際の問題に適用した時にどのような挙動になるのか見てみましょう。 今回は問題として、一次元の衝撃波管問題を考えました。 CFDの世界では衝撃波管 問題が頻繁に用いられるのですが、その理由は

  • 解析解がわかっているから
  • 衝撃波、膨張波、接触面という圧縮性衝撃波で重要となる要素が全て含まれているから として、頻繁にスキームのテスト問題として用いられます。

-- 2020/11/2 修正 ここから --
以下の三つのグラフが解析解、FDS、 FVSの密度、速度、圧力の結果になります。

f:id:kane-please:20201102163849p:plain
一次元衝撃波管問題の速度グラフ

f:id:kane-please:20201102164034p:plain
一次元衝撃波管問題の圧力グラフ

f:id:kane-please:20201102164048p:plain
一次元衝撃波管問題の圧力グラフ

さて、FDSはやはり特徴として衝撃波付近で鈍っていますね。これはそもそもスキームの段階で数値粘性がそのように働くよう作ったので狙い通りです。 それでFVSはっと... うーん...FDSはまだそんなものかと思えますが、FVSの結果が異常に悪いですね...もしかしたら実装を間違えたかもしれないので、そこだけ申し訳ありません。

コードミスを修正した結果、FVSの方もFDSとかなり近い結果が得られています。やはり初期のFDS, FVSでは実用上では大きな違いがなかったように見られます。

一応、FVSは論文に衝撃波管問題への適用結果があったので、載せておきます。

f:id:kane-please:20181205191714p:plain
FVMの衝撃波管問題への適用結果（paperより）
いや、これも相当ひどいな（笑） 今回私たちが学んだのは本当に最初のFVSです。このあと、FVSというジャンルの中には様々に新しいスキームが誕生しました。新しいスキームが出るってことは、この 最初のFVSは使い物にならなかったからで、もしかしたら、自分の実装は間違ってなかったかもしれません（笑）

-- 2020/11/2 修正 ここまで --

終わりに

どうですか？意外と圧縮性のCFDってそんなに難しくないんですよね。 考え方はとても簡単なので皆さんに理解してもらえたなら幸いです。 ところで、FDSとFVS、世界では今どちらが最も使われていると思いますか？ ...答えはFVSなんです！なんででしょうね？今圧縮性で頻繁に用いられているのはFVSから派生したAUSM系スキームと呼ばれるものです。（参考に入れておきます） 最初はあんなに結果が悪くても、のちに誰かが発展させて広まっていくのは素晴らしいなぁと思いました。

誰か、自分の研究を代わりに発展させてくれないかなぁと思った今日この頃でした。

参考文献

【メモ】【Python】リスト内包表記っぽいgeneratorに関して

generatorは次の形でも表現できる

generatorは大抵defとyieldの組み合わせで表現される。 しかし、それ以外にもgeneratorはリスト内包表記的に生み出すこともでき、この方が簡潔な気がした。

import numpy as np

# 2行4列のアレイを作る
arr1 = np.array([[1,2,3,4],[5,6,7,8]])
print("arr1 == \n", arr1)

b = (i ** 2 for j in arr1 for i in j)
print(b.__next__())
print(b.__next__())
print(b.__next__())
print(b.__next__())
print(b.__next__())
print(b.__next__())
print(b.__next__())
print(b.__next__())

結果

arr1 == 
 [[1 2 3 4]
 [5 6 7 8]]
1
4
9
16
25
36
49
64