745 日本機械学会論文集 (B 編 ) 78 巻 788 号 (-4) 原著論文 No.-JBR- * akas FRSAWA *, Kotaro MAKINO, Sator YAMAMOO and Kazro NAKAHASHI * Det. of Comter and Matematal Senes, ook nversty, 6-6- Aramak Aza Aoba, Aoba-k, Senda-s, Myag, 98-8579 Jaan Preondtonng metod develoed by te ators was oled wt Bldng-Cbe Metod(BCM) w an smlate flows arond a omlex-saed body. In ts stdy, ts metod s aled to te smlaton of termal onvetve flows arond an arbtrary-saed body. Frst, fored and natral onvetons arond a ylnder are allated and te reslts are omared wt te exerments and te revos nmeral reslts to valdate te allated reslts. Seond, natral onveton arond a fn-saed body s allated and te avalablty of ts metod to ratal roblems nvolvng a termal onvetve flow s dsssed. Key Words : Nmeral Smlaton, ermal Conveton, Preondtonng Metod, BCM 古澤卓 *, 牧野幸太郎 *, 山本悟 * *3, 中橋和博 Smlaton of ermal Conveton Based on Preondtonng Metod and BCM. 緒言 PC やスマートフォンなど I 機器の高速化やクラウド志向によるデータセンターの大規模化に伴い,CP やメモリなどから発生する熱の増大が表面化している. これらエネルギー機器や I 機器から発生する熱を効率的に制御するためには, その温度差に伴って発生する熱対流を正確に把握することが重要である. しかしながら, これらの機器は概して複雑な形状を有し, かつ発電設備などマクロスケールから半導体チップなどマイクロスケールまで広範囲の空間スケールに亘っている. したがって, これらの熱対流をほぼ同じ手間と時間で解析できる CFD 手法を開発すれば極めて有用なツールになり得る. 熱対流の計算手法として,Savlle ら ( ) の境界層方程式に基づく解法や,Ken ら ( ) のナビエ ストークス方程式を解く手法などこれまでいくつかの手法がすでに提案されている. それらに共通している問題は, 非圧縮性流れを仮定した支配方程式を解いているため密度変化に伴う浮力を直接考慮することができない点であり, これをいわゆる Bssnesq 近似により近似的に評価している. しかしながら Gray ら (3 ) は, 空気では温度差 8.6 以下, 水では. 5 以下でなければこの近似手法が正確ではないことを報告している. したがって, 温度差がより大きな熱対流問題を解析するためには, 本来の密度変化が考慮できる圧縮性流れを仮定した数値解法が必要であることが示唆される. 一方, 圧縮性流れの数値解法は自然対流などの遅い熱対流にそのまま適用すると, 解の硬直 (Stffness) が生じて解の収束性が極端に悪くなることが知られている.Co ら (4) はこれを解決するため前処理法を提案した. 前処理法は, 遅い流れでは擬似密度を導入した非圧縮性ナビエ ストークス方程式に, 速い流れでは圧縮性ナビエ ストークス方程式に支配方程式が切り替わる方法である.Yamamoto ら ( 5 ) は, 前処理法に基づく圧縮性熱対流の数値解法を提案して,Bssnesq 近似を用いずに様々な物質の熱対流問題を計算する手法を開発した. 本数値解法 ( 6 ) -( 8 はその後, 凝縮流体や超臨界流体の数値解法 ) に拡張されている. しかしながら, 一般曲線座標系の構造格子を用いており, 複雑な形状の格子生成が容易ではないという問題点があった. 本論文の目的は, 著者らがこれまでに開発してきた前処理法に基づく任意物質のための熱対流の数値解法に, 中橋らが提案した Bldng-Cbe Metod(BCM) ( 9 ) を組み込むことにより, 任意形状物体周りの熱対流が計算できる数値解法を構築して, その妥当性を検証することである. また, 本数値解法では熱物性値を PROPAH () によって算出しており, * 原稿受付 年 月 5 日 * 正員, 東北大学大学院情報科学研究科 ( 98-8579 宮城県仙台市青葉区荒巻字青葉 6-6-) * 東北大学大学院情報科学研究科 *3 正員, フェロー, 東北大学大学院工学研究科 E-mal frsawa@aero.me.took.a.j 39
様々な物質の熱対流をそれぞれの物質の熱物性値に基づいて計算することができる.BCM は直交格子に基づく差分計算手法である. 計算領域は Cbe の集合体からなり, 物体近傍ではより細かな Cbe, 物体から離れるに従って Cbe を粗くする. 各 Cbe は三次元の場合, さらに n 3 格子点数 ( ただし,n は Cbe の辺を構成する格子点数 ) からになり, Cbe 毎は完全に並列計算できる.. 基礎方程式及び数値解法前処理を施した二次元圧縮性ナビエ ストークス方程式は, 一般曲線座標系で次のように定義される. = =, Γ S F F t Q v () ただし,,, = = e x x J F J Q = =, g g J S x x J F s s k jk j j j v κ τ τ τ,,,,, e, τ j, κ ならびに J は, それぞれ, 圧力, 物理速度成分, 温度, 反変速度成分, 密度, 全内部エネルギー, 粘性応力テンソル, 熱伝導率, ならびに変換のヤコビアンである.S は浮力項であり, s, g は基準密度, 重力加速度である. また, 前処理行列 Γ は次のように定義される. Γ = () ただし, = (e) / であり, Τ, ならびに はそれぞれ密度の温度偏微分, の温度偏微分, の圧力偏微分である. また, は次式のように定義される. r = (3) r は Wess ら () の定義より次のように選択される. = r ε < < < < ε ε (4) ここで, は音速であり,ε は極めて小さな値, = である. 空間差分, 時間積分には Yamamoto (6) によって提案された前処理型 Roe スキームと 次精度 MSCL スキーム, 前処理型 L-SGS スキームを用いた. 粘性項の計算には 次精度中心差分を用いた. 密度, 音速, 熱伝導率や分子粘性係数の値などの熱物性値は,PROPAH () に定義された関数により温度, 圧力を引数として計算される. 例えば水の場合には IAPWS IF97 で定義された値, 二酸化炭素は IPAC に基づいた Vral 型状態方程式に基づき計算される. また, PROPAH では物質ごとのライブラリファイルとして提供されているために物質の変更が容易に可能である. 各イタレーションごとに格子点の圧力および温度を用いて密度, 音速, 熱伝導率や分子粘性係数の値が計算される. 746 4
747 3. 計算結果 3 計算格子計算格子には中橋らが提案した Bldng-Cbe Metod(BCM) ( 9) を用いる.BCM はマルチブロック構造を用いた直交格子法である.BCM では計算空間が多数の Cbe に分割され, それぞれの Cbe で独立に計算が行われる.Cbe のサイズを局所的に変更することにより, 等間隔直交格子の単純さを保ちつつ局所的に解像度を上げることを可能にしている. 円柱周りの Cbe 配置例を図 に示す. また, ソルバの単純性を保つため, 図 に示すように物体表面は多数の格子点で階段状に表現されている. それぞれの Cbe は等間隔直交格子を持つため,Cbe ごとの計算負荷はほぼ等しく並列化が容易である. また, 各 Cbe は隣接 Cbe とのオーバラップ領域を持ち, 各計算ステップで未知変数のデータやりとりを行う. 本論文では二次元 BCM にて格子を生成して, 前処理法に基づく数値解法を用いて計算を行った. Fg. Semat of be dstrbton on BCM. Fg. Grds near te ylnder. 3 円柱周りの二次元強制対流本数値解法の強制対流における妥当性を検証するため, まず初めに円柱周りの二次元強制対流問題を計算した. 計算条件として, 対流が定常になるレイノルズ数が ならびに 4 の場合の低レイノルズ数流れについて計算し, 得られた双子渦の長さを既存の計算結果と比較した. 計算格子として 個の Cbe を用い, 総計算格子数は約.8 万点である. また, 円柱後方で Cbe を細分化している. 流体は の水を仮定し, 各イタレーションごとに格子点の圧力および温度を用いて密度, 音速, 熱伝導率や分子粘性係数の値が計算される. また, 円柱の直径を mm とした. 図 3(a),(b) にそれぞれ計算により得られた流線分布を示す. 円柱後方にそれぞれ面対称渦が捕獲されており, ノイノルズ数が 4 の場合にはより長い双子渦が形成されている.Re=4 の場合の双子渦の長さを既存の計算結果と比較した結果を表 に示す. 本数値解法による長さはおおよそ他の結果と一致しており, 強制対流において本数値解法によって妥当な解が得られることを示している. (a) Re= Fg.3 Stream lnes near te ylnder. (b) Re=4 4
748 able Vortex lengt bend te ylnder (Re=4). L/D seng et al. (). Ye et al. (3).7 Present Metod.7 3 3 円柱周りの二次元自然対流 本数値解法の自然対流における妥当性を検証するため, 一定温度の雰囲気中に置かれた高温円柱周りの自然対流の計 算を行い, レイリー数が 5 の場合の温度分布を Ken らの実験結果 () と比較した. 流体は水蒸気を仮定し, 雰囲気温 度を 4K, 円柱表面温度を 45K とした. 計算格子として 368 個および 78 個の Cbe を用い, 総計算格子数はそれぞ れ約 9.8 万点および約 万点である.368 個の Cbe を用いた計算によって得られた温度の等値線図を Ken の実験結 果と共に図 4 に示す. 円柱近傍で熱せられた流体が浮力の影響で円柱表面に沿って上昇しており, 円柱上部に高温領域 プルームが形成されている. これらは計算結果でも同様の様子が再現されている. それぞれのグリッドでの円柱から法線方向の温度分布の比較を図 5 に示す. ただし, φ = ( ) ( w ), Y * 4 = YRa D である. ここで,, w, Y, D はそれぞれ雰囲気温度, 円柱壁面温度, 円柱からの距離, 円柱直径である. また, 鉛直下向きを とした. 円柱近傍において,78 個の Cbe を用いた計算では実験とのずれがあるものの,368 個の Cbe を用いた計算では実験とよく一致しており, 十分に Cbe を細分化することで円柱周りの自然対流の温度分布を本数値解法で再現できることを示した. 次に, 異なる円柱直径におけるレイリー数が ~ 6 での加熱円柱周りの自然対流の計算を 368 個の Cbe を用いて行い, それぞれの結果を比較した. 図 6 にそれぞれの温度分布を示す.Ra= の場合には円柱直径に対して温度境界層は厚くなっている. 一方で, レイリー数が大きくなるにつれて温度境界層は薄くなり, 円柱上部の高温領域のプルームも細くなった. これらの結果は Ken らの結果と同様の傾向を示している. Fg.4 Comarson of exermental and nmeral soterms at Ra= 5 (left: resent metod, rgt: exerment by Ken () ). 4
749..9.8 Present Metod by 368bes (9 degree ) Present Metod by 78bes (9 degree ) Exerments by Ken (9 degree) Present Metod by 368bes (8 degree ) Present Metod by 78bes (8 degree ) Exerments by Ken (8 degree).7 φ.6.5.4.3... 3 4 5 6 7 Y * Fg.5 Comarson of exermental and nmeral temeratre dstrbtons. (a) Ra= (b) Ra= () Ra= 4 (d) Ra= 6 Fg.6 emeratre ontors near te ylnder. 43
75 3 4 複雑形状物体周りの二次元自然対流 BCM では複雑な形状でも容易に格子生成が可能である. さらに本研究では熱物性値算出に PROPAH を用いているために, 様々な物質の熱対流の数値計算が可能である. ここでは任意形状の計算例として図 7 に示すような凸凹を持つフィン形状周りの熱対流を計算した. 凹凸間隔を 4.6mm, 外寸 mm とし, 作動流体は水, 二酸化炭素および窒素の雰囲気温度を 93.5K, 壁面温度を 35.5K とした場合の計算を行った. 計算格子として 495 個の Cbe を用い総計算格子数は約.6 万点である. 図 7(a)(b) に水の自然対流における温度分布, 密度分布を示す. 凹部の内側にも温度分布があり, 円柱の結果と同様に上部に高温のプルームが発生している. しかし, 凹凸の影響によってプルームがやや右側に偏るという興味深い結果となった. また, 密度分布も温度分布と同様に壁面近傍で変化しており, 凹部の内側に分布が見られる. ここでの作動流体は水であるが, 従来の Bssnesq 近似に基づく自然対流の計算では水は非圧縮性流体と仮定されて計算される. 一方, 本数値解法では水の密度は厳密に考慮して計算するので, 図 7(b) で示されているように温度に伴う水の微小な密度差を捕獲することができ, 実際の熱物理に忠実に計算される. 図 8(a) にはフィン形状先端近傍における速度ベクトルを可視化する. 温度分布で示されるプルームに沿う形で熱対流が発生していることが示されている. 図 8(b) にはさらに先端部分を拡大して可視化した図を示す. 本解法では BCM を用いていることから物体境界は階段状になっている.CG の都合で図中の物体はなめらかな直線になっているが, 実際には階段状の境界である. また, 図 9 に外寸で無次元化した物体近傍の温度分布を示す. 物体近傍から滑らかな分布を持ち, 今回計算した自然対流問題では速度ベクトル, 温度分布共においてその影響は確認できなかった. 図 (a),(b) に同じ条件における二酸化炭素および窒素の温度分布, 密度分布を示す. 水の場合と比較すると熱伝導が支配的となっている. また, それぞれの物質において密度が異なる. 本数値解法を用いることで任意形状における様々な物質の熱対流が計算可能である. しかし, それぞれの物質においてさらなる検証が必要である. (a) emeratre ontors [K] (b) Densty ontors [kg/m 3 ] Fg.7 Natral onveton arond te fn sae. (a) Near te to regon (b) Fosed on te to Fg.8 Veloty vetors near te to regon of fn sae 44
75 Fg.9 emeratre dstrbton near te fn bondary. emeratre ontors [K] Densty ontors[kg/m 3 ] (a) Carbon doxde emeratre ontors [K] Densty ontors[kg/m 3 ] (b) Ntrogen Fg. Natral onveton of arbon doxde and ntrogen arond te fn sae. 45
75 4. 結論前処理法と BCM に基づく数値解法を構築して熱対流問題の計算に応用しその妥当性を検証した. まず, 円柱周りの低レイノルズ数流れの計算では, 計算により得られた双子渦の長さを既知の計算結果と比較したところ, ほぼ同様の値が得られることが示された. 次に, 円柱周りにおける水蒸気の自然対流を計算して得られた温度分布を実験値と比較したところ, 良い一致が示された. さらに, レイリー数を変化させて計算したところ, 既知の計算結果と同程度の温度分布が得られることがわかり, 広範囲のレイリー数に対して本数値解法が適用できることが示された. 最後に, 任意形状物体への応用として, フィン形状周りにおける水, 二酸化炭素, 窒素の自然対流を計算した. それぞれの物質の熱物性値の違いにより温度分布, 密度分布が異なり, 本研究で提案した数値解法が様々な物質における任意形状物体周りの自然対流計算に適用できることを示した. 謝 辞 本研究の一部は, 科研費基盤 S( 課題番号 :68, 研究代表者 : 中橋和博 ) ならびに科研費特別研究員奨励費 ( 課題番号 :576, 研究代表者 : 古澤卓 ) により実施された. 文 献 () Savlle, D.A. and Crll, S.W., Lamnar Free Conveton n Bondary Layers near Horzontal Cylnders and Vertal Axsymmetr Bodes, e Jornal of Fld Means, Vol.9, (967),.39-399. () Ken,.H. and Goldsten, R.J., Nmeral Solton to te Naver-Stokes Eqatons for Lamnar Natral Conveton abot A Horzontal Isotermal Crlar Cylnder, Internatonal Jornal of Heat and Mass ransfer, Vol.3, (98),.97-979. (3) Gray, D.D. and Goggn, A., e Valdty of te Bossnesq Aroxmaton for Lqd and Gases, Internatonal Jornal of Heat and Mass ransfer, Vol.9, (976),.545-55. (4) Co, Y.-H. and Merkle, C.L., e Alaton of Preondtonng n Vsos Flows, Jornal of Comtatonal Pyss, Vol.5, (993),.7-3. (5) Yamamoto, S., Nyama, D., and Sn. B.-R., A Nmeral Metod for Natral Conveton and Heat Condton arond and n A Horzontal Crlar Pe, Internatonal Jornal of Heat and Mass ransfer, Vol.47, (4),.578-579. (6) Yamamoto, S., Preondtonng Metod for Condensate Fld and Sold Colng Problems n General Crvlnear Coordnates, Jornal of Comtatonal Pyss, Vol.7, (5),.4-6. (7) Yamamoto, S., Frsawa,. and Matszawa, R., Nmeral Smlaton of Serrtal Carbon Doxde Flows aross Crtal Pont, Internatonal Jornal of Heat and Mass ransfer, Vol.54, No.4(),.774-78. (8) 古澤卓, 山本悟, 相変化を伴う超臨界水熱合成の数値シミュレーション, 日本機械学会論文集 B 編,Vol.764, No.765, (),.77-777. (9) Nakaas, K., Hg-Densty Mes Flow Comtatons wt Pre-/Post-Data Comressons, AIAA aer 5-4876, (5). () A Program Pakage for ermoysal Proertes of Flds, Ver.., PROPAH GROP. () Wess, J.M. and Smt, W.A., Preondtonng Aled to Varable and Constant Densty Flows, AIAA Jornal, Vol.33, (995),. 5-56. () seng, Y.-H. and Ferzger, J.H., A gost-ell mmersed bondary metod for flow n omlex geometry, Jornal of Comtatonal Pyss,Vol.9, (3),.593 63. (3) Ye,., Mttal, R., daykmar, H. S. and Syy, W., An Arate Cartesan Grd Metod for Vsos Inomressble Flows wt Comlex Immersed Bondares, Jornal of Comtatonal Pyss, Vol.56, (999),.9 4. 46