非圧縮性定常流
ゆく河の流れは絶えずしてしかももとの水にあらず。鴨長明の代表作「方丈記」の有名な冒頭です。哲学的解説は岡部さんに任せるとして、ここでは皆さんには全く馴染みがないであろう「流体力学」のお話をします。私はこの流体力学が大好きです。このコラムでその魅力をお伝えできれば幸いです。※前には伝熱学が大好きと言っていたような？
方丈記のこの一文は見事に「非圧縮性定常流」という流体力学計算の基本的な状態を表した言葉です。川の流れをじ～と眺めてみてください。水はどんどん入れ替わるのに、全く同じ形の流れが常に続きますよね。流体力学とは、流体(気体・液体)が流れる現象を計算で解明するものですが、その際最もシンプルな形がこの「非圧縮性定常流」であり、以下その前提でお話を展開します。
数値解析 (コンピュータシミュレーション)
普通に生活しているとまずお目にかかることのない、流体力学。世の中のどこに役立っているのでしょう？
一番身近なのは天気予報ではないかと思います。昔の天気予報ってあてにならなかったですが、最近の天気はほとんど外れないと思いませんか？これはコンピュータの発達による流体力学の計算精度向上のお陰です。そのほか、車や飛行機のような機械の設計をする際にも、周囲の空気の流れを計算しないといけません。面白いところではイベントの際の人の流れなんかも流体力学でシミュレーションできたりします。
流体のコンピュータシミュレーションとはどういうことでしょうか？例えば天気予報を考えてみましょう。地球というひとつのボールを想像してみてください。このボールから例えば高度1万メートルまでの空間の空気の流れを計算してみます。そうすると、地上の空気の層は、ちょうどイチゴ大福のおもちの部分みたいな感じでしょうか。
この「おもち」の部分を、三次元的に細かく切り刻みます。経線方向、緯線方向、高さ方向という感じです。こうやってできたひとつのバームクーヘンを取り出してみてください。これを計算セルと呼びましょう。
このセルの中に、流れを決定するためのどういった未知数(変数)があるかといえば、
○ 速度 （３成分）
○ 圧力 （一つ）
です。このうち「速度」は大きさだけでなく向きも決める必要があります。こういうものを「ベクトル量」と言って、x,y,zの三方向の大きさを知る必要があります。すなわち、速度を決定するために3つの変数があるということになります。圧力はベクトルではありません(スカラーと言います)ので、結果的に合計４つの変数が決まれば、このセルの流れは決まることになります。※実際は温度も影響し、伝熱学の計算も必要ですがここでは無視します。
気象でいえば、これらは、
○ 速度 （３成分） → 風向きと風速
○ 圧力 （一つ） → 気圧
ということですね。
ここでもう一度、イチゴ大福を思い出してください。細かく切り刻んだセルのひとつひとつの速度(風向きと風速)と圧力(気圧)がわかれば、天気図が描けると思いませんか？
ところがこれらの計算は教科書通りに簡単に答えが出るものではなく、複雑な計算式をコンピュータによって何度も何度も繰り返しを行う必要があります(※)。そして、計算セルは細かく切り刻むほど計算の精度は上がっていきます。したがって、コンピュータは発達するほど天気予報の精度が上がるというわけです。
※例えば、以下の方程式を解いてみましょう。
X2－540X＋60131＝０
この式を(X－157)(X－383)＝０という風に因数分解できれば答えは簡単にわかりますよね。しかし、因数分解の答えってコンピュータの暗号に使われるくらい、実は見つけるのが難しいんですね。じゃあ、どうやって解くかと言えば、Xにひとつずつ数字を入れるわけです。まず1で試すと、これは小さすぎ。次に500で試すと、これは大きすぎ。という感じで何度も何度も繰り返します。そして右辺が0.1より小さくなれば、それを答えにする。といった感じです。こういう考え方を数値解析と言います。数値解析のためには繰り返し計算が必要ですから、コンピュータの性能が求められるわけですね。閑話休題。
流体を計算する
繰り返しますが流体を計算するためには、４つの変数を決定する必要があります。計算セルの中心にある点を想像してください。その点における、
- X方向の速度： Ux
- Y方向の速度： Uy
- Z方向の速度： Uz
- 圧力： P
を決定するためには、４つの方程式が必要です。※変数の数と同じだけ方程式が必要ですよね。
それが、連続の式、ナビエ・ストークス方程式と言われるものです。まずは、書いてみます。
○連続の式
○ナビエ・ストークス方程式
さあ、意味が分かりませんね！大丈夫です。
・６をひっくり返したような記号は変数ではなく数学的記号(偏微分記号)なので気にしないでください。
・x,y,zってのは計算セルの位置を表すもので、これらは変数ではないので気にしないでください。
・Uの右下にある小さなx,y,z(添え字)は方向を表すもので、気にしないでください。
・ρってのは流体の比重を表す数字(例えば水なら１ですね)ですので、これも決まっていて変数ではありません(非圧縮性)。
そうすると、これらの中で変数は、
X方向の速度： Ux
Y方向の速度： Uy
Z方向の速度： Uz
圧力： P
動粘性係数： ν (victoryのvに見えるかもですがギリシャ文字のニューです)
の５つということになります。
４つの方程式で、ここにある５つの変数を見つけてあげれば、あるセルにおける流体の状態を知ることが出来ます。厄介なやつが出来てきましたね！動粘性係数νです。あとで話します。
連続の式とは？
これこそ、冒頭の「ゆく河の流れは絶えずして しかももとの水にあらず」を表した式です。
この数式をどうやって説明しようかと思いましたが、諦めました。これは、セルの中の流体の量が常に一緒ということを表した式です。入ってくる流体の量＝出てゆく流体の量、という式です。
したがってこれは「質量保存の式」とも言われます。
ナビエ・ストークス方程式とは？
いや、実に美しい式だと思いませんか！！私だけ？(汗)
これらについても式の中身の一つ一つの説明はあきらめます。これらは「運動量保存則」とも言われます。で、運動量っていったい何？高校の物理で習ったような？？？
この「運動量」とは、動いている物体の「勢い」を表すもので、質量ｘ速度で表されます。セルの中にある流体の勢いは前後の圧力差とまわりとの摩擦力で決まる、という式だと解釈してください。ちょっとイメージはしづらいかもですね。
簡単にするためセルの形をバームクーヘンではなくキューブだと考えます。このX方向だけに着目してみます。そうすると、この流体の勢い(運動量)を決めるのは、まず左右の圧力の違いです。これが押す力になりますね。
それから、流体には粘性(ねばり)というものがあり、このキューブは周囲の流体にひっぱられます。隣のキューブが自分よりも早ければ推進力となり、遅ければ抵抗力となります。先ほど出てきた、νというのがこの粘り気を表す数値と言うことになります。
※今単純化しましたが、実際は、速度はx,y,zの三方向(ベクトル)ですし、粘性力は一つの速度に対して6面分(前後の面も引っ張りますので)、その差分は全部で9成分ある(これをテンソルという)ことになります。
ナビエ・ストークス方程式を見てみると、
○左辺→ 運動量
○右辺のPが入ってるやつ→ 圧力による推進力
○右辺のνが入ってるやつ→ 粘性による推進力・抵抗力
ということです。例によって式の数学的な説明は割愛します。
しかし、とにかくいい感じのところまで来ました。４つの変数に４つの方程式となれば、あとは数値解析でなんとか答えは出そうですが、まだ変数が５つあります。さあいよいよあとは、動粘性係数νですね。
k-εモデル
動粘性係数νって何でしょうか？これは、流体の持つ粘っこさのことです。お水はサラサラですが、水あめはねばねばですよね。この粘っこさが、流体の運動量＝勢いを決める一つの要因になります。であれば、動粘性係数νは流体の種類によって簡単に決まるんじゃないの？と考えたくなりますが、そうはいきません。なぜか？
それは、ほとんどの流体では「渦」が発生するからなのです。渦によって流体の見かけ上のねばねば感が変わってしまいます。これを、「渦粘性」と言います。この渦粘性は大変複雑で単純なモデルで表すことが出来ません。この渦に対して、またまたたくさんの方程式が必要になるわけです。
こうした渦粘性を計算する代表的な手法が、k-ε方程式(k-εモデル)と言われるものです。またややこしくなってきましたね。
これも細かい説明はしませんが、式中のkは渦のエネルギーで、εは渦が消えてゆく割合です。これらの式の中で出てくる色々な未知数には、実験で得られた数字が入ります。これらもナビエ・ストークス方程式に似た複雑な形をしていますね。
これらk(渦エネルギー)の方程式とε(渦が消えてゆく割合)の方程式の連立方程式を数値解析することにより、ようやくkが決定し、その結果νt(渦粘性係数)が決まります。上では省略した書き方をしていますが、kの方程式も、εの方程式も、いずれも三方向に対してあります。本当にたくさんの方程式です。
これらをさらに、ナビエ・ストークス方程式に戻して数値解析すれば、めでたく以下の値が決まります！
X方向の速度： Ux
Y方向の速度： Uy
Z方向の速度： Uz
圧力： P
しかしながら、これらの手続きは、ある一つのセルの速度と圧力が計算出来たに過ぎません。これをすべてのセルに対してやっていく必要があります。
そして、そして、一通り計算出来たらどうなるか？さっきのセルの周囲の環境が変わっているではないですか！前後の圧力もさっきと違うし、左右から受ける粘性の力も変わってるじゃん！
ということで、その状態からまた一からやり直しです・・・。何度計算しても変わらない状態になった！これこそ「ゆく河の流れ」の状態ですね。これでひとまず流れが決定されるわけです。
しかしながら、天気図は刻一刻と変わりますね。そうですね、天気予報の計算は流れが時間とともに変わる「非定常流」なわけですね。こうなると、時間を進めて同じ処理を行わないといけません。当然、時間間隔を短くすれば制度は上がりますが、その分コンピュータの性能が必要になります。
なので、天気予報の計算手順としては、
- 計算する空間をメッシュに区切り、たくさんのセルをつくる。
- 全セルに適当な圧力、速度をセットする。
※ここですでに分かっている条件を固定値としてセット。これを境界条件という。 - あるセルを選ぶ。
- εの方程式を数値解析。
- kの方程式を数値解析。
- ナビエ・ストークス方程式を数値解析。
- 隣のセルに移る。
- ④～⑦を全セルに繰り返す。
- 全セルに新しい圧力、速度がセットされる。
- ③～⑨を繰り返す→前回の計算結果との差がどんどん小さくなる。
- それが十分小さくなる。→終了
- 時間を進める
- ①～⑫を繰り返す。
- 計算したい時間に到達。→終了
という、気が遠くなるような計算を繰り返すことで、天気予報が可能になるわけですね。こうなると、コンピュータ性能がものを言うため、「京」や「富岳」のようなスーパーコンピュータ性能が必要になるわけですね。
眞鍋淑郎先生
今や愛媛の星である眞鍋淑郎先生。彼は50年以上前に非常に単純化した気候モデルをつくり、当時のコンピュータを使って解析してCO2上昇による地表の気温変化をとても精度よく予言しました。
1950年代、近代コンピュータの生みの親ともいえるジョン・フォン・ノイマンにあこがれた眞鍋先生は、日本(東京大学)を飛び出し、ノイマンが指揮するプリンストン高等研究所で気象予測の研究に従事します。それは、ノイマンが気象学こそが最新コンピュータが最も活躍できる分野だと考えていたからです。
上で説明した通り、今では複雑な気象モデルをたくさんの三次元空間セルに対して、何度も繰り返し数値計算が出来ます。しかし、眞鍋先生の当時は最新コンピュータといえどもこうした計算は到底出来ませんでした。
そこで眞鍋先生がやったモデルは、まず一次元モデルを作ること。地球上のある一点に着目し、Z方向(高さ方向)のみにセルを重ねてゆく手法です。そして、k-εのような複雑な渦モデルは出来ませんので、理論と予測に基づきシンプルな伝熱・流体モデルを使ったところ、とてもよく実際の現象と合ったとのことです。
1960年代に入りコンピュータの発達に伴って、ここで紹介したような三次元の流体計算を行うことが出来るようになります。いろいろなパタンで地球の気象予測を計算している中で、眞鍋先生は「CO2の濃度が二倍になったら地球の平均気温が2℃上がる」という計算結果を導きます。そしてその30年後、実際にCO2濃度が1.5倍になると気温が0.8℃上昇したことが現実となり、この計算モデルの精度がとても高かったことが証明されました。地球温暖化・気候変動が世界的な問題となったためにノーベル賞につながるというのは、ある意味で皮肉なものですね。
結びに
今回も完全に自分の世界に入ってしまいましたね。反省していま・・・せん。時期的にもコラム大賞の対象外！来年の今頃はどうせみんな忘れてる。であれば、書き買いことを書くぞ！大衆には迎合しない！・・・と、強がってます。次のコラム大賞では「サイエンス部門」作ってください！え、たてヨコ・イグノーベル賞ですか？
TATE NAGAFUCHI
