流体最強の方程式
※この記事は航空宇宙工学科アドベントカレンダー用記事です.
初めに
皆様いかがお過ごしでしょうか.kanepleaseです.実は記事の投稿予定よりめっちゃ遅くなりました.ごめんなさい.
ところで流体,お好きでしょうか?もしかしたら嫌いっていう方いるかもですね.その理由としてはNavier–Stokes方程式(NS方程式)
がやたら長くて覚えにくいとかですかね.よく見るのはこんな感じ.
\begin{equation}
\frac{\partial \boldsymbol{v}}{\partial t}+(\boldsymbol{v}\cdot\nabla)\boldsymbol{v}=-\frac{1}{\rho}\nabla p+\nu \nabla^2\boldsymbol{v}+\boldsymbol{g}
\end{equation}
こりゃ長いよね()ということでNS方程式を忘れても導出できるようにしよう!
さらにどうせならEulerの方程式も導出できるようにしよう!しかも覚えるのは面倒だから一つの式から導出したい!
そんなあなたにお勧めなのがこの方程式
\begin{equation}
\frac{\partial f}{\partial t}+\boldsymbol{v}\cdot\frac{\partial f}{\partial \boldsymbol{r}}+\frac{\boldsymbol{F}}{m}\cdot\frac{\partial f}{\partial \boldsymbol{v}}=\int_{\boldsymbol{v}_{1}}\int_{\Omega '}(ff'_1-ff_1)V σ d\Omega ' d^3\boldsymbol{v}_1
\end{equation}
名を「ボルツマン(Boltzmann)方程式」と言います.この方程式をうまく操作すれば大体の流体の式は導出できます.あれ?NS方程式より長くない?と思ったらそれは幻覚です.この方程式,流体力学のなかでは(おそらく)最強の方程式です.何が最強かって,使用条件が一番少ないんですよね.こいつを完全に解析できたらもう流体力学者いらないと思います.この後はこの方程式についてもう少し話します.
ボルツマン方程式ってすっごーい!
希薄気体解析できるよ!
そもそもこのボルツマン方程式は希薄気体に対する解析でよく用いられます.希薄気体では流体は連続体としてみなすことができなくなり, いわゆる非圧縮のNS方程式が成り立たなくなります.そこで代わりにボルツマン方程式が出番となるのです. 出番っていうのは,この方程式から希薄気体に対する様々なCFD手法が生まれている,ってことです.(例えばDirect Simulation Monte Carlo法)
いろんな方程式生まれるよ!
先ほども言いましたが,ボルツマン方程式からは様々な方程式を生みだすことができます.その一例としてNS方程式やEuler方程式があるんですよね.
それでは今から華麗に導出しましょう!
...と言いたいところですが,その方法は結構複雑なので今回は方針だけ述べておきます.気が向いたら次とかに具体的な計算とか示しますね.
まず,クヌーセン数()という無次元数を導入します.\lambdaを平均自由行程(分子一個が他の分子に衝突するまでに進む平均距離),を
物体の代表長として以下のようにして与えられます.
\begin{equation}
Kn = \frac{\lambda}{L}
\end{equation}
このクヌーセン数というのは流体の希薄度を示していて,例えばよく言われているのがの時には連続体としてみなせるのでNS方程式が
使えるようになるとかです.
次にこのクヌーセン数周りに先ほどのボルツマン方程式を展開(?!)し,右辺と左辺で係数を比較するとオイラー方程式やらNS方程式やらが出てきます.展開ってのは要するにテイラー展開とほぼ同じです.しかしおそらく意味不明ですよね,なのでこれは後日説明ということで...
このように展開してボルツマン方程式から種々の方程式を導出する方法は昔いろんな人が好き勝手に示して,それらを曽根さんという方が「漸近理論」によってまとめ上げました.これによって流体力学がすっきりしたことは間違いなしです.
導出しよう!
いかがでしょうか?ボルツマン方程式のすばらしさ,わかってもらえましたでしょうか?
今回の記事の本題はここからです.ボルツマン方程式の強さが分かったら,そもそもどのような経路でこの方程式が生まれたのか知りたくありませんか?
ぜひ導出しましょう!思ったより簡単なので!
無衝突ボルツマン方程式
まずは簡単な場合を考えて,分子同士の衝突がないものとして方程式を立ててみましょう. 分子の速度の分布関数を以下のように定義しましょう. \begin{equation} f(\boldsymbol{r},\boldsymbol{v},t)d^3\boldsymbol{r}d^3\boldsymbol{v} = \begin{pmatrix} \mbox{時刻}t\mbox{において質量中心が}\boldsymbol{r}\mbox{と}\boldsymbol{r}+d\boldsymbol{r}\mbox{の間にあって}\\ \boldsymbol{v}\mbox{と}\boldsymbol{v}+d\boldsymbol{v}\mbox{の間の速度を持っている分子の平均の個数}\end{pmatrix} \end{equation} それではこの分布関数がどのような関係式を満たすか考えてみましょう.それぞれの分子一つ一つがを受けるものとします. 時刻で位置と速度がとの近くの領域の中の値を持っている分子を考えます. 微小時間後のにこれらの分子が力のもとで運動する結果,それぞれの近くの 領域の中の位置と速度を持つことになります.つまり,図にするとこんな感じ.
で,このことを式で表すと下のような感じ \begin{equation} f(\boldsymbol{r'},\boldsymbol{v'},t')d^3\boldsymbol{r'}d^3\boldsymbol{v'}=f(\boldsymbol{r},\boldsymbol{v},t)d^3\boldsymbol{r}d^3\boldsymbol{v} \label{f'f} \end{equation} なお時間後に関しては以下の二式が自明に成立する. \begin{eqnarray} \boldsymbol{r'}&=&\boldsymbol{r}+\dot{\boldsymbol{r}}dt=\boldsymbol{r}+\boldsymbol{v}dt\\ \boldsymbol{v'}&=&\boldsymbol{v}+\dot{\boldsymbol{v}}dt=\boldsymbol{v}+\frac{1}{m}\boldsymbol{F}dt \end{eqnarray} またここで,ヤコビアンを計算すると(ここは割愛) が成立することより,式\eqref{f'f}から \begin{equation} f(\boldsymbol{r'},\boldsymbol{v'},t')=f(\boldsymbol{r},\boldsymbol{v},t') \label{fddfdd} \end{equation} よって, \begin{equation} f(\boldsymbol{r}+\dot{\boldsymbol{r}}dt,\boldsymbol{v}+\dot{\boldsymbol{v}}dt,t+dt)-f(\boldsymbol{r},\boldsymbol{v},t)=0 \end{equation} これを偏微分で表現すると, \begin{align} Df &\equiv \frac{\partial f}{\partial t}+\dot{\boldsymbol{r}}\cdot\frac{\partial f}{\partial \boldsymbol{r}} +\dot{\boldsymbol{v}}\cdot\frac{\partial f}{\partial \boldsymbol{v}}\\ &= \frac{\partial f}{\partial t}+\boldsymbol{v}\cdot\frac{\partial f}{\partial \boldsymbol{r}} +\frac{\boldsymbol{F}}{m}\cdot\frac{\partial f}{\partial \boldsymbol{v}}=0 \label{nocol} \end{align} という感じになります.これが無衝突ボルツマン方程式になります.ここで,このような形で物質微分的な項が出てきてるのは面白いですね. ボルツマン方程式がNS方程式につながっていることを考えると,この物質微分っていうのは数学的には当ったり前だったんだと今更理解しました()
はい,ひとまず無衝突ボルツマン方程式はここまで,簡単ですね.
二体衝突の記述
ここからは分子の衝突に関して考察しましょう. まず,「分子が単分子でないときでも分子の内部運動の状態,例えば回転や振動などは衝突によっては影響を受けない」 と仮定します. したがって,考えている二つの粒子はそれぞれ質量,位置ベクトル,速度 を持つ単なる粒子と考えることができます. ま,要するにただの古典力学の二体衝突にしてしまおうってこと.この仮定ではなく本気で散乱の話をしながら取り組むのもごくまれにありますが, そこまで難しくする必要はあまりないですね.ここからの二体衝突の話は調べればすぐわかるので,割愛して結果だけ.
相対速度を,質量中心速を,換算質量をとする. すると,を以下のように表現することができる. \begin{align} \begin{cases} \boldsymbol{v}_1 = \boldsymbol{c}+\cfrac{\mu}{m_1}\boldsymbol{V} & \\ \boldsymbol{v}_2 = \boldsymbol{c}-\cfrac{\mu}{m_2}\boldsymbol{V} \end{cases} \end{align}
散乱断面積
2個の粒子の衝突は「散乱断面積」によって記述されます. 衝突による散乱がお互いどのくらいの面積の範囲内で有効になるかということですね. 先ほどまでと同様に,それぞれの粒子を粒子1,粒子2とする. 「標的」粒子2が静止している基準の座標系で,標的粒子2に向かって単位面積単位時間当たり一様な 粒子束で入射する1の型の粒子を考えることにしましょう.
散乱過程の結果,標的粒子から遠く離れたところで1の型の粒子がとの間にある最終速度
をもって単位時間当たり個出てくる.
そこで,散乱されたビームの方向周りの小さい立体角を定める.
この個数は入射粒子束と立体角に比例する.つまり
\begin{equation}
d\mathfrak{N} = \mathfrak{F}_1 σ d\Omega '
\end{equation}
と書くことができます.ここの比例因子σのことを「微分散乱断面積」と呼びます.微分散乱断面積は一般に入射粒子の相対速度の大きさと
入射方向に相対的な散乱ビームの特定方向に依存します.
さらにここで,散乱過程の結果を記述するために,微分散乱断面積とは別に以下のような記号を導入します.
\begin{equation}
σ '(\boldsymbol{v}_1,\boldsymbol{v}_2 \to \boldsymbol{v}'_1,\boldsymbol{v}'_2)\,d^3\boldsymbol{v}'_1d^3\boldsymbol{v}'_2
\equiv\begin{pmatrix}
相対速度\boldsymbol{V}をもって2の型の分子に入射する\\
1の型の分子の単位粒束当たり,\\
散乱後それぞれの最終の速度が\boldsymbol{v}'_1と\boldsymbol{v}'_1+d\boldsymbol{v}'_1の間,\\
および\boldsymbol{v}'_2と\boldsymbol{v}'_2+d\boldsymbol{v}'_2の間の値を\\
もって単位時間当たりに出てくる分子の数
\end{pmatrix}
\end{equation}
ここでも式(\ref{v1v2})が成立している.また,運動量とエネルギーの保存により, である.つまり,例えば と がこれらの条件を満たすものでなければσ'は零になります.
また,散乱後のはに対する天頂角と方位角
によって完全に指定されます.
したがって,すでに式(\ref{N_sigma})で導入してある微分散乱断面積を以下のように再び定義しなおすことができます.
\begin{equation}
σ (\boldsymbol{V'})d\Omega ' \equiv
\begin{pmatrix}
散乱後角\theta 'と\phi 'の周りの立体角d\Omega 'の \\
領域内の方向をもつ最終の相対速度\boldsymbol{V'}\\
をもって単位時間あたり相対速度\boldsymbol{V}で\\
2の型の分子に流\mbox{入}する\\
1の型の分子の単位流束当たり\\
出てくる分子数
\end{pmatrix}
\end{equation}
そして,これら二つの定義を結びつける必要があります.
全ての散乱する粒子の個数に関して以下のように式を立てることができます.
\begin{equation}
\int_{\Omega '}σ (\boldsymbol{V'}) d\Omega ' = \int_{\boldsymbol{c'}}\int_{\boldsymbol{V'}}
σ ' (\boldsymbol{v}_1,\boldsymbol{v}_2 \to \boldsymbol{v}'_1
,\boldsymbol{v}'_2)\, d^3\boldsymbol{v}'_1 d^3\boldsymbol{v}'_2 \label{omegav1v2}
\end{equation}
ただし,積分はとのすべての値について行う.
続いて先ほどの式と同じようにヤコビアンの計算を行うことにより以下の二つの式を導くことができます. \begin{align} d^3\boldsymbol{v}_1 d^3\boldsymbol{v}_2 &= d^3\boldsymbol{c}\,d^3\boldsymbol{V} \label{v1v2} \\ d^3\boldsymbol{v'}_1 d^3\boldsymbol{v'}_2 &= d^3\boldsymbol{v}_1\,d^3\boldsymbol{v}_2 \end{align}
散乱断面積の対称性
特に剛体球を考えれば自明なんですが,散乱断面積は対称性を多く含んでいます.
運動方程式は時間の符号をと反転しても不変でなければならない. \begin{equation} σ '(\boldsymbol{v}_1,\boldsymbol{v}_2 \to \boldsymbol{v'}_1,\boldsymbol{v'}_2) =σ '(-\boldsymbol{v'}_1,-\boldsymbol{v'}_2 \to -\boldsymbol{v}_1,-\boldsymbol{v}_2) \end{equation}
運動方程式はすべての空間座標の符号を逆にする変換とする変換に対して不変でなければならない. \begin{equation} σ '(\boldsymbol{v}_1,\boldsymbol{v}_2 \to \boldsymbol{v'}_1,\boldsymbol{v'}_2) =σ '(-\boldsymbol{v}_1,-\boldsymbol{v}_2 \to -\boldsymbol{v'}_1,-\boldsymbol{v'}_2) \end{equation}
運動方程式は反転衝突に対して不変でなければならない.反転衝突っていうのは,もとの衝突では速度と で衝突し,速度とで出てくるのが, 反転衝突においては速度とで衝突し,速度とで出てくる 衝突のことです. 反転衝突はもとの衝突から(1)時間反転,(2)空間反転の順に操作を行うことで得られる. \begin{equation} σ '(\boldsymbol{v}_1,\boldsymbol{v}_2 \to \boldsymbol{v'}_1,\boldsymbol{v'}_2) =σ '(\boldsymbol{v'}_1,\boldsymbol{v'}_2 \to \boldsymbol{v}_1,\boldsymbol{v}_2) \end{equation}
ボルツマン方程式の導出
さて,長い前置きがありましたが,ここからボルツマン方程式の本格的な導出に入ります.
といっても考え方はとてもシンプルなものです.
もし衝突が無ければ\eqref{nocol}のような関係が成立するのですが,衝突があると,そうはいきません.
なぜなら,衝突の結果もともとこの領域
の中には入らない位置と速度を持っていた分子が散乱によってこの領域に放り込まれ,逆にもともとこの領域の中に入っていた
分子が散乱によってはじき出されることがあるからです.(位置は変わらなくても,速度は大きく変化することがあり得るから)
このような衝突の結果領域の中の分子数が単位時間当たりに増加する正味も個数を
で表すことにする.そうすると,\eqref{fddfdd}は以下のように書き換えることができます.
\begin{equation}
\label{boltzeasy}
Df=D_Cf
\end{equation}
続いて,このに対する具体的な表式を求める.衝突によって生じるの変化率を計算するために 次のような仮定を設ける.
- 気体は十分希薄で,二体衝突を考えるだけでよい.
- 衝突断面積の大きさに与える可能性のある外力\boldsymbol{F}の影響はすべて無視できる
- 分布関数は分子が衝突している程度の時間内では問題になるほど変化せず, また,分子間力の及ぶ領域程度の空間内では問題になるほど変化しない.
- 2個の分子間の衝突を考えるとき,衝突する前の初速度の間にあるかもしれない相関を無視することができる.
以上の基本的な近似を「分子の無秩序性」と呼ぶことにします. (この仮定って結構適当な気がするんですよね,だってボルツマン方程式からNS方程式出てくるんだよ?おかしくない?)
それでは正味の変化を導出します. 次のように考えましょう.
- もともとの中に入っていた分子がの速度領域から衝突で出ていく分を,
- もともとの中に入っていて,の速度領域に入っていなかったが,衝突でこの速度領域内に入ってくるために起こる分子数の増加を
としましょう.つまり, \begin{equation} \label{dc} D_Cf = -D_C^{(-)}f + D_C^{(+)}f \end{equation} と表されます.
衝突項の具体的な計算
を計算するために,体積要素の中に入っていて,近くの速度を持つ分子(これらを分子と呼ぶ)
が同じ体積要素の中に入っていて,ある速度を持つ他の分子(それらを分子と呼ぶ)との[tex: dt}時間内に発生した衝突によって散乱されて,
この速度領域から出ていく場合を考えます.
分子がその速度をからに変え,一方で分子はその速度を
からに変えるような衝突が生じる確率は散乱確率によって記述される.
ここで衝突によって減少する分子の総数
を求めるには,まず分子に入射する分子の流束
に散乱確率を掛け,
さらにこのような散乱をする分子の数を掛け合わせなければならない.
そして,それを分子のすべての可能な初速度について足し合わせ,さらに散乱された分子と分子の
全ての可能な終速度とについて足し合わせる.このようにしては
以下のように記述することができる
\begin{multline} D_C^{(-)}f(\boldsymbol{r},\boldsymbol{v},t)d^3\boldsymbol{r}d^3\boldsymbol{v}dt = \\ \int_ {\boldsymbol{v'}_1} \int_ {\boldsymbol{v'}} \int_ {\boldsymbol{v}_1} \left[|\boldsymbol{v}-\boldsymbol{v}_1|f(\boldsymbol{r},\boldsymbol{v},t)d^3\boldsymbol{v}\right] \left[f(\boldsymbol{r},\boldsymbol{v}_1,t)d^3\boldsymbol{r}\,d^3{v}_1\right] \left[σ '(\boldsymbol{v},\boldsymbol{v}_1 \to \boldsymbol{v'},\boldsymbol{v'}_1)d^3\boldsymbol{v'}\,d^3\boldsymbol{v'}_1\right] \label{dcminus} \end{multline}
を計算するために,同じように体積要素を考える.衝突を繰り返したのち最後に何個の分子が と\boldsymbol{v}+d\boldsymbol{v}]の間の領域内の速度を持つことになるか考える. すなわち,内にある分子のなかで任意の初速度とを持っていて,衝突後には一方の分子が ここで注目しているとの間の領域内の速度をとり,他方の分子はと の間の速度をとるようなすべての分子を考えるのである. この散乱過程は散乱確率によって記述される.の近くの 初速度をもつ分子の相対的な流束はであり,これらの分子は の近くの速度をもつ個の分子によって散乱される. したがって,時間内にの中にあってとの間の速度をもつ分子数の 全増分に対する表式は
\begin{multline} D_C^{(+)}f(\boldsymbol{r},\boldsymbol{v},t)d^3\boldsymbol{r}d^3\boldsymbol{v}dt =\\ \int_ {\boldsymbol{v}_1} \int_ {\boldsymbol{v'}_1} \int_ {\boldsymbol{v'}} \left[|\boldsymbol{v'}-\boldsymbol{v'}_1|f(\boldsymbol{r},\boldsymbol{v'},t)d^3\boldsymbol{v'}\right] \left[f(\boldsymbol{r},\boldsymbol{v'}_1,t)d^3\boldsymbol{r}\,d^3{v'}_1\right] \left[σ '(\boldsymbol{v'},\boldsymbol{v'}_1 \to \boldsymbol{v},\boldsymbol{v}_1)d^3\boldsymbol{v}\,d^3\boldsymbol{v}_1\right] \end{multline} と書くことができる.ただし,積分は分子のすべての可能な初速度と及びの近くの問題 にしている領域内の速度で終わらない他の分子のすべての可能な終速度についてとる.
ここで,反転衝突の確率は元の衝突と確率が等しいことに注意しておく.つまり, \begin{equation} σ '(\boldsymbol{v'},\boldsymbol{v'}_1 \to \boldsymbol{v},\boldsymbol{v}_1) =σ '(\boldsymbol{v},\boldsymbol{v}_1 \to \boldsymbol{v'},\boldsymbol{v'}_1) \end{equation} ですね.
さらに,相対速度 \begin{align} \boldsymbol{V} &\equiv \boldsymbol{v} - \boldsymbol{v}_1 & \boldsymbol{V'} &\equiv \boldsymbol{v'} - \boldsymbol{v'}_1 \end{align} を導入することで,弾性衝突に対するエネルギー保存は \begin{equation} |\boldsymbol{V'}| = |\boldsymbol{V}| \equiv V \end{equation} と表現できる.さらに,ここで略号として \begin{align} f &\equiv f(\boldsymbol{r},\boldsymbol{v},t) & f_1 &\equiv f(\boldsymbol{r},\boldsymbol{v}_1,t) \\ f' &\equiv f(\boldsymbol{r},\boldsymbol{v'},t) & f'_1 &\equiv f(\boldsymbol{r},\boldsymbol{v'}_1,t) \end{align} を導入すると式(\ref{dc})は \begin{equation} \label{DCdefine} D_Cf = \int_ {\boldsymbol{v}_1} \int_ {\boldsymbol{v'}_1} \int_ {\boldsymbol{v'}} (f'f'_1-ff_1)Vσ '(\boldsymbol{v},\boldsymbol{v}_1 \to \boldsymbol{v'},\boldsymbol{v'}_1)d^3\boldsymbol{v}_1d^3\boldsymbol{v'}\boldsymbol{v'}_1 \end{equation} となる.さらに,式\eqref{omegav1v2},\eqref{v1v2}の変換を用いて,この結果をとこの速度の周りの立体角 によって表すことができる.よって,ボルツマン方程式\eqref{boltzeasy}は以下のような具体的な形 \begin{equation} \frac{\partial f}{\partial t}+\boldsymbol{v}\cdot\frac{\partial f}{\partial \boldsymbol{r}} +\frac{\boldsymbol{F}}{m}\cdot\frac{\partial f}{\partial \boldsymbol{v}} = \int_{\boldsymbol{v}_1}\int_{\Omega '}(ff'_1-ff_1)Vσ d\Omega ' d^3\boldsymbol{v}_1 \end{equation} に書くことができる.ただし,である.
お疲れ様でした
これがボルツマン方程式の導出の概要です. ごめんなさい,結構いくつか割愛することになってしまいましたが,なんとなくでも形はつかめたでしょうか? おそらく「ボルツマン方程式」と検索すると,とても難しい数式で語るものが多いです.ですが工学的にはこれくらいで 十分な理解といえると思います.
この後,このボルツマン方程式からは流体力学的極限を用いてNS方程式を導出をしたり,DSMC法を導くことができます. さらに,流体力学的極限からは,「幽霊効果」と呼ばれる不思議な現象が起こったりしています. 興味のある方はぜひ調べてみてください.
それでは!
追記:今回初めてブログ書いてみたのですが,これ結構大変なんですね.知りませんでした.今回こんなに遅れてしまったのは 数式の打ち込みにとても苦労して,時間がTeXの10倍ぐらい時間がかかったからなんですよね.と同じ感覚で入力してるとMarkdownとぶつかったりして,σとかτとか表示するにはMarkdownを使わないといけないとは...(きちんと[tex:]と書いていればこの問題は発生しない, ただ,毎回これ書くの,だるすぎ) しかも,pmatrix内では「入」という漢字が入力できない?とかいう謎のバグにあったり(たぶん原因は別にあるのだろうけど)matrix内で表示されるときによくわからん空白が入ったりして,(これは半角が混ざるとだめという原因だった)本当に意味不明...とりあえずはてなブログで数式書きまくるのはやめた方がいいのかも. ウサギさんの言う通り自分でMathJaxを埋め込んだ方がよさそうです.