円筒座標系の熱流体方程式と熱ロスビー波

　大気大循環に伴う波動発生のメカニズムは、例えば北極と見做した中心を氷等で冷やした回転二重円筒間の流れ(熱ロスビー波)で模擬される事が知られています。そこで今回はこの実験に倣って、半径方向に温度勾配を与えられながら、一定の回転角速度で回転し続ける二重円筒形状容器を考え、その内部に充填された試験流体の熱流動の数値シミュレーションを独自に開発して波動の再現を試みました。

　解析に用いる基礎方程式は、次のような連立非線形偏微分方程式です。

【円筒座標系で表現される強制対流のNavier-Stokes方程式】
　（∂ｖ/∂ｔ）＋（ｖ・∇）ｖ＋２ω×ｖ＝－∇ｐ＋（１/Ｒｅ）∇²ｖ＋（Ｇｒ/Ｒｅ²）Ｔδ_i3

【円筒座標系で表現される熱エネルギー方程式】
　（∂Ｔ/∂ｔ）＋（ｖ・∇）Ｔ＋２ω×ｖ＝（１/Ｐｒ・Ｒｅ）∇²Ｔ

【微分演算子とベクトル】
　∇＝（∂/∂ｒ）＋（∂/ｒ∂θ）＋（∂/∂ｚ），ｖ＝（ｖ_ｒ，ｖ_θ，ｖ_ｚ）

　ここで、ｖ_ｒ、ｖ_θ、ｖ_ｚは各々半径方向、円周方向、鉛直方向の速度成分であり、ｐ、Ｔは各々圧力、温度である。また、Ｒｅ、Ｇｒ、Ｐｒは各々Reynolds数、Grashof数、Prandtl数である。

　今回想定する実験装置(数値モデル)の概要は次のようなイメージです。

　

　実験装置は半径r_０及びr_１の二重円筒で構成されており、内側の円筒を冷源T_ｃとし、外側の容器内に試験流体が充填されるものとします。更に外側の容器の周囲には一様な熱源T_ｈが分布しているものとします。

　そして、この実験装置全体を一定角速度ω_０で反時計回りに回転させます。

　ここで、寸法はΔr= r_１－r_０を基準として、r_０＝Δr、r_１＝2Δr,z＝Δr/2に設定しました。 従って、Reynolds数は試験流体の動粘度係数νを用いてRe=ω_０(Δr)^２ /νで表され、Grashof数は体膨張率β_０、重力加速度gを用いる事により、Gr＝gβ_０(T_ｈ－T_ｃ)×(Δr)^３ /ν^２で表されます。

　続いて、初期条件・境界条件の概要は次のような形です。

　壁面摩擦を受けない場合v_θは単純にv_θ＝ｒωに従うため、最外周で極大となります。しかし、実際には壁面摩擦を受けるので内壁面および最外周壁面上ではθ方向の速度は0にならなければいけません。初期状態においてはこの境界条件を満足する発達した流れを想定し、v_θ(r)はr方向に放物状の速度分布を成すものと仮定します。

　また、温度に関する境界条件は、内側円筒に相当するr≦r_０は冷源であるため T=T_ｃ、最外周壁面に相当するr=r_１は熱源であるため T=T_ｈです。初期条件は、r_０≦r≦r_０+Δr/2では試験流体が冷却されているため T=T_０とし、r_０+Δr/2≦r≦r_１では試験流体が加熱されているためT=T_０＋ΔTとします。

　今回はGr＝1.0×10^２,Re＝3.0×10^２,Pr＝0.71の条件を設定し,３次元熱流動の数値シミュレーションを行いました。その結果は・・・。

　こちらは、円筒の高さの半分に相当するZ=Δr/4における水平面上の温度と流れの様子を解析した結果です。波数5の熱ロスビー波が形成されているのが分かります。波が中心に向かって盛り上がる位相(リッジ)の右側では時計回り、左側では反時計周りの渦が形成されています。一方、中心から遠ざかる位相(トラフ)ではこれとは反対の構造になっています。

　続いてこちらは、この水平面上における鉛直渦度ζと流れの様子を解析した結果です。

　正渦度域(ζ>0)と負渦度域(ζ<0)が交互に形成されている事が分かります。従って、トラフの右側(リッジの左側)では正渦度循環、リッジの右側(トラフの左側)では負渦度循環が形成される現象が再現されました。