第3章 粘性流体の運動
先の2つの章にて、弾性理論における応力・歪みテンソル、構成則、そして運動方程式を説明した。最も 重要な式を復習の意味も兼ねて、以下にまとめる(符号などの説明は省略): これらの関係式を少し変更しただけで、流体運動における方程式を簡単に導けることを、この章では示す。 ただし、流体力学は一般には非線形性があり、問題の状況によって扱い方が大きく異なる。気象学や海洋 物理学などに必要な視点からの流体力学は別の詳しい講義があるので、これらに関係することは簡単に触 れるだけにして、固体地球惑星科学において関係する流体力学の視点から考える。つまり、マントル対流、 マグマ流体、post glacial rebound 等に関係する「粘性があって」「運動が比較的ゆっくりした」現象に対し ての流体運動を対象とする。物質科学的な考察も含めた内容は次の章で考える。3.1
粘性流体の構成則と粘性率 Viscosity and Newtonian fluid
ばねにつながれた質点では fi= −kxiというフックの法則があり、弾性体ではこれが のよ うに、力が「応力テンソル」σijに、変位が「歪みテンソル」eklに対応するフックの法則を示した。では、 気体や液体のような流体の場合はどうなるか。一般の流体力学の本では粘性のない理想流体(ideal fluid) から始まり、大部分が粘性のない理想流体の説明で占められている。理想気体ではせん断応力は発生せず、 式 (1.9) の で応力が表現される。この場合には運動方程式そのものを扱わず に解析できる例が多い。気象学や海洋物理学でもこのように扱う例は多い。 これに対して、固体が長い時間かけて流動するような現象では粘性が重要な役割を持ち、いわゆる として扱わないといけない。粘性は摩擦力と同じくエネルギー散逸過程であり、熱 力学的には不可逆過程(irreversible process)に当たる。質点の力学では単振動の際に摩擦が入って減衰振 動となった現象に対応する。すなわち、静止している場合には力はかからず、運動している、つまり速度に 比例して力がかかる形式となる(例:ランダウ・リフシッツ「力学」第25章参照)。弾性体では応力テン ソルは歪みテンソル eij、すなわち変位 ukの空間微分に比例したが、粘性によって生じる力は変位ではな くて速度に比例するから、歪みテンソルの時間微分である「歪み速度テンソル」あるいは「歪み率テンソ ル」(strain rate tensor)に比例する:
(3.1) ここで、vi≡ ∂ui/∂tという速度成分である。(流体力学では歪みテンソルそのものは関係ないので、eijの ような表式を上の定義として採用している場合も多いので注意)。 理想流体にはせん断応力成分はないと前の章では述べたが、粘性流体では速度に比例した力がかかり、さ らに実質的に重要な成分はせん断成分のみである。厳密な導出は専門書に譲り(例:ランダウ・リフシッツ 「流体力学」第15章)、ここでは図3−1の模式的な場合の結果から、一般的な粘性流体の構成則の一般 的なテンソル形式の表現を示す。図3−1のように、平らな面の上に流体が面に沿って流れているとする。 粘性がない理想流体ではせん断応力がないので、水平に一様な速度 v で流れることができ、せん断応力成 分に伴う力はかからない。粘性があると、面(および流れ)に沿ったせん断応力 F が生じて、面上では速
度がゼロで、面からの距離 D によって流れる速度 v が次第に大きくなっていく。D が小さい場合、実験的 にはせん断応力 F と面の距離からの速度の増加率 v/D がほぼ比例すると見なしてよい: F ∝Dv −→ F ≡ ηDv (3.2)
v
(a) 図3- 1 (b)=const.
理想流体D
粘性流体v
F
F=0
このような比例関係にある粘性流体をニュートン流体(Newtonian fluid)と呼んでいる。上の関係は、あ るせん断応力成分 σ と歪みの時間変化(歪み率とも呼ぶ) ˙e が比例関係にあるので、σ = η ˙e ということに なる。一般の成分をテンソル表現すると、フックの法則である式 (2.2) の歪みテンソルを歪み速度テンソル ˙eklと置き換えた関係となるので、 σij = Dijkl˙ekl といった比例定数 Dijklで表現できる。ここでは等方的な粘性流体を考えるので(層状構造をしている固体 の流動は異方性があるなど、最近は異方的粘性流体も地球惑星科学では対象になりつつある)、上の一般的 なテンソル表現ではなくて、等方媒質でのフックの法則 (2.14) の歪みテンソルを歪み速度テンソル ˙eijに置 き換えた表現を考える: (3.3) しかし、慣例的に上式の第2項であるせん断弾性率 µ に対応するような粘性流体でのせん断応力とせん 断歪み速度との比例係数は η という符号で表わし、これを「粘性率」(viscosity)と呼ぶ。また、この式で 引っ張り・圧縮を表わす第1項については流体であるから「圧力」p が重要なパラメターと見なされる。式 (1.9)より等方媒質では応力テンソルの対角成分が −p であるので、これからのずれが理想流体にはなくて粘 性流体のみに生ずる力として、その係数を λ とする。よって、圧力 p と粘性率 η を用いた慣例的な表現は、 (3.4) となる。等方媒質では2つしか独立な変数がないが、圧力を導入したために3つの係数(λ, η, p)になって いる。よって一つは消去できる。応力テンソルの対角和はσxx+ σyy+ σzz= 3λ ˙ell+ 2η ( ˙exx+ ˙eyy+ ˙ezz) − 3p = (3λ + 2η) ( ˙exx+ ˙eyy+ ˙ezz) − 3p
となり、左辺は式 (1.9) から −3p なので、
(3.5) という関係式があることがわかり、結果的に(等方的でニュートン流体である)粘性流体の構成則は以下の ように p と η という二つの係数で表現できる:
(3.6) ここで、 ˙θ ≡ ˙exx+ ˙eyy+ ˙ezzであり、式 (1.25) から「体積変化速度(体積変化率)」に対応する。
実際の流体では、引っ張り・圧縮の応力に対してはほとんど変形しないのに対して、かなり強い粘性があっ てもせん断応力には弱い力でも変形する(流体はせん断弾性率 µ = 0 に対応することに注意)。よって、上式 (3.6)に示す構成則の体積変化率 ˙θ に関する第2項は他の項に比べて著しく小さいことがわかる。これをゼロ
とみなすと体積変化しない (または密度が変化しない) 流体に対応するので と呼ばれ、地球惑星科学ではほとんどこの仮定で十分である。後で示す「連続の式」でも非圧縮流体 はこの値がゼロとなることが示される: (3.7) 以後は、応用上では十分な以下のような形式の構成則が成り立つ粘性流体を考える: (3.8) 粘性率についての物質科学的な視点からの簡単な考察は、第4章で触れる。 [問題3-1] 粘性率 η の単位を求めよ。また、理科年表やインターネットなどを調べて、以下の物質の粘性 率を示せ。水、空気、ハチミツ、摂氏1000度のカンラン石、有珠山のマグマ、ハワイ島の火山の溶岩 (赤く流れる)、固体地球全体の平均値
3.2
流体での微分:ラグランジュ微分対オイラー微分 Lagrangian versus Eulerian
derivatives
前節の粘性流体の構成則である式(3.6)、または (3.8) を弾性運動方程式である式(2.23)に代入すれば、 基本的には粘性流体の運動方程式が得られるはずである。しかし、流体力学の講義や本のまず最初に学習す るように、弾性体とは異なり(特にせん断応力に対して)大変形が可能な流体では、式 (2.23) の左辺の慣 性項の時間微分に特別な注意が必要となる。すなわち、物質に固定した座標系でのラグランジュ微分 D/Dt と、空間に固定した座標系でのオイラー微分 ∂/∂t を厳密に区別する必要がある。詳しい解説は流体力学等 の講義に譲り、ここでは簡単に復習する。 図3−2のように、時刻 t に点 P の位置 X にあった流体の小片が、変位 u によって点 ¯P の位置 x に移 動したとする。つまり、 x (X, t) = X + u (3.9) という関係にある。この時に、流体の変位ベクトル u を X と x のどちらの位置によって考えるかによって、 大変形の際には差が生じることが本質的である:1. u(X, t): ラグランジュ流表現(Lagragian description)
「元の位置、またはやって来た位置としての関数」であり、物質粒子にくっついた座標系となる。速 度を例に取ると、流体に粒子を落としてその動きを写真やビデオカメラで測定した値に対応する。 2. u(x, t): オイラー流表現(Eulerian description)
「今いる位置としての関数」であり、空間に固定した座標系となる。速度を計測する装置を位置 x に 固定して測定した場合の値に対応する。 X P O O x P u(X) or u(x) ? P 図3- 2 X さて、ある物理量(ここではスカラーとするが、ベクトルなどについてもそれぞれの成分と見なせば同じ である)q が空間と時間の関数とする。元の位置 X が基準となるから、普通ならば q(X, t) の形になり、こ れを時間微分したものが「ラグランジュ微分」と定義される。本によっていくつかの表現があるが、ここで
は通常の時間微分と明確に区別するために、Dq/Dt と表示する。一方で、q の値そのものは同じでも変形 した後の位置 x の関数を q(x(X, t), t) として、この場合の時間微分は x を固定した偏微分である ∂/∂t であ り、「オイラー微分」と呼ぶ。ラグランジュ微分を x を関数とするオイラー微分として見ると、 Dq Dt ≡ !∂q ∂t " X −→ dq(x(X, t), t) dt = !∂q ∂t " x + 3 # j=1 ∂q ∂xj ∂xj ∂t = ∂q ∂t + vj ∂q ∂xj (3.10) ここで、最後の速度は vj ≡ ∂xj/∂tと定義されるように、ラグランジュ流の X についての速度ではなくて、 オイラー流の x の時間微分である速度であることに注意。 上の式から、ラグランジュ微分 D/Dt とオイラー微分 ∂/∂t (と流体力学では時間の偏微分は暗黙のうち に x と t の関数として x を固定しての微分とみなす)との以下の関係が示される: (3.11) 速度ベクトル v の時間微分は加速度ベクトルであるが、ラグランジュ流でこれを表すと、 (3.12) 次の節では流体力学の運動方程式(通称、ナビア・ストークス方程式)を示すが、既に第2章で弾性運動 方程式を導いたので、大きな違いは上の関係式に注意する点だけである。しかし、式 (3.12) の右辺の第2 項は速度 v について二乗の形になっており、つまり「非線形微分方程式」という、これまでの初等の物理 (や数学)で学習した線形微分方程式とは根本的に異なった性質を持つ(そして、解析的には難解な)方程 式である点のみ、まずここで指摘しておく。
3.3
粘性流体の運動方程式:ナビア・ストークス方程式: Navier-Stokes equation
第2章で弾性体における運動方程式である式 (2.23) を導入した: ρ∂ 2u i ∂t2 = σji,j+ fi 3.1節において粘性流体の構成則である式(3.8)(圧縮効果が大きい特殊な流体では (3,6))を導き、3.2 節 では流体のように大きな変形をもたらす媒質での時間微分の大きな注意点を指摘した。これで、粘性流体 の運動方程式の導入の準備はすべて整った。 上の式 (2.23) を導く過程を振り返ると、媒質内の微小な直方体にかかる力を考えたわけで、左辺の時間 微分は厳密には空間に固定したものでなく、「その部分の物質に固定した座標」についての微分であること がわかる。よって、厳密には である。弾性理論ではフックの法則によって応力テ ンソルが「歪みテンソル」に比例した、つまり静止していても変位があれば力が生じる形式になっていた。 しかし、流体力学(特に粘性流体)では静止していたら力はかからず(正しくは 2.2 節で説明した静水圧平 衡による圧力のみかかる)、運動していることで応力が初めて生じるので、変位 u ではなくて速度 v で表現 するべきである。以上から、粘性流体の運動方程式は式(2.23)を若干変更して以下のようになる: (3.13) ここで、弾性方程式 (2.23) の場合に D/Dt と ∂/∂t を区別しなかった理由を考える。上に述べたように正式 には式 (2.23) の左辺の慣性項も D/Dt を用いるべきである。しかし、弾性理論で扱う問題は「微小変形」 という条件が基本であることは、第2章のフックの法則に関してや、実際の変形の大きさの見積もり例など で示した。つまり、変位や速度の大きさが微小 |u | , |v | ≪ 1 (3.14)なので、例えば、式 (3.12) の左辺の第2項の大きさは O(v2)と速度の二乗で、他の項は O(v) と一乗なの で、より高次の微少量となり、省略しても構わない: Dv Dt = ∂v ∂t + v · ∇v ≃ ∂v ∂t (3.15) のように、ラグランジュ微分とオイラー微分を弾性理論では区別する必要がないことがわかる。このよう に、式 (2.23) では時間微分はあいまいに扱って構わないのに対して、大変形する流体運動では式(3.13)の ように2種類の時間微分を厳密に区別する必要がある。 式 (3.14) の条件があるために、弾性運動方程式 (2.23) は u について線形であるのに対して、式 (3.13) に は v · ∇v という v について二乗の項が含まれることになり、これは「非線形方程式」であり、初等の数学 や物理で学習した方程式とは全く異なる。非線形方程式は足し合わせができないなどの重要な特徴があり、 解析解を得ることそのものが難しいので、30年前くらいまでは驚くことに基礎的な研究すらほとんどさ れてこなかった。計算機による数値解法による解の振る舞いが詳しく研究されるにつれ、カオスや複雑系な どの近年の数学分野も含めた発展は著しいが、中でも古くから知られ、また室内実験などで検証も可能な ナビア・ストークス方程式は、カオス概念の誕生を例として、非線形方程式の研究の発展において中心的役 割を演じてきた。 非圧縮で粘性率 η が一定とする粘性流体の運動方程式は、構成則である式 (3.8) を式 (3.13) に代入して、 ρDvi Dt = ρ !∂v i ∂t + vj ∂vi ∂xj " = fi− ∂p ∂xi + η∇ 2v i (3.16) となる。ベクトル表示では以下のようになる: ρDv Dt = ρ ! ∂v ∂t + (v · ∇) v " = f − ∇p + η∇2v (3.17) これらの形式をナビア・ストークス方程式と呼ぶ場合も多い。 [問題3−2] 式(3.16)と (3.17) を導け。 運動方程式と並んで流体力学で重要な方程式に あるいは「質量 保存の式」がある。これは図2−4のような小さな直方体を考え、この直方体の質量の変化を、その密度 ρ の変化と6つの面(外側の法線方向を向く単位ベクトルを n とする)を通して流れ込む(負の値ならば流 れ出す)質量が釣り合っていることから導く。直方体の体積 V と6つの面 S を用いた積分形では (3.18) となる。右辺の面積分にガウスの定理を用いて体積積分に変換しても求められるし、弾性運動方程式と同 じように、右辺の被積分関数について図2−4で示したと同じような表と裏の面の差を ∆x などをゼロにす る極限操作などをして体積積分の形にしても求められる。結果として、次の微分形式が求められ、通常はこ れを「連続の式」と呼ぶ: (3.19) [問題3−3] 式(3.18)を導け。またガウスの定理を用いて右辺を体積積分の形にして、式 (3.19) を導け。 [問題3−4] 図2−4の応力テンソルのように、表と裏の面を通しての質量の流れの差 ρv · n を用いて、式 (3.19)を導け。
連続の式は、式 (3.11) のようにラグランジュ微分を用いると、以下のようにも表せる: (3.20) [問題3−5] 式(3.19)から式 (3.20) を導け。 ここで式 (3.7) で示した非圧縮流体の式を、式 (3.20) から再確認する。非圧縮流体では、流体の密度が変 化しないわけだが、これは「ある物質の小片については時間変化しない」意味であり、空間的に密度が一様 でない流体では例え非圧縮であっても流れによって別の位置に移動するから、「空間に固定したある点での 密度は時間変化がないとは限らない」。つまり、非圧縮流体ではラグランジュ流で密度の時間変化がゼロで あるが、 。つまり、 (3.21) のように、速度ベクトル v の発散がゼロであることが、非圧縮流体の条件となる。
3.4
流体問題と無次元化: Non-dimensionalization of fluid dynamic equations
前節で導入した粘性流体の運動方程式であるナビア・ストークス方程式、式 (3.13) や (3.17) は、非線形 な方程式であることが本質的であることを指摘した。このような非線形方程式では足し合わせの原理が成 立しないので、大きなスケールの運動と小さなスケールの運動の間の相互作用や、非定常な複雑な時間変 化が起こることがわかっている。この方程式をそのまま解析的に解いていくことが難しいばかりか、数値的 に解いていくにも安定性などの判定が極めて難しい。多くの流体運動においては、ナビア・ストークス方 程式のすべての項が重要でなく、その中で2、3の項だけが他より圧倒的に大きくて、そのような項だけ を取り出して考えてよい場合が多い。そうすれば、解析的にも数値的にも扱いやすい。ナビア・ストーク ス方程式の各項の大きさを見積もるためには、それぞれの変数や座標、時間などを以下に示すように、方 程式を「無次元化」した形にする。その際に導入される と呼ばれる 無次元の値のみで、慣性項と粘性による応力項の大きさを見積もり、無視してもよいような項を推定して、 流体運動の様式を分類することで、より詳しい解の振る舞いを解析的にしても数値的にしても求める手順 が一般的である。 ここで、用いるのは体積力 f を省略したナビア・ストークス方程式 (3.13) と、粘性流体(ニュートン流 体)の構成則である式 (3.8) である: ρDvi Dt = ∂σji ∂xj , σij = −pδij+ 2η ˙eji= −pδij+ η !∂v i ∂xj +∂vj ∂xi " 変数や空間・時間を表す、速度と座標と時間を以下のように代表的な大きさを示す大文字の量と、ダッシュ をつけた変数の積とみなす。後者はだいたい1程度の大きさの無次元量と考える: (3.22) 圧力の項はここでは無視し、応力テンソルについても v′ iや x′iで表現した値をやはり1程度の無次元量で 表す。すると、上の2つの式からナビア・ストークス方程式は以下のようにダッシュの付いた変数などとそ れぞれのスケールを表す量で分離できる: ρV T Dv′i Dt′ = ηV L2 ! ∂2v′ i ∂x′j∂x′j + ∂2v′ j ∂x′j∂x′i " −→ ReDv ′ i Dt′ = ∂2v′ i ∂x′j∂x′j + ∂2v′ j ∂x′j∂x′i (3.23) ここで Re は以下で定義されるレイノルズ数(Reynolds number)と呼ばれる無次元の値である: (3.24)
ここで、ν = η/ρ は「動粘性率」(kinematic viscosity)と呼ばれる量であり、粘性率 η の代わりにこれを 粘性率に代えて用いる研究者もいるので注意すること。 式 (3.23) のダッシュがついている部分は、定義 (3.22) より大きさが1程度の値である。よってレイノル ズ数は、右辺の粘性項と左辺の慣性項の大きさの比を表している。2つの極端な例を考える: 1. 2. 平行な流れの中に円柱の障害物を置いた場合に、長さ L は円柱の直径で、流れの速度を V とする。密 度 ρ と粘性率 η の決まった流体を用いて、流れの速度 V をだんだん大きくしていけば、レイノルズ数を しだいに大きくしていった場合の流れの変化がわかる。図3−3はそうした実内実験の結果の一例である (Batchelor, ”An Introduction to Fluid Dynamics”, Cambridge University Press, 2000 より引用)。
R = 32 R = 55 R = 73 R = 65 R = 161 R = 102 (a) (b) (e) (c) (d) (f) 図3- 3 図3−3 (a) のように Re が1より小さいと、流れが円筒の上流側と下流側で対称の形であり、この流れ のパターンは定常状態で時間変化しない(ごくわずかはしているはずだが、測定の精度の範囲内ではそう である)。流体の速度を上げていき、(b) のように Re が1を超えると円筒の下流側の裏側に細かな渦のよ うな乱れが生じる。この部分では流れが時間とともに大きく変化している。しかし、大部分の領域での流れ はまだ定常的である。速度をさらに大きくして、(c)、そして (d) と Re が大きくなっていくと、円筒の裏側 には明瞭な渦を生じる。さらに Re > 10 になると、(e) のように円筒の裏側で渦ができたりはがれて流され て行ったりと、激しく時間変化していく。渦の大きさも最初は似ているが、次第に大小さまざまな様相の乱 れが不規則に生じ、非定常性が極めて大きくなる。究極的には、Re がもっと大きくなると、(e) のように 見かけは完全にばらばらな波の乱れ、いわゆる乱流状態へと移行する。左図にはより広範囲の流れのレイ ノルズ数による流れの変化の様子を示す。 (野球のナックルボールやサッカーの無回転シュートでは、「ボールが不規則にふらつく」が、これはボー ルの後ろに図3−3 (d) のように渦がときどきはがれるために、その反作用で逆にボールがその度に反対側 に動くのだと言われている。逆に、野球の投手が投げる直球でスピードが上がるためにはボールに十分な
回転があって、それによって乱流を防ぎ、空気抵抗が少なくなるからとも言われている。ただし、専門家に よるともっと複雑な現象らしいので、興味がある者は調べられたい。) このようにレイノルズ数は流れの様式や重要な項を見積もる場合に有効である。さらに、式(3.23)の形 式より Re さえ同じであれば、全く異なった流体運動でも全く同じように扱うことができることを示してい る。例えば、地球の北半球の中緯度から高緯度にかけてはジェット気流で蛇行した偏西風が定常的に流れて いる。中緯度で地球一周だからだいたい長さ L が 2 万キロほどである。風速 V は大きい所で 100 m/s くら いある。一方、流体の室内実験は長さが 5 m 程度で流速はせいぜい 1 m/s 程度である。よって、この二つ の現象の V L の比は 4 × 108: 1程度にもなるが、実験室の液体の密度が空気よりもはるかに大きく、さら に粘性率がずっと小さい液体を用いれば、ρ/η がこの値と逆の比に設定できることも不可能ではない。そう すれば、二つの現象のレイノルズ数は同じだから、室内実験で偏西風のパターンが再現できることになる。 流体力学では無次元化することで、時間、空間、速度のスケールが全く異なる現象でも同じ様式であると 見なして解析することが、実験的にも数値的にも可能であり、応用上で極めて重要な概念である。ここで は、ナビア・ストークス方程式に関係するレイノルズ数について詳しく解説したが、流体運動については 他にも関係する方程式をそのまま解いていくのではなく、無次元化して、そこに現れるパラメターの大き さを見積もって、その違いという視点から扱うことが有効である。例えば、マントル対流は下から暖めら れて軽くなった物質が上に上がり、上面(つまり地表面)で冷えて密度が大きくなり下に沈むという「熱 対流」(thermal convection)が重要で、このような現象では対流運動の強さを示すレーリー数(Rayleigh number)Ra や、熱伝導による熱輸送に対して実際の熱輸送の量(物質の対流によって熱が運ばれるため大 きくなる)の比であるヌッセルト数(Nusselt number)Nu などが重要となる。詳しくは、本多「マントル ダイナミクス II - 力学」、岩波講座・地球惑星科学 10, 1997 から学び始めるとよい。 [問題3−6] 以下のそれぞれの現象について、レイノルズ数がいくらだか調べ、その流れの様式との関係 を考えよ プロ野球の投手の直球、サッカーでのシュート、平泳ぎの一流選手、川を遡上するサケ、竜巻、台風、黒潮、 マグマだまりでの流れ(有珠山とハワイ等の火山では大きく異なる)、マントル対流、木星の大赤班
3.5
2次元非圧縮流体の手法: Two-dimensional incompressible fluid
ナビア・ストークス方程式は非線形なので、そのままでは非常に限られた場合しか解析的に扱えず、数値 的解法が必須となる。それでも何らかの制約条件があると、数値的な扱いをするにしてもかなり便利ない くつかの方法が考案されている。ここでは、固体地球惑星科学の流体運動として応用が広い、2次元でかつ 非圧縮な粘性流体の扱いを取り上げる。 非圧縮流体では式 (3.21) の ∇ · v = 0 となる。実際はすべての流体は圧縮されるが、多くの場合にはこの 効果が無視できる。さらに、例えばマントル対流などの問題では熱対流なので密度変化が本質的だが、圧力 効果よりも温度による熱膨張・収縮の方が効果的なので、熱輸送の方程式を別に計算することで流体運動方 程式の方はこのような仮定の下で扱える(詳しくは、前出の本多を参考のこと)。さらに、2次元の形状な らば流体運動の本質が理解できる場合も多い。2次元流体の速度場を v = (u, v, 0) とすると、式 (3.21) は つまり、速度場の2つの成分である u と v は独立でないことがわかる。簡単な洞察から、以下のような「流 れ関数(stream function)」という x と y を変数とするスカラー関数 ψ(x, y) を導入すれば、上の式は自動 的に満たされる: (3.25) つまり、流れ関数 ψ(x, y) のみを考えれば十分であることがわかる。
では、流れ関数はどんな物理的な意味があるだろうか。ψ が一定となる曲線に沿って考えると、 ψ(x, y) = const. −→ dψ = ∂ψ ∂xdx + ∂ψ ∂ydy =−vdx + udy = 0 より、 dx u = dy v (3.26) となる。つまり、図3−4 (a) のように、ψ が一定の線は速度ベクトル v = (u, v) の方向に対応するので、 この線は流体の流れの軌跡(流線 stream line とも呼ぶ)を示す、故に「流れ関数」と呼ばれる。 (a) 図3- 4
ds
dy
dx
u
v
v
(b)x
O
y
C
P
ds
n
^
(c)ψ+Δψ
ψ−Δψ
ψ
A
B
さらに流れ関数の等高線を描くと、その間隔が速度を示すことにもなる。図3−4 (b) のように、点 O と点 P の流れ関数の差を考えて、その経路 C に沿って積分すると、 ∆ψ =$ C dψ = $ C (∂ψ ∂xdx + ∂ψ ∂ydy) = $ C(udy − vdx) = $ C ! u∂y ∂s− v ∂x ∂s " ds ここで、図3−4 (b) のように線素 ds の法線単位ベクトルを ˆn とすると、 ds ds = ! ∂x ∂s, ∂y ∂s " , ˆn = ! ∂y ∂s,− ∂x ∂s " より、 (3.27) となる。右辺の被積分関数である u · ˆn は点 O と P の間の積分路 C を抜けていく単位時間あたりの流量で あることがわかる。 よって、∆ψ ずつ異なる3本の流れ関数の等高線が図3−4 (c) のようにあれば、隣接する2本の流線の 間を流れる単位時間あたりの流量は一定となる。そして、この間隔が大きい A のような部分は速度が小さ く(速度と間隔の積が一定のため)、逆に B のような細い部分は流れが大きくなり、直感的な流れの場の理 解に有用になる。 ベクトル解析の rotation(または curl や ∇×)がゼロの場合には、力学では外力はスカラーポテンシャル Uの gradient、∇U、で表せることを学んだはずである。同様に、渦がない2次元非圧縮流体の速度場は、 以下の「速度ポテンシャル」φ(x, y) で表現できる: 流れ関数 ψ(x, y) と速度場の関係式 (3.25) は、z = x + iy という複素変数を用いた複素関数 f(z) の実部を 速度ポテンシャル φ(x, y)、虚部を流れ関数 ψ(x, y)。つまり f (z)≡ φ(x, y) + iψ(x, y) とみなすと、φ と ψ は独立でなく、u と v の表現からコーシー・リーマンの微分方程式と関係で結ばれる。 そして、f(z) は「正則関数」という z で表現できる形式となり、解析的な扱いがいろいろと可能になる。こ のように複素関数による表現は、2次元流体運動についての解析的な扱いと密接な関係にある。(複素関数 に関する簡単なまとめは、蓬田 (2007) の第1章を参照のこと。)流れ関数を用いた数値計算の一例として、図 3-5 にマントル対流の様子を示す(Fowler(2005) より引用)。 上の図が温度分布に当たる等温線の分布で、下の図が流れ関数の等高線を示す。つまり下の図により速度場 が理解できる。詳細は省くが、下面の高温と上面の低温を固定した場合にいくつかの条件を変えていくと、 上昇流と下降流の間隔や形状が異なり、これらと地震波トモグラフィーなどから得られた観測結果を比較す ることで、地球内部のダイナミックな挙動が推測できる。 図3- 5 2次元非圧縮流体の運動のうち、定常状態、すなわち流れの場が時間によらず一定の場合には、流れ関数 による解析的な扱いが可能になる。定常状態ではナビア・ストークス方程式である式(3.13)の左辺の時間 微分の項がゼロになるので(定常状態でなくても非常にゆっくりと動くようなレイノルズ数が小さい場合 でも適用できる)、外力の項が無視できると となる。簡単のために粘性率 η を一定とすると、2次元非圧縮流体を考えているので、式 (3.8) の応力テン ソルの3つの成分は以下のように表せる: σxx= −p + 2η ∂u ∂x σyy = −p + 2η ∂v ∂y σxy= σyx = η ! ∂v ∂x+ ∂u ∂y " これらを上の式に代入すると、i = x の場合には 0 = ∂σxx ∂x + ∂σyx ∂y = − ∂p ∂x+ 2η ∂2u ∂x2+ η ∂2v ∂x∂y + η ∂2u ∂y2 ここで非圧縮の式 (3.21) より左辺の第3項は ∂2v/∂x∂y =−∂2u/∂x2となるので、結果的に以下のように なる: 0 = −∂p∂x+ η!∂2u ∂x2 + ∂2u ∂y2 " (3.28) 同様に、i = y については以下のような式が導ける: 0 = −∂p ∂y+ η !∂2v ∂x2 + ∂2v ∂y2 " (3.29) [問題3−7] 式 (3.29) を導け。 圧力 p の項を消去するために、式 (3.28) に −∂/∂y をかけ、式 (3.29) に ∂/∂x をかけてから、両辺を足せ ばよい。さらに u, v を流れ関数 ψ で表現すると、以下のようになる: (3.30) ここでは2次元なので、∇2= ∂2/∂x2+ ∂2/∂y2である。ラプラス方程式である ∇2ψ = 0の形は調和関数
つまり、この式を解いて ψ を求めれば、定常(または非常にゆっくり変化する)な2次元非圧縮流体の運 度が求められる。簡単な境界条件ならば、解析的な解を得る事ができる。 [問題3−8] 式 (3.30) を導け。 例えば氷河期には多くの地表に分厚い氷河で覆われており、その過重によって地殻は凹んでいた。気候 が温暖になり氷河がなくなると、過重が取れるのでアイソスタシーが成立するように凹んでいた地面が次 第に盛り上がっていく。スカンジナビア半島ではこのような原因で土地が毎年ゆっくりと上昇しており、 post-glacial reboundと呼ばれる。観測された上昇率から固体地球の粘性率を推定することができる。
その他の解析解が求められる一例として、1次元の層状の流れ場を考える(Turcotte and Schubert, ”Geo-dynamics”, 2nd ed., Cambridge Univ. Press, 2001, 第6章から引用)。図3−6 (a) のように、厚さ h の 層の x 方向の流れを考える。流れ関数を使わず、運動方程式などから導く手順を以下に述べる。
図3- 6
(a)
(b)
(c)
(d)
図3−6 (a) の点線のように長さ ℓ、厚さが無減少の δy の部分に働く力を考える。ここでは流れは1次元 的なので、x 成分の力のみを考えればよい。下流と上流の圧力を p0, p1とすると、縦の2つの面にかかる力 は (p1− p0)δyとなる。上面には −σxy(y)という単位長さあたりにせん断応力がかかり、下面には σxy(y + δ) がかかる。よって上下2面にかかる力は、δy が微小とすると、 −σxy(y)ℓ + σxy(y + δy)ℓ ≃ dσxy dy δyℓ となるので、長さ ℓ をゼロの極限を取ると、圧力による力も含めて以下のようになる: (p1− p0)δy + dσxy dy δyℓ = 0 −→ dσxy dy = − p1− p0 ℓ ≃ dp dx (3.31) ここで最後の圧力 p の微分では、上流側が下流に比べて高い(だから流れができる)から、p1> p0という 制約から符号を決めていて、dp/dx < 0 である。この1次元の流れは、図3−6 (b) のように流れ方向の位 置 x には寄らないはずなので、dp/dx は一定だと考えてよい。 一方、応力テンソル σxyは式 (3.8) で、速度場は u しかなく、かつ y のみの関数なので、 σxy= 2η ˙exy = η !∂u ∂y + ∂v ∂x " = ηdu dy (3.32) となり、よって、式 (3.31) と (3.32) をまとめると dσxy dy = η d2u dy2 = dp dx (3.33) dp/dxが一定であることから、この方程式の解は未定定数 c1, c2を用いると、以下のようになる: u = 1 2η dp dxy 2+ c 1y + c2 (3.34)2つの未定定数は、上下の壁の境界条件によって決まる。流体での境界条件としてよく用いられる”no-slip boundary condition”は、境界で壁の速度と一致するということである。図3−6 (c) のように、下面 y = h は止まっていて、上面 y = 0 は速度 u0で動いているとすると、 u = 0 at y = h, u = u0 at y = 0 を上の解に代入すると、2つの未定定数が決まり、解は u = u0 ' 1 − yh( (3.35) このような壁からの距離に速度が比例する単純な流れを”Couette flow”と呼ぶ。もし上面も静止していた ら、y = 0 での境界条件を u = 0 とするだけなので、図3−6 (d) のように解は u = 1 2η dp dxy(y− h) (3.36) のように、放物線の形の速度プロファイルになる。 [問題3−9] 解 (3.35) と (3.36) を導出せよ。 [問題3−10] 流れ関数 ψ(x, y) が重調和関数(式 (3.30))を満たすことを利用して、一次元の層状の流れ の一般解である式 (3.34) になることを示せ。また、上に示した2つの境界条件での流れ関数の形も求めよ。