« 磁石と鉄球の引力(誤解答の例) | トップページ | 数論の演習問題 »

2006年9月25日 (月)

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

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

 

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

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

 

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

 

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

 

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

 

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

 "質量保存の方程式=連続の方程式"は,(∂ρ/∂t)+div(ρ)= 0 ですが,これをLagrange微分で表すとDρ/Dt+ρdiv=0 です。

 

 Lagrange微分は,D/Dt≡∂/∂t+∇です。

 

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

 div=0 となります。

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

 

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

 /Dt=-∇P/ρ-g+ν∇2v  を出発点にします。

 

 ρ=ρ0+Δρ,P=P0+P'とおき,∇P0=-ρ0として運動方程式の右辺の-∇P/ρを展開して2次の微小量を無視します。

 

 すると,-∇P/ρ=-∇(P0+P')/(ρ0+Δρ)=-(1/ρ0)[1-(Δρ/ρ0)]∇(P0+P')=[1-(Δρ/ρ0)]g-∇P'/ρ0です。

 

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

 そこで,運動方程式は,

 D/Dt=-∇P'/ρ0(Δρ/ρ0)g+ν∇2となります。

 

 これがBoussinesq近似です。

 

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

 

 ρ0,P0に対応する温度をT0とし,T=T0+ΔTとおくと,圧縮率をαとしてρ=ρ0(1-αΔT)となります。

 

 Δρ/ρ0=-αΔTですから,運動方程式は,

 D/Dt=-∇P'/ρ0+αgΔT+ν∇2v と変形できます。

 最後に,エネルギー方程式です。
 これは温度の方程式に直すことができて,粘性による散逸エネルギーを無視すればDT/Dt=κ∇2Tで与えられます。
 
 ここに,κ=k/(ρCp)は温度拡散係数です。
 kは熱伝導係数であり,Cpは定圧比熱です。

 以下では,(u,v,w)として話を進めます。

 

 まず最初の設定としては領域全体でu=v=w=0 とします。

 

 そして,温度は床z=0 ではT=T0,天井z=hではT=T1で時間的に変わらないとします。この境界条件を満たす定常熱伝導の方程式∇2T=0 の解はT(z)=T0-βz=T0{(T1-T0)/h}zですね。

 

 ここで,β>0 ,すなわちT0>T1とします。

 何故なら,

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

 

      ベナール対流の模式図↓

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

 

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

 そして,圧力はP(z)=-ρgz+const.とします。

 

 =(u,v,w)を,改めてu=v=w=0 の平衡状態からのずれ:摂動流速と見なし,温度もT≡T(z)+θ≡T-βz+θ,圧力もP≡P(z)+P'として摂動温度θ,摂動圧力P'を取ります。

 

 連続方程式は,∂u/∂x+∂v/∂y+∂w/∂z=0 です。

 

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

 

 Du/Dt=-(∂P'/∂x)/ρ0+ν∇2,

 Dv/Dt=-(∂P'/∂y)/ρ0+ν∇2,

 Dw/Dt=-(∂P'/∂z)/ρ0+αgθ+ν∇2

 

 となります。 

 温度方程式は,DT/Dt=κ∇2Tですが,DT/Dtの中の移流項∇Tは主流を取って,w∂T/∂z=-βwのみですから, 

 ∂θ/∂t=βw+κ∇2θ となります。 

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

 

 ζ≡∂v/∂x-∂u/∂yとすると,運動方程式から

 ∂ζ/∂t=ν∇2ζが得られます。

 

 また,(∂/∂t)∇2w=αg(∇2-∂2/∂z2)θ+ν∇4

なる式も得られます。

 

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

 さらに,x,y,tの独立変数をzから分離するため,変数分離形を仮定して,w≡W(z)exp[i(kxx+ky)+pt],θ≡Θ(z)exp[i(kxx+ky)+pt],ζ≡Ζ(z)exp[i(kxx+ky)+pt] と置きます。

 

 これは方程式系において,∂/∂t=p,∇2-∂2/∂z2=-k2≡-(kx2+ky2),∇2=d2/dz2-k2と置き換えることに相当します。 

 すると,

 

 p(d2/dz2-k2)W=-αgk2Θ+ν(d2/dz2-k2)2W,

 pΘ=βW+κ(d2/dz2-k2)Θ,

 pΖ=ν(d2/dz2-k2

 

 となります。

 

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

 そして,床と天井での境界条件はz=0 とz=hでΘ=W=Ζ=0 ,dW/dz=0 です。

 

 ここで,これらの方程式を無次元化するために単位を[L]=h,[T]=h2/νに取ることにして,a≡kh,σ≡ph2/νとします。

 

 x,y,zの単位もhとした後に,D≡d/dzと定義し,温度拡散のプラントル数(Prandtl number)をPr≡ν/κとします。

 すると,3つのうち前の2つの方程式は,

 (D2-a2)(D2-a2-σ)W=(αgh2/ν)a2Θ,および

 (D2-a2-Pr・σ)Θ=-(βH2/κ)W となります。

 

 そして境界条件は,z=0 とz=1 においてΘ=W=Ζ=0 ,

 DW=0 です。

 

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

 Θを消去すると,

 (D2-a2)(D2-a2-σ)(D2-a2-Pr・σ)W=-Ra2

 となります。

 

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

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

 

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

 

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

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

 

基本となる微分方程式は(D2-a2)3W=-Ra2Wであり,

境界条件はW=DW=(D2-a2)2W=0 atz=±1/2です。

 

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

W≡eqzとして,(q2-a2)3=-Ra2です。

 

Ra2=τ36と置いて解くと,

2=-a2(τ-1),q2=a2[1+(1/2)τ(1±i√3)],

または平方根を取って,±iq0,±q,±q* (q0=a(τ-1)1/2 etc.)

ですね。 

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

 

 偶関数解としては,

 W=A0cos(q0z)+Acosh(qz)+A*cosh(q*z),

 奇関数解としては,

 W=A0sin(q0z)+Asinh(qz)+A*sinh(q*z)

 が得られます。

 後は,境界z=±1/2でW­=DW=(D2-a2)2W=0 である,という条件から0,A,A*を決めればいいわけです。

 

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

 例えば偶関数解ではこの行列式=0 の条件は,Im[(i+√3)qtanh(q/2)]+0tanh(q0/2)=0 となります。

 

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

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

 

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

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

 

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

|

« 磁石と鉄球の引力(誤解答の例) | トップページ | 数論の演習問題 »

308. 微分方程式」カテゴリの記事

108. 連続体・流体力学」カテゴリの記事

コメント

コメントを書く



(ウェブ上には掲載しません)




トラックバック


この記事へのトラックバック一覧です: ベナール対流の安定性とレイリー数:

« 磁石と鉄球の引力(誤解答の例) | トップページ | 数論の演習問題 »