流体の運動方程式

油断してないか? 前回出てきたガウスの定理は今回も使うよ!

[
前の記事へ]  [流体力学の目次へ]  [次の記事へ]


ニュートンの運動方程式を当てはめる

 流体の速度さえ分かればいいと思っていたのに、流体の各点での速度を表す 3 つの関数に加えて、各点での密度を表す関数の合計 4 つの未知関数が出てきてしまった。
\[ \begin{align*} v_x\,&(x,\, y,\, z,\, t) \\ v_x\,&(x,\, y,\, z,\, t) \\ v_x\,&(x,\, y,\, z,\, t) \\ \rho\,&(x,\, y,\, z,\, t) \end{align*} \]
 これらを制限するための方程式は今のところ「連続の方程式」一つだけである。他の条件をもっと追加して現実の振る舞いに近付けないといけない。

 前から言っているように、流体力学はニュートン力学の応用なのであるから、運動方程式を当てはめることにしよう。空間の各成分について、同じ形の方程式が 3 つ作れるはずだ。

 しかしオイラーの方法というのは流体の一部分に注目してその行き先を時々刻々追いかけて見ているわけではないのだった。普通の加速度の概念を当てはめることは得意ではない。\( F = ma \)という形の運動方程式では表せないのである。

 心配は要らない。ニュートン力学はもう少し違う形でも表せる。運動量の時間変化が力であるという考え方が出来るのだった。

 さて、ここでちょっと困ったことがある。運動量は物理では\( \Vec{p} \)で表すことが多いのだが、圧力も\( p \)で表すことが多い。「pressure」の頭文字の\( p \)だろう。流体力学では圧力が良く出てくるのでこちらを優先させたい。仕方ないので今回は大文字の\( \Vec{P} \)を使って運動量を表すことにしよう。これで今言った関係が次のように表せる。
\[ \begin{align*} \Vec{F} \ =\ \dif{\Vec{P}}{t} \end{align*} \]
 この形に書かれていてもニュートンの運動方程式であることに変わりない。


全運動量の時間変化

 流体の持つ運動量なら、オイラーの方法による速度の表現を使っても簡単に表すことが出来そうだ。ある固定した領域を考えて、その内部の運動量の変化を考えてやることにしよう。

 運動量とは質量と速度を掛けたものであった。微小体積\( \diff V \)内の質量は\( \rho \diff V \)だから、それに速度を掛けた\( \rho \Vec{v} \diff V \)というのが微小体積内の運動量である。これを領域全体で足し合わせればその領域内に含まれる全運動量になる。
\[ \begin{align*} \Vec{P}(t) \ =\ \int_V \ \rho(x,y,z,t) \, \Vec{v}(x,y,z,t) \diff V \tag{1} \end{align*} \]
 ベクトルの積分なんてしたことがないという人もいるだろうが、例えば\( P_x \)を知りたければ\( v_x \)を使って計算すればいいというだけである。つまり、この (1) 式は次のような 3 つの式を一つにまとめて書いたものだと解釈して欲しい。
\[ \begin{align*} P_x \ &=\ \int_V \ \rho \, v_x \diff V \\ P_y \ &=\ \int_V \ \rho \, v_y \diff V \\ P_z \ &=\ \int_V \ \rho \, v_z \diff V \end{align*} \]
 こんな具合にして毎回 3 つの式を書き下すのは面倒だし、(1) 式ではどれか一つの成分だけに注目して話が出来ないのでいまいち不便である。そこで次のようにして添え字\( i \)を使って区別する表し方をしてみよう。
\[ \begin{align*} P_i \ =\ \int_V \ \rho \, v_i \diff V \end{align*} \]
 領域内の全運動量\( \Vec{P} \)\( i \)方向成分はこのように表されるということだ。この式を時間で微分してやれば領域内の全運動量のそれぞれの成分の時間変化が表せる。
\[ \begin{align*} \dif{}{t} P_i \ &=\ \dif{}{t} \int_V \ \rho \, v_i \diff V \\[3pt] &=\ \int_V \ \pdif{(\rho \, v_i)}{t} \diff V \tag{2} \end{align*} \]
 また大した検討もしないで気軽に微分と積分の順序の交換を行っているが、あまり気にしなくても大丈夫だろう。


解説の途中ですが広告です


体積力

 次に、領域内の運動量\( \Vec{P} \)がこのような変化を起こす原因について考えよう。まず一つは、領域内の流体の全てに直接働く力の存在が考えられる。これを「体積力」と呼んだりするが、具体的には流体に掛かる重力などのことである。重力以外にそんな力があるのかというとちょっと今のところ想像が付かないが、後で理論をいじりやすいように今は具体的には書かずに、とりあえず単位体積あたりに\( \Vec{f}(x,y,z,t) \)という力が働くという表現にしておこう。領域全体に働く体積力を\( \Vec{F}\sub{V} \)と表すことにすると、次のように表現できる。
\[ \begin{align*} \Vec{F}\sub{V}(t) \ =\ \int_V \ \Vec{f}(x,y,z,t) \, \diff V \tag{3} \end{align*} \]
 もし重力を考えたいのなら、具体的に使う段階でこの\( \Vec{f} \)のところに
\[ \begin{align*} \Vec{f} \ =\ (\ 0,\ 0,\ \rho g \ ) \end{align*} \]
というものを代入して使えばいいのである。(\( g \)は重力加速度で、この場合には下方向を正としていることになる。)

 この (3) 式も成分ごとに分けて話をしやすいように添字を使って表すことにしよう。
\[ \begin{align*} (F\sub{V})_i \ =\ \int_V \ f_i \, \diff V \tag{4} \end{align*} \]
 以前に、加速中の物体が流体中を駆け抜けるときの様子を表せないのではないかということを気にしていたが、その問題は解決しそうである。「一方向に加速する流体」が「静止した物体」の周りを駆け抜けることにすれば等価であり、そのような条件をこのようにして理論に取り込むことが出来るからだ。


面積力

 力は他にもある。領域の表面を通して働く力で「表面力」と呼ばれる。具体的には流体どうしの間に働いている圧力や、粘性によって生じる接線応力のことである。最初の方で話したように、圧力というのは流体内に何らかの面を想定したときに、その面に対して垂直に働く単位面積あたりの力のことであった。静止流体ではそれ以外には力は生じないということだったが、動いている流体では話は別で、粘性による力も生じるのである。粘性についてはまだ詳しく話していないが、これも圧力と同じように、面を通して働く力として扱われる。面の接線方向に働く、単位面積あたりの力として定義されることになる。

 圧力は法線応力、粘性力は接線応力というので、実はどちらも向きが異なるだけの「応力」であり、「応力テンソル」という形にひとまとめにして表すことができる。しかしここで応力テンソルの説明を始めるとめんどくさいことになるので、詳しくは次回やることにしよう。ここではごく簡単な理解をしておくだけで十分である。応力テンソルを\( T_{ij} \)という形で表記する。この\( i \)\( j \)にはそれぞれ\( x \)\( y \)\( z \)のいずれかが入ることになる。
\[ \begin{align*} T \ =\ \left( \begin{array}{lll} T_{xx} & T_{xy} & T_{xz} \\ T_{yx} & T_{yy} & T_{yz} \\ T_{zx} & T_{zy} & T_{zz} \end{array} \right) \end{align*} \]
 これらの成分の意味は単純である。例えば法線ベクトルが\( x \)軸方向であるような面を考えたときに、その面を通して掛かる応力の\( y \)軸成分を\( T_{xy} \)とする、という具合に定義されている。
 二つの添字の意味を逆にして定義している教科書もあるが、実はこのテンソルは対称になっており、添字を入れ替えても同じ値になるのであまり気にしなくても大丈夫である。 どうやら物理の理論では逆の定義を採用している教科書が多くて、流体力学などの工学に近い分野ではこの定義が多いような気がする。 定義の仕方の違いが気になる人はこのあとの説明の応力テンソルの添字を入れ替えて読んでもらえばいいだろう。 なぜ対称になっているのかは次回説明する予定である。
 すると、\( T_{xx} \)というのは\( x \)軸方向を向いた面に掛かる\( x \)方向の力であるからそれは圧力\( -p \)のことである。マイナスが要るのは、ただそういう習慣なのである。物理や工学では物体を両側から引っ張るような「引っ張り応力」を正の値として、両側から押すときの「圧縮応力(圧力)」を負の値で表すことになっている。
 土木、建築分野では逆になっていて、圧縮応力の方を正の値として表す習慣である。
 同様に\( T_{yy} \)\( T_{zz} \)も圧力を意味している。流体中ではどの方向にも同じ圧力が掛かるのだったから、これらは全く同じ値である。粘性が 0 であるような「完全流体」の場合にはこれ以外の応力は存在しないので、応力テンソルは次のようになるだろう。
\[ \begin{align*} T \ =\ \left( \begin{array}{rrr} -p & 0 & 0 \\[3pt] 0 & -p & 0 \\[3pt] 0 & 0 & -p \end{array} \right) \end{align*} \]
 しかし今はこのような具体的な形を使わない方が後で理論の修正や拡張が容易であるし、今から形式的な式変形を行うのにも都合がいいので、テンソルのままで話を続けよう。

 このような応力テンソルは、考えている面が\( x \)軸や\( y \)軸や\( z \)軸の方をまっすぐ向いているような状況しか表せていないのではないかと思うかもしれない。しかしこの応力テンソルさえあれば面がどの方向を向いた場合にでも対応できるのである。法線ベクトルが\( \Vec{n} = (n_x,\, n_y,\, n_z) \)であるような面に働く応力の\( x \)軸成分\( T_{nx} \)が幾つになるかというと、次のように意外と単純な形が成り立つ。
\[ \begin{align*} T_{nx} \ =\ T_{xx} \, n_x \ +\ T_{yx} \, n_y \ +\ T_{zx} \, n_z \end{align*} \]
 なぜこういうことが成り立つかについては今は掘り下げないようにしておこう。とりあえずこの知識を使って計算を続けることにする。領域の表面に掛かる応力を合計してやりたいのだが、ベクトルで書こうとするとややこしいことになるので、とりあえず表面力\( \Vec{F}\sub{S} \)\( x \)成分である\( (F\sub{S})_x \)だけを表してみよう。(S は Surface の頭文字という意味で付けてある。)領域の表面の微小な面\( \diff S \)の法線ベクトルを\( \Vec{n} \)とすると、次のように書けるだろう。
\[ \begin{align*} (F\sub{S})_x \ &=\ \int_S \ T_{nx} \diff S \\[3pt] &=\ \int_S \ (T_{xx} \, n_x \ +\ T_{yx} \, n_y \ +\ T_{zx} \, n_z) \diff S \\ \end{align*} \]
 ところで気付いているだろうか? 法線ベクトル\( \Vec{n} \)というのは領域の外側を向いているという設定であり、領域を外側から押してくる圧力\( p \)というのはそれとは逆向きに働く力である。応力テンソルの対角成分として圧力を負の値で表すのは単なる習慣だと書いたが、今回のような考え方をするときに非常に都合がいいのでわざとそうしているのだろう。

 さて、この式をよく見ると、積分内は\( \Vec{n} \)と応力テンソルの成分で作ったベクトル的な何らかの量との内積である。あまりこれを表すための万国共通の表記はないのだが、ここでは仮に次のようなものを定義してみよう。
\[ \begin{align*} \Vec{T}_x \ \equiv\ ( T_{xx} ,\, T_{yx} ,\, T_{zx} ) \end{align*} \]
 これを使えば次のように書けるし、前回も出てきた「ガウスの定理」が当てはめられるので、面積分から体積分へと書き換えてやろう。
\[ \begin{align*} (F\sub{S})_x \ &=\ \int_S \ \Vec{T}_x \cdot \Vec{n} \diff S \\[3pt] &=\ \int_V \ \Div \, \Vec{T}_x \diff V \\[3pt] &=\ \int_V \ \left( \pdif{T_{xx}}{x} + \pdif{T_{yx}}{y} + \pdif{T_{zx}}{z} \right) \diff V \end{align*} \]
 この最後の式は\( \Div \)の定義に従って書き下しただけである。このように積分内の項の数が増えてくると毎回書き下すのが面倒なので次のようにまとめてしまうことが多い。
\[ \begin{align*} (F\sub{S})_x \ =\ \int_V \ \sum_{j=1}^3 \pdif{T_{jx}}{x_j} \diff V \end{align*} \]
 添字の\( j = 1 \sim 3 \)\( x,y,z \)に対応させて使うのである。分母にある\( x_j \)というのは、\( j \)の値によって\( x\sub{1} \)\( x\sub{2} \)\( x\sub{3} \)となるが、それぞれ\( x \)\( y \)\( z \)の意味である。

 さて、式をできるだけ簡略化して表そうという努力はまだ続く。和の記号\( \sum \)も省略して書くという「アインシュタインの省略」または「アインシュタインの縮約記法」などと呼ばれるルールを適用することも珍しくない。
\[ \begin{align*} (F\sub{S})_x \ &=\ \int_V \ \pdif{T_{jx}}{x_j} \diff V \end{align*} \]
 添字に同じ文字が使われているときには、その文字についての和の記号が省略されていると考えて欲しい。これでかなりすっきりして見えるだろう。これは一般相対性理論では当たり前に使われる省略記法であるし、それ以外の分野でも良く使うものなので、是非慣れて欲しい。もちろん流体力学でも良く使われている。

 ここまで、表面力の\( x \)成分だけを求めてみたわけだが、\( y \)成分や\( z \)成分についても同じ考えで進めることが出来て、次のように表せる。
\[ \begin{align*} (F\sub{S})_y \ &=\ \int_V \ \pdif{T_{jy}}{x_j} \diff V \\[3pt] (F\sub{S})_z \ &=\ \int_V \ \pdif{T_{jz}}{x_j} \diff V \end{align*} \]
 3 つの成分をひとまとめにして表したい時は添字\( i \)を導入して次のように書けばいいだろう。
\[ \begin{align*} (F\sub{S})_i \ &=\ \int_V \ \pdif{T_{ji}}{x_j} \diff V \tag{5} \end{align*} \]
 少々苦労したが面積力も比較的単純な形で表すことに成功した。繰り返すが、和の記号が省略されていることは忘れてはいけない。
 こうして流体に働く力を考えていくと、前にパスカルの原理について考えたときに少し悩んだことを思い出してしまう。 外部から加えられた圧力も重力によって生じる水圧も、どちらも圧力であるという点では同じなのに、 どうして別々に分けて考えないといけないのだろうという点に引っかかったのだった。 外部から加えられた圧力は流体を介して伝わるのに対して、重力は流体に直接働くという本質的な違いがあると気付いて納得したのだった。 そういう話がここでも少し形を変えて「体積力」「面積力」として出てきているわけだ。


解説の途中ですが広告です


運動量を持った流体の出入り

 さて、これで全て終わりだろうか? いや、まだ重要な要素を見落としている。運動量を持った流体が領域の表面から出入りしていて、運動量を持ち逃げしたり運び入れたりしているのである。それも (2) 式のような運動量の時間変化の原因の一つとして数えられる。
 私はこのことでかなり悩んだのだった。 「運動量の時間変化が力に等しい」という考えに固執していたせいである。 領域内の運動量変化は (2) 式ですでに正確に表されているはずだから、あとは力だけを考えればいいはずだという誤りに陥ったのであった。 領域内の運動量の時間変化は (2) 式だけですでに正確であるという考えについては問題ない。 その変化の原因として力以外のものがあるという点を見落としていたのである。

 教科書でラグランジュの方法やラグランジュ微分などを先に見てしまって中途半端に理解していたのも誤解を誘った原因であった。 後で説明する予定だが、ラグランジュ微分を使えば運動量の項をひとまとめに書けるようになる。 その形式を盲信してしまったせいで、今から計算する要素が運動量の補正項のように見えてしまっていたのである。 (2) 式はもう十分に正しくて補正なんて必要ないはずなのに、なぜオイラーの方法でもこんなことを考えなくてはならないのだろうという勘違いをしてしまっていた。
 境界面の付近に速度\( \Vec{v} \)の流体があるとする。境界面の法線ベクトルを\( \Vec{n} \)とすると、この流体は境界面に向かって\( \Vec{v} \cdot \Vec{n} \)の速さで進んでいることになるだろう。つまり、微小面積\( \diff S \)の境界を通って、単位時間あたりに、\( \Vec{v} \cdot \Vec{n} \diff S \)という体積の流体が出ていっていることになる。

 法線ベクトル\( \Vec{n} \)は境界の外側を向いているという設定なのだから、これが正の値なら出ていっていることになるのである。

 この流体は単位体積あたり\( \rho \, \Vec{v} \)という運動量を持っているから、結局、単位時間あたりに、微小面積\( \diff S \)を通って\( \rho \, \Vec{v} \, (\Vec{v} \cdot \Vec{n}) \diff S \)という運動量が流出していると言える。領域全体では次のように書けるだろう。
\[ \begin{align*} \Vec{Q} \ =\ \int_{S} \ \rho \, \Vec{v} \, (\Vec{v} \cdot \Vec{n}) \diff S \end{align*} \]
 この\( \Vec{Q} \)は単位時間に領域から失われる運動量ということなので、運動量の時間微分や力と同じ次元の量である。それをどんな記号で表したらいいものか迷ったがとりあえず\( \Vec{Q} \)を使うことにした。この式もガウスの定理で変形できそうだし是非そうしておきたいのだが、慣れていないと混乱を招きそうである。\( \Vec{Q} \)\( i \)番目の成分という意味で次のように書き直せば少しは分かりやすいだろう。
\[ \begin{align*} Q_i \ &=\ \int_S \ \rho \, v_i \, (\Vec{v} \cdot \Vec{n}) \diff S \\[3pt] &=\ \int_S \ (\rho \, v_i \, \Vec{v}) \cdot \Vec{n} \diff S \\[3pt] &=\ \int_V \ \Div (\rho \, v_i \, \Vec{v}) \diff V \\[3pt] &=\ \int_V \ \left( \pdif{(\rho \, v_i \, v_x)}{x} + \pdif{(\rho \, v_i \, v_y)}{y} + \pdif{(\rho \, v_i \, v_z)}{z} \right) \diff V \\[3pt] &=\ \int_V \ \pdif{(\rho \, v_i \, v_j)}{x_j} \diff V \tag{6} \end{align*} \]
 今回もアインシュタインの省略の記法を使って和の記号を省略している。


解説の途中ですが広告です


まとめて一つの式にする

 ここまで別々に計算してきた要素を一つにしてみよう。(2) (4) (5) (6) 式は次のような関係になっている。
\[ \begin{align*} \dif{P_i}{t} \ =\ (F\sub{V})_i \ +\ (F\sub{S})_i \ -\ Q_i \end{align*} \]
 力が加わればその方向に運動量が増加するからプラスでいいが、\( Q_i \)は運動量が流出する量を意味していたから、このようにマイナスを付ける必要がある。

 それぞれの項を具体的に書けば次のようになる。
\[ \begin{align*} \int_V \ \pdif{(\rho \, v_i)}{t} \diff V \ =\ \int_V \ f_i \, \diff V \ +\ \int_V \ \pdif{T_{ji}}{x_j} \diff V \ -\ \int_V \ \pdif{(\rho \, v_i \, v_j)}{x_j} \diff V \end{align*} \]
 どれも積分範囲は同じなので次のようにまとめられるだろう。
\[ \begin{align*} \int_V \left( \pdif{(\rho \, v_i)}{t} \ -\ f_i \ -\ \pdif{T_{ji}}{x_j} \ +\ \pdif{(\rho \, v_i \, v_j)}{x_j} \right) \diff V \ =\ 0 \end{align*} \]
 この関係が、領域をどのように選んでも成り立っているのだから積分の中の式が常に 0 になると考えていいだろう。
\[ \begin{align*} \pdif{(\rho \, v_i)}{t} \ -\ f_i \ -\ \pdif{T_{ji}}{x_j} \ +\ \pdif{(\rho \, v_i \, v_j)}{x_j} \ =\ 0 \end{align*} \]
 力に関係する項を右辺に集めれば負号も消えてかなりすっきりする。
\[ \begin{align*} \pdif{(\rho \, v_i)}{t} \ +\ \pdif{(\rho \, v_i \, v_j)}{x_j} \ =\ f_i \ +\ \pdif{T_{ji}}{x_j} \tag{7} \end{align*} \]
 もし\( T_{ji} \)となっている部分が美しくないと感じるなら、応力テンソルは添字の入れ替えをしても値が変わらない性質があるので、こっそり\( T_{ij} \)に書き換えておいても構わないだろう。

 この (7) 式が流体の運動方程式だ。名前はまだ無い。この\( f_i \)\( T_{ij} \)のところに流体の性質に関する具体的な条件を入れることでこれから色々な名前で呼ばれることになる。とりあえずそれまでは「流体の運動方程式」とでも呼んでおこうか。あるいは「流体の運動量保存則」と呼んでもいいかもしれない。
追記:この記事を書いた後に色んな教科書を読んでいて「コーシーの運動方程式」という名前がついていることを知った。 他の分野で馴染みのある科学者の名前が流体力学のあちこちに出てくるので驚かされる。 みんな流体力学をやってたのか! 他に「運動量方程式」と呼んでいる教科書もあった。


目的は達成したのか

 (7) 式の\( i \)のところを\( x \)\( y \)\( z \)に変えることで合計 3 つの方程式が得られることになる。我々は新たに 3 つの条件式を得たのである。

 これで 4 つの未知関数に対して 4 つの方程式が与えられたことになるだろうか? いや、残念なことが起きてしまっている。

 先ほども少し見たように、応力テンソルには圧力\( p \)が必ず含まれることになる。この圧力\( p \)は場所と時間の関数\( p(x,y,z,t) \)であり、新たな未知関数がまた一つ加わってしまったことになるのだ。1 つの未知関数と引き換えに 3 つの条件式を得たのだからよくやったと言えるのではないだろうか。

 条件によっては圧力の他に粘性力も加わることになるが、粘性力は流体の速度\( \Vec{v} \)と定数を使って表すことができるので、これ以上の未知関数が増える心配はないだろう。

 我々は今、5 つの未知関数と、4 つの条件式を持っている。もう一つの条件式はどこから持ってこればいいだろうか?

 その話をする前に、次回は応力テンソルについて、もう少し詳しく説明しておくことにしよう。