ベナール対流の安定性とレイリー数

　今日は上下にそれぞれ一様な温度の床と天井があるところに閉じ込められた流体があって,元々はその高さに比例した温度変化と圧力変化を持って静止していた流体の,上下についての対流の発生,および安定性とレイリー数(Rayleigh number)との関係について述べたいと思います。

　これは,いわゆるベナール対流(Benard's convection)の問題です。

　そのために,まず流体方程式系の非圧縮近似について考察するところから始めましょう。

　実在の流体は流体速度が音速に比較して非常に小さいときは,非圧縮性流体:Ｄρ/Ｄｔ＝0 で近似することができます。

　特に重力による圧力傾度と釣り合った静力学平衡:∂Ｐ/∂ｚ＝－ρｇを満たすような密度ρ＝ρ₀ からの温度分布によるずれ:Δρだけを圧縮性として考慮し,それ以外は非圧縮性流体として扱う近似が正当化されます。

　これは,気象力学などでブジネスク(Boussinesq)近似としてよく利用される近似です。

　準備として,まずその定式化をすることにします。

　"質量保存の方程式＝連続の方程式"は,(∂ρ/∂ｔ)＋div(ρｖ)＝ 0 ですが,これをLagrange微分で表すとＤρ/Ｄｔ＋ρdivｖ＝0 です。

　Lagrange微分は,Ｄ/Ｄｔ≡∂/∂ｔ＋ｖ∇です。

　そこで,非圧縮性流体Ｄρ/Ｄｔ＝0 の近似を行なうと連続の方程式は,

　divｖ＝0 となります。

　次に,運動方程式は初めから圧力項を除く部分では非圧縮 divｖ＝0 を仮定し,外力は重力だけであるとするナビエ・ストークス方程式(Navier-Stokes equation)とします。

　すなわち,運動方程式は,圧力項以外,非圧縮のNavier Stokes方程式:

 Ｄｖ/Ｄｔ＝－∇Ｐ/ρ－ｇｋ＋ν∇²ｖ 　を出発点にします。

　ρ＝ρ₀＋Δρ,Ｐ＝Ｐ₀＋Ｐ'とおき,∇Ｐ₀＝－ρ₀ｇｋとして運動方程式の右辺の－∇Ｐ/ρを展開して２次の微小量を無視します。

　すると,－∇Ｐ/ρ＝－∇(Ｐ₀＋Ｐ')/(ρ₀＋Δρ)＝－(1/ρ₀)[１－(Δρ/ρ₀)]∇(Ｐ₀＋Ｐ')＝[1－(Δρ/ρ₀)]ｇｋ－∇Ｐ'/ρ₀です。

　ただしｇは重力の加速度,νは動粘性係数＝η/ρです。

　そこで,運動方程式は,

　Ｄｖ/Ｄｔ＝－∇Ｐ'/ρ₀－(Δρ/ρ₀)ｇｋ＋ν∇²ｖとなります。

　これがBoussinesq近似です。

　温度Ｔは密度と圧力で決まるので温度での表現も与えておきます。

 ρ₀,Ｐ₀に対応する温度をＴ₀とし,Ｔ＝Ｔ₀＋ΔＴとおくと,圧縮率をαとしてρ＝ρ₀(1－αΔＴ)となります。

　Δρ/ρ₀＝－αΔＴですから,運動方程式は,

　Ｄｖ/Ｄｔ＝－∇Ｐ'/ρ₀＋αｇΔＴｋ＋ν∇²ｖ　と変形できます。

　最後に,エネルギー方程式です。

　これは温度の方程式に直すことができて,粘性による散逸エネルギーを無視すればＤＴ/Ｄｔ＝κ∇²Ｔで与えられます。

　ここに,κ＝ｋ/(ρＣ_p)は温度拡散係数です。

　ｋは熱伝導係数であり,Ｃ_pは定圧比熱です。

　以下では,ｖ＝(ｕ,ｖ,ｗ)として話を進めます。

　まず最初の設定としては領域全体でｕ＝ｖ＝ｗ＝0 とします。

　そして,温度は床ｚ＝0 ではＴ＝Ｔ₀,天井ｚ＝ｈではＴ＝Ｔ₁で時間的に変わらないとします。この境界条件を満たす定常熱伝導の方程式∇²Ｔ＝0 の解はＴ(ｚ)＝Ｔ₀－βｚ＝Ｔ₀＋{(Ｔ₁－Ｔ₀)/ｈ}ｚですね。

　ここで,β＞0 ,すなわちＴ₀＞Ｔ₁とします。

　何故なら,

　さもないと地球上の重力場の中では対流が起きないからです。

　　    ベナール対流の模式図↓

 (PDF:熱対流現象(山中透:京大)から勝手に引用)

　そればかりでなく|β|＝－βの値が,いわゆる断熱温度減率(＝湿度によって変化します)より大きくないと対流は起きません。(例えばランダウ著「流体力学１」(東京図書)参照） 

　そして,圧力はＰ(ｚ)＝－ρｇｚ＋const.とします。

　ｖ＝(ｕ,ｖ,ｗ)を,改めてｕ＝ｖ＝ｗ＝0 の平衡状態からのずれ:摂動流速と見なし,温度もＴ≡Ｔ(ｚ)＋θ≡Ｔ－βｚ＋θ,圧力もＰ≡Ｐ(ｚ)＋Ｐ'として摂動温度θ,摂動圧力Ｐ'を取ります。

　連続方程式は,∂ｕ/∂ｘ＋∂ｖ/∂ｙ＋∂ｗ/∂ｚ＝0 です。

　またBoussinesq近似での運動方程式は,

　Ｄｕ/Ｄｔ＝－(∂Ｐ'/∂ｘ)/ρ₀＋ν∇²ｕ,

　Ｄｖ/Ｄｔ＝－(∂Ｐ'/∂ｙ)/ρ₀＋ν∇²ｖ,

　Ｄｗ/Ｄｔ＝－(∂Ｐ'/∂ｚ)/ρ₀＋αｇθ＋ν∇²ｗ

　となります。 

　温度方程式は,ＤＴ/Ｄｔ＝κ∇²Ｔですが,ＤＴ/Ｄｔの中の移流項ｖ∇Ｔは主流を取って,ｗ∂Ｔ/∂ｚ＝－βｗのみですから, 

　∂θ/∂ｔ＝βｗ＋κ∇²θ　となります。 

　ここで,運動方程式からｕとｖを消去するために,ｖ＝(ｕ,ｖ,ｗ)の"回転＝渦"(vortex):ω＝rotｖを取ります。

　ζ≡∂ｖ/∂ｘ－∂ｕ/∂ｙとすると,運動方程式から

　∂ζ/∂ｔ＝ν∇²ζが得られます。

　また,(∂/∂ｔ)∇²ｗ＝αｇ(∇²－∂²/∂ｚ²)θ＋ν∇⁴ｗ

なる式も得られます。

　こうして基本方程式系としては,ｗ,θ,ζの３つを未知数とする３つの微分方程式が得られたことになります。

　さらに,ｘ,ｙ,ｔの独立変数をｚから分離するため,変数分離形を仮定して,ｗ≡Ｗ(ｚ)exp[i(ｋ_xｘ＋ｋ_yｙ)＋ｐｔ],θ≡Θ(ｚ)exp[i(ｋ_xｘ＋ｋ_yｙ)＋ｐｔ],ζ≡Ζ(ｚ)exp[i(ｋ_xｘ＋ｋ_yｙ)＋ｐｔ]　と置きます。

　これは方程式系において,∂/∂ｔ＝ｐ,∇²－∂²/∂ｚ²＝－ｋ²≡－(ｋ_x²＋ｋ_y²),∇²＝ｄ²/ｄｚ²－ｋ²と置き換えることに相当します。 

　すると,

　ｐ(ｄ²/ｄｚ²－ｋ²)Ｗ＝－αｇｋ²Θ＋ν(ｄ²/ｄｚ²－ｋ²)²Ｗ,

　ｐΘ＝βＷ＋κ(ｄ²/ｄｚ²－ｋ²)Θ,

　ｐΖ＝ν(ｄ²/ｄｚ²－ｋ²)Ζ

　となります。

　これは３つめの方程式とは別に,前の２つのＷとΘだけで独立しているのでこれら２つだけでも解ける形になっています。

　そして,床と天井での境界条件はｚ＝0 とｚ＝ｈでΘ＝Ｗ＝Ζ＝0 ,ｄＷ/ｄｚ＝0 です。

　ここで,これらの方程式を無次元化するために単位を[Ｌ]＝ｈ,[Ｔ]＝ｈ²/νに取ることにして,ａ≡ｋｈ,σ≡ｐｈ²/νとします。

　ｘ,ｙ,ｚの単位もｈとした後に,Ｄ≡ｄ/ｄｚと定義し,温度拡散のプラントル数(Prandtl number)をＰｒ≡ν/κとします。

　すると,３つのうち前の２つの方程式は,

　(Ｄ²－ａ²)(Ｄ²－ａ²－σ)Ｗ＝(αｇｈ²/ν)ａ²Θ,および

　(Ｄ²－ａ²－Ｐｒ・σ)Θ＝－(βＨ²/κ)Ｗ　となります。

　そして境界条件は,ｚ＝0 とｚ＝1 においてΘ＝Ｗ＝Ζ＝0 ,

　ＤＷ＝0 です。

　ただし,ＷとΘは通常の単位であり無次元化していません。

　Θを消去すると,

　(Ｄ²－ａ²)(Ｄ²－ａ²－σ)(Ｄ²－ａ²－Ｐｒ・σ)Ｗ＝－Ｒａ²Ｗ

　となります。

　Ｒはレイリー数(Rayleigh number)で,これはＲ≡ｇαβｈ⁴/(κν)で定義される無次元数であり方程式を支配する特徴的な数です。

　安定性を調べるためには,σ＝0 (つまりｐ＝0)の特別なケースを調べればよいと思われます。

　何故なら,時間的に減衰するか増幅するかの限界は,exp(ｐｔ)という因子でのｐ＝0 を臨界点として生じるからです。

　そして,その境界におけるRayleigh数に最初に達したときに流れは初めて不安定になると考えられるからです。

そして対称性を利用するために,ｚの原点を移動して,ｚ＝±1/2を境界に取ります。

基本となる微分方程式は(Ｄ²－ａ²)³Ｗ＝－Ｒａ²Ｗであり,

境界条件はＷ＝ＤＷ＝(Ｄ²－ａ²)²Ｗ＝0 atｚ＝±1/2です。

この方程式の特性方程式は,

Ｗ≡ｅ^qzとして,(ｑ²－ａ²)³＝－Ｒａ²です。

Ｒａ²＝τ³ａ⁶と置いて解くと,

ｑ²＝－ａ²(τ－1),ｑ²＝ａ²[1＋(1/2)τ(1±ｉ√3)],

または平方根を取って,±iｑ₀,±ｑ,±ｑ^* (ｑ₀＝ａ(τ－1)^1/2 etc.)

ですね。 

　あとの詳細は簡単なので本質的な部分以外は省略しますが,解は偶関数解と奇関数解の２種類に分かれ,

　偶関数解としては,

　Ｗ＝Ａ₀cos(ｑ₀ｚ)＋Ａcosh(ｑｚ)＋Ａ^*cosh(ｑ^*ｚ),

 奇関数解としては,

　Ｗ＝Ａ₀sin(ｑ₀ｚ)＋Ａsinh(ｑｚ)＋Ａ^*sinh(ｑ^*ｚ)

　が得られます。

　後は,境界ｚ＝±1/2でＷ­＝ＤＷ＝(Ｄ²－ａ²)²Ｗ＝0 である,という条件からＡ₀,Ａ,Ａ^*を決めればいいわけです。

　これらが自明でない解を持つための条件は境界条件を満足する係数ベクトル:Ａ₀,Ａ,Ａ^*の係数行列の行列式がゼロでなければならないという1つの固有値方程式になります。

　例えば偶関数解ではこの行列式＝0 の条件は,Im[(i＋√3)ｑtanh(ｑ/2)]＋ｑ₀tanh(ｑ₀/2)＝0 となります。

　こうした固有値問題を満足するＲを求めると,流れが不安定になり熱伝導だけでは熱を運ぶことが不可能となって物質で熱を運ぶ対流が激しくなる臨界点が求まり,最後には乱流へと移行することになる境界の臨界レイリー数が求められます。

　私がかつて数値的に解いた偶関数の臨界レイリー数Ｒはたとえばａ＝3.0でＲ＝1710 etc.であり文献と一致しています。

　つまり,ａ＝3.0ではこのＲを境として不安定流となり対流が激化していくことになるわけですね。

　参考文献：S.Chandrasekhar 「Hydromagnetic and Hydromagnetic stability」(Dover Books),モーニン,ヤグロム(A.S.Monin&A.M.Yaglom)著(山田豊一訳)「統計流体力学1」(総合図書),

(追加):神部勉,Ｐ.Ｇ.ドレイジン「流体力学の安定性と乱流」(東京大学出版会)