ナビエ・ストークス方程式

Size: px

Start display at page:

Download "ナビエ・ストークス方程式"

えみこのえ
5 years ago
Views:

Vol. 4 ナビエストークス方程式序 Dt = F 1 ρ gradp 1 3 νgradθνδv 会員堀城之これが, ナビエストークス方程式である主に流体力学において用いられる土木工学, 建築工学, 機械工学, 航空工学, 船舶工学, 物理学を履修した方は学んだかもしれない流体に係る発明も少なくないであろうしかし, 大学院の流体研究室或いは,

ナビエストークス方程式の生みの親土木技術者であったアンリナビエが粘性を考慮した方程式を記載した論文を 18 年にフランス学士院に提出したしかし, 技術者の論文は無視されるのはいつの時代でも世の常で, 提出から 0 年以上経過した 1845 年に物理数学者のストークスがナビエとは別個に粘性流体方程式を導いたそこで, 二人の名前をつけてナビエストークス方程式とした

1 Vol. 4 ナビエストークス方程式序 Dt = F 1 ρ gradp 1 3 νgradθνδv 会員堀城之これが, ナビエストークス方程式である主に流体力学において用いられる土木工学, 建築工学, 機械工学, 航空工学, 船舶工学, 物理学を履修した方は学んだかもしれない流体に係る発明も少なくないであろうしかし, 大学院の流体研究室或いは, 当該数学研究室に入らない限り圧縮項 ( 右辺第 3 項 ) までは導かないと思われるそこで, 本稿ではナビエストークス方程式について圧縮項まで含めて説明するなお, 特別な境界条件の場合を除いて未だに誰も解けていない方程式である 1. ナビエストークス方程式の生みの親土木技術者であったアンリナビエが粘性を考慮した方程式を記載した論文を 18 年にフランス学士院に提出したしかし, 技術者の論文は無視されるのはいつの時代でも世の常で, 提出から 0 年以上経過した 1845 年に物理数学者のストークスがナビエとは別個に粘性流体方程式を導いたそこで, 二人の名前をつけてナビエストークス方程式としたナビエとストークスとの間にが付いているのはその為である第一発見者の名前をとってナビエ方程式の方が正しいような気がする因みにナビエは, その後, エリート養成学校であるエコールポリテクニークの教授になっているアンリナビエ from Wikipedia. ナビエストークス方程式が導かれた理由ダランベールの背理右辺第 3 項 ( 圧縮項 ) 及び第 4 項 ( 粘性項 ) を除いた式は, オイラーの運動方程式であるこの方程式は, ナビエストークス方程式の約 100 年前に, あの大数学者オイラーによって, ニュートンの第法則から導かれたここで, 問題です感覚的に分かってもらえばよいので, 細かい条件等は抜きにして大雑把に書きます Q: 一様に定常的な流れるプールにあなたが飛び込んでプールの底に立つとするあなたは水平方向にどのような力を受けますか? あなたはこう答えるでしょう A: 下流に流される力を受ける流れが強いと下流に流される答えはですオイラーの運動方程式に支配される流れでは, あなたは下流方向への力を受けない無風状態で打たれたゴルフボールは, 重力以外に力を受けず空気中を飛んでいくのである非常に奇妙な話であるあなたも流れるプールや川で流されたかもしれない流体中に存在するものは上流から下流に流されるように力を受けることを感覚的に理解しているゴルフをやる人なら流体抵抗が無かったらどんなに楽なことかと思うかもしれないしかし, オイラーの運動方程式で扱う流体では下流方向に流される力を受けないのであるオイラーの運動方程式に支配される流体を完全流体というそして, このパラドックスをダランベールの背理という当該パラドックスを解くために考え出されたのが, ナビエストークス方程式なのであるナビエは, 土木橋梁技術者であり, 橋脚に受ける力を研究した結果, ナビエストークス方程式を導き出したパテント

2 3. ナビエストークス方程式の導出ナビエストークス方程式自体を導くことはそれほど難しくはない高校の数学, 物理の知識があれば導くことが出来る (1) オイラーの運動方程式まず, ナビエストークス方程式の圧縮項 ( 筆者が勝手に名付けた学部学科で名称が異なるようである ) 及び粘性項 ( 右辺第 4 項 ) が無い方程式を導く当該方程式をオイラーの運動方程式という x Dt = F 1 ρ gradp grad =~ y z オイラーの運動方程式は, 右辺は単位質量当たりの力, 左辺は質量力と圧力であるゆえに, オイラーの運動方程式は, ニュートンの第法則 F=Mcから導くことが出来る x x,y y,z z 図 1 図 1に示す, 流体中における微少流体の速度成分は, 場所と時間の関数, u = u(x, y, z, t) v = v(x, y, z, t) w = w(x, y, z, t) 1 代表される時間経過後の移動距離は, x = u(x, y, z, t) y = v(x, y, z, t) z = w(x, y, z, t) 座標は, x x y y z z 3 3を1に代入すると, u = u(x x, y y, z z, t ) v = v(x x, y y, z z, t ) w = w(x x, y y, z z, t ) 4 微少時間経過後の速度の変化量をu として,4を以下のように書き換える u(x, y, z, t)u v(x, y, z, t)v w(x, y, z, t)w 5 5を多変数テイラー展開 ( ) すると, u(x, y, z, t)u = u(x, y, z, t)( x y y z )u(x, y, z, t) ( x y y z ) u(x, y, z, t) 3 乗, n 乗 v(x, y, z, t)v = v(x, y, z, t)( x y y z )u(x, y, z, t) ( x y y z ) u(x, y, z, t) 3 乗, n 乗 w(x, y, z, t)w = w(x, y, z, t)( x y y z )u(x, y, z, t) ( x y y z ) u(x, y, z, t) 3 乗 n 乗影響が小さい高次項を無視し,を代入すると, u(x, y, z, t)u = u(x, y, z, t)( u y v w ) v(x, y, z, t)u = v(x, y, z, t)( v u v y v w v ) w(x, y, z, t)w = w(x, y, z, t)( w u w y v w w w ) 6 u lim 0 として6を整理すると, du dt = u y v w dv dt =v u v y v v w dv dt =w u w y v w w 7 左辺を見れば分かるように, 所謂, 加速度であるを時間的加速度, u y v w を場所的加速度, という前者は理解できると思うが, 後者がなぜ加速度と思うかもしれないつまり場所で微分してなぜ加速度なのかと思うかもしれない場所的加速度は, ある場所から他の場所に移動したときの速度の変化率 ( 加速度 ) だと考えて頂ければよい以上で加速度が求められたつまり, オイラー方程式の左辺である次に, 右辺, すなわち, 力について考える完全流体に加わる力は, 質量力と圧力のみである質量力は, 流体工学の本に出てくる用語なのだが, 筆者が訳すと慣性力流体中に微少な立方体を考える 89 パテント 013

1 粘性係数 μ 剪断係数 τとの関係は, τ=μ v で表されるこれをニュートンの摩擦法則という ( ニュートンの粘性法則ともいう ) 応力テンソル図図に示す立方体の質量は,ρxyz である =v Fx 各面に加わる単位重量当たりの力を F Fy Fz =v Fx とすると, 質量力 = Fρxyz Fy Fz ρxyz 8 次に, 微少流体の中心を (x,y,z) とすると,

3 1 粘性係数 μ 剪断係数 τとの関係は, τ=μ v で表されるこれをニュートンの摩擦法則という ( ニュートンの粘性法則ともいう ) 応力テンソル図図に示す立方体の質量は,ρxyz である =v Fx 各面に加わる単位重量当たりの力を F Fy Fz =v Fx とすると, 質量力 = Fρxyz Fy Fz ρxyz 8 次に, 微少流体の中心を (x,y,z) とすると, 微少流体に作用する圧力は, 例えば x 軸に垂直な面に作用する圧力差であるから, { P (x dx dx, y, z) P (x, y, z)}yz = P xyz 9 他の面に作用する圧力も同様に求めると, P y xyz, P xyz 10,11 7 9をニュートンの第法則に代入し, 質量で除せば u y v w = F x 1 ρ v p u v y v v w v = F y 1 p ρ y w u w y v w w w p = F z 1 ρ D Dt (= u y v w ) を用微分演算子いて書き換えると, ρ Dt = F y 1 p ρ y Dw Dt = F z 1 p ρ grad を用いて記載し直せば, Dt = F 1 ρ gradp これで, オイラーの運動方程式が導けた () 粘性流体における基礎方程式 1 μ 圧縮項 3 ρ gradθ, 粘性項 νδv のオイラー運動方程式への導入ダランベールの背理が生じるのは粘性を考慮していないからであるまた, 圧縮性も考慮されていない因みに, 水は非圧縮性流体である図 3 σは着目面に垂直に作用する応力 ( 材料力学で言えば面外応力 ),τ は着目面に作用する剪断応力 ( 同面内応力 ) である最初の添え字は着目面を, 番目のそれは作用方向を示す以上を行列で表したものが応力テンソル ( 略 ) であるこの図をしっかり覚えて頂ければ終わったようなものです着面に作用する剪断応力は, 係るつの剪断応力の和で表されるので, ニュートンの粘性法則から, τ xy =τ yx =μ v y τ yz =τ zy =μ w y v τ xy =τ yz =μ w 1 また, 着面に作用する面外応力は, 圧力と, 係るつの剪断応力の和と, 圧力による歪みの項との和で表されるので, σ xx = p μ λθ σ yy = p μ v y λθ σ zz = p μ w λθ 13 μを第 1 粘性係数,λを第粘性係数という θ= div v = v y w 非圧縮性流体では歪み 0 故,θ= 0 次いで, 式 1,13を使って, 図 3に示す立方体に働く力を考えるパテント

4 x 軸方向の力は対面の差で表されるから x 面では σ xx x yz y 面では τ yx y xz y Z 面では τ zx z xy であり, これらの総和を単位質量で表すと, 1 σ xx ρ τ yx y τ zx 14 14を 3.(1) で求めたオイラーの運動方程式に代入する Dt = F x 1 ρ σ xx τ yx y τ zx 圧力項が消えていると思われるかもしれないが, 式 13にあるように,σの中に入っている式 15に式 1,13の当該式を代入してやれば, ρ 1 3 νθ ν u u u y ν= μ ρ ν: 動粘性係数 y 面,z 面も同様に求めると, ρ y 1 3 νθ y ν u u u y Dw ρ 1 3 νθ ν u u u y これらの式を grad,δを用いて表すと Dt = F 1 ρ gradp 1 3 νgradθνδv Δ= = y ' 15 = y 非圧縮性流体では, Dt = F 1 ρ gradp νδv 多変数テイラー展開 f(x m )=6 1 n=0n [(x m Y m ) m] n f(y m ) from Wikipedia イギリスの数学者ブルックテイラーが導入した級数式であるある関数の 1 点から少し離れたところの値を級数を使って真値に近づけるものである高校の授業で 1 変数テイラー展開 ( マクローリン展開 (a = 0)) を使って円周率や三角関数を単純な数の和として求めた方もいるかもしれない工学, 物理で用いられる便利な式である重要なので3 次元の場合を導くナビエストークス方程式は 3 次元なので 3 変数場合は, 関数 f(x, y, z) を座標 a, b, c 周りでテイラー展開すると, f(x, y, z)=σ n=0 1 h n! k y j n = 1 h 1! k y j 1 1 h! k y j 1 h 3! k y j 3 1 h n! k y j n R n1 R n1 = 1 h (n 1)! k y j n1 f(a θ h, b θ k, c θ j) 0 <θ< 1 : 多変数ラグランジュ剰余項ここで, x = a ht,y = b kt,z = c jt つまり, 流線上の 1 点 x, y, z から微少時間経過 ( 微小変化 ) した場合 ( 故に 1 次 ) の予測値 ( 座標 ) を表している (f(xx, yy, zz)= f(x, y, z)f(x, y, z)) 導出関数 f(x, y, z) の座標における x, y, z 軸方向の傾きに単位時間当たりの各軸方向における変化量をそれぞれ掛けてやれば変化量 f(x, y, z) が出てくるから f(x, y, z) dw = h k y j dt f(x, y, z)f(x, y, z)= h k y j h k y j の偏微分係数も変化するので, それを考慮すると時間経過後の偏微分係数の変化量は, = f(x, y, z) f(x, y, z) y f(x, y, z) y = f(x, y, z) y f(x, y, z) y f(x, y, z) y = f(x, y, z) f(x, y, z) y f(x, y, z) を二つに分割すると, 前半の後の変化量は, 1 h 1! k y j 1 後半の後の変化量は, 91 パテント 013

5 1 h 1! k y j z 1 y = 1 h 1! k y j z 1 f(x, y, z) f(x, y, z) f(x, y, z) f(x, y, z) y y f(x, y, z) f(x, y, z) f(x, y, z) y y f(x, y, z) f(x, y, z) y 時間経過後の変化量は, 前半と後半との和で表されるから, f(x, y, z)= 1 h 1! k y j 1 1 h 1! k y j 1 f(x, y, z) f(x, y, z) y f(x, y, z) f(x, y, z) f(x, y, z) y y f(x, y, z) f(x, y, z)= 1 1! h k y j 1 y f(x, y, z)= 1 1! h k y j 1 1! h k y j f(x, y, z)f(x, y, z)= h k y j 1! h k y j n =0, 1, についての導出終わり高次項は無視されるのでこれで打ち止めなお, ラグランジュ剰余項を含めて, ロルの定理等でテイラー級数を導入することが出来る Fin 4. 賞金 $100 万ドルナビエストークス方程式は, 上記の通り, 階非線形偏微分方程式である以前, パテントに掲載した振動方程式は階非線形偏微分方程式であり ( 線材は 4 階 ), ナビエストークス方程式と同じである断面変化等の場合を除き, 時間微分のみであるナビエストークス方程式では時間微分と場所微分とが混在しているよって, 振動方程式よりも解きにくい特定の条件下で境界条件を設定して解くことは可能である振動方程式では減衰を無視すれば階線形偏微分方程式となり解くことが出来るナビエストークス方程式では, 例えばマッハのような高速流体で且つ非圧縮性であれば粘性項及び圧縮項を無視できるので1 階線形微分方程式となり解くことは可能であるしかし, 上記の如き特定の場合を除き, 数学的には, 厳密解が求められていない求めることができるのかも分かっていないそこで,CMI( クレイ数学研究所 ) は賞金を懸けている CMI は, 解けないという証明をしても賞金をもらえるそうであるところで,CMI は,006 年にポアンカレ予想を証明できたとして, グリゴリーヤコヴレヴィチペレルマンに賞金を出した ( 本人は受け取り拒否 ) 先日 NHK で数学者はキノコ狩りの夢を見るという題名で放送していたので見た方もいるかもしれないなお, 見逃した人は, オンデマンドで見ることが出来る eman.jp/ ポアンカレ予想とは, 位相幾何学の問題で, 紐の基端を手で持ち先端をロケットに結びつけて飛ばし, 宇宙の隅々まで飛行し地球に帰還した後, 紐の基端と先端とをたぐり寄せると紐の全てを回収できるという予想である ( 星はツルツルで紐は星に引っかからない紐は燃えない, 切れない ) ペレルマンの論文は, 例えば,Cornel University Library でダウンロードすることが出来るペレルマンは, ウィリアムサーストンの幾何化予想を, 微分幾何学に物理工学的アプローチも加えて解き, ポアンカレ予想を解明したサーストンの幾何化予想とは, 宇宙では最大 8つの図形に分類されるというものであるナビエストークス方程式も純粋に数学的アプローチではなく, 他のアプローチにより解けるかもしれないたとえば, 生みの親であるナビエが土木技術者だから土木工学的アプローチも面白そうである土木工学的に, 流体を類型化し, それぞれ圧縮項, 粘性項が 0 とならない境界条件を定めてプロットしていくと, 解法のヒントが見つかるかもしれない以上参考書 : 日野幹雄著流体力学株式会社朝倉書店本間仁著標準水理学丸善株式会社パテント 013 9

パソコンシミュレータの現状

パソコンシミュレータの現状第 2 章微分偏微分, 写像豊橋技術科学大学森謙一郎 2. 連続関数と微分工学において物理現象を支配する方程式は微分方程式で表されていることが多く, 有限要素法も微分方程式を解く数値解析法であり, 定式化においては微分積分が一般的に用いられており. 数学の基礎知識が必要になる. 図 2. に示すように, 微分は連続な関数 f() の傾きを求めることであり, 微小なに対して傾きを表し, を無限に