はじめに 現在，多くの工業製品の設計，開発プロセスに おいて，数値シミュレーションの重要性が非常に 高くなっていることは誰もが認めるところであろ う．特に，流れと熱との連成現象の解析は，熱交 換器をはじめとする熱機器はもちろん，近年ます ます小型化され，発熱密度が増大しているノート 型PC，サーバーを含む電子機器内部などの設計， 開発においても極めて重要である． よく知られているように，偏微分方程式の数値 計算手法は，計算に用いる格子のタイプにより大 きく 2 つに分類することができる．構造格子を用 いる代表的なものは有限差分法である．規則性を 持たない非構造格子は，有限要素法や有限体積法 で用いられる．ただし，有限体積法では構造格子 を用いた計算例も多い． 異なる計算手法について，本当の意味での長所， 短所を議論するのは難しいが，計算量，使用メモ リ量，さらに高精度化が容易であるという点にお いては構造格子を用いる方がよい．一計算点あた りの計算時間を考えるなら，有限差分法が最もパ フォーマンスが高いと言ってよいだろう．一方， 任意形状への適合性を考えた場合には非構造格子 が有利である．しかし，格子が規則性を持つ必要 がないという大きな利点を持つ反面，構造格子の 利点として挙げた事柄が，そのまま非構造格子を 用いる際の不利な点となることに注意が必要であ る． もともと工業製品の設計，開発プロセスにおけ る解析形状は複雑であることが多いが，小型化・ コンパクト化の要求により，ますます複雑化して いるのが現状である．このような複雑な解析対象 に対しては，CADデータからメッシュジェネレ ータを用いて非構造格子を生成し，非構造格子を 使って数値解析を行うというアプローチが一般的 である．しかし，この過程ではデータ形式の変換 やCADデータが持つ不整合に対して修正を行う 必要がある．また，解析対象の形状によっては， 部位ごとに解析形状に見合った要素サイズを設定 しないと格子生成がうまくいかないこともあり， 労力と時間がかかる作業であることは間違いな い．さらに，格子が生成されたとしても，作成者 の熟練度，あるいは解析する現象に対する理解度 によって，その格子の品質が異なるということも 起こる． このような，手作業の介入，格子生成を行う者 ごとのばらつきを排除し，格子生成段階に要する 労力と時間を削減するアプローチとして，最近，
工学的な設計のための流れと熱の数値シミュレーション
Introduction of Computational Simulation Methods of Flow and
Heat Transfer for Engineering Design
＊理化学研究所
白 崎 実
Minoru SHIRAZAKI
＊理化学研究所岩 田 正 子
Masako IWATA
＊理化学研究所姫 野 龍太郎
Ryutaro HIMENO
ながれ 22（2003）247−259. ＊〒351-0198 和光市広沢 2-1 †E-mail : shirazak@postman.riken.go.jp 〔特集〕熱・流体・構造連成問題 †工学的な設計のための流れと熱の数値シミュレーション 直交格子法がクローズアップされている．直交格 子法にはいくつか近似レベルがあるが，基本的に は直交格子における多数の直方体（ボリュームセ ル，Voxel）を用いて解析領域の形状（つまり流 体部分や固体部分）を表現する方法である．完全 自動化が可能であり，計算に必要なメモリ量も少 なくて済むという長所がある反面，境界層の挙動 や物体表面における熱伝達の取り扱いなどは苦手 である． これら 2 つのアプローチは適材適所で，それぞ れ適する対象に使い分けるものである．本稿では， 流れと熱の連成解析に関して，まず第 1 節では並 列計算を行うことを前提に，非構造格子を用いて 流れと熱を統一的に扱う計算手法について述べ， 第2節においては，計算結果に応じて格子サイズ の分布を変化させながら解析を進める解適合型直 交格子による計算手法について紹介する． 1．並列有限要素法による流れと熱の連成解析 流体中に置かれた物体と流体との間の熱伝達に 関する研究は，従来から数多く行われてきており， 実験はもちろん数値計算による解析も多く報告さ れている1）．しかし，物体には表面温度一定，あ るいは表面熱流束一定などの条件が課されたり2）， 実験によって得られた熱伝達係数を与えて計算が 行われたりすること3, 4）が多い．「流体−固体」間 あるいは「流体−固体−流体」間の連成を考える 場合には，このような境界条件は適切とは言えな い．熱的挙動を正確に見積もり，適切な熱設計を 行うためには，熱流体の運動，固体内部の熱伝導， 固体と流体の間の熱伝達を正しく計算する必要が ある． 固体内部の熱伝導と流体との熱伝達を扱った例 としては次のようなものがある．水俣ら5）は，有 限要素法にもとづいた定式化において，流体およ び固体を含めた解析領域を考え，その界面と要素 境界を一致させるようにとれば，流体部分と固体 部分とで異なる物性値を指定することが可能であ ることを示した．その上で，固体部分を有する2 次元キャビティ−流れを解き，実験と比較して妥 当な解を得ている．増田ら6）は，3 次元非構造格 子を用い，有限体積法により熱流体と固体の熱交 換を解いている．ひとつの計算セルの中に流体と 固体が存在する場合も計算できると述べている． 大西ら7, 8）は，単一あるいは複数の円管からなる プレートフィン付き管熱交換器ユニットついて， 管外の流体と管壁との間の熱のやり取りを考慮 し，有限差分法を用いて流動，伝熱特性について 調べている．また，奥田ら9）は，六面体要素を用 いた有限要素法により，流体と固体との熱伝達を 統一的に扱い，データ並列化を行っている． 複雑な解析形状への対応を考えると，非構造格 子を用いた解析が有効である．しかし，すでに述 べたように，非構造格子を用いた解析は，計算量 やメモリ使用量の点で不利であるだけでなく，一 般には，メモリアクセスが連続的にならないため さらに不利となる．解析対象の精密なモデル化や 解析精度の向上といった要求に応えるために，解 析規模はますます大きくなる傾向にある．これら のことから，非構造格子を用いた解析を行う場合 には，並列計算は必須であるといえる．並列計算 を行う場合には，通常，使用するプロセッサ数に 見合った計算速度の向上（並列化効率）が求めら れる．ここでは，並列計算を念頭においた有限要 素法による流れと熱の統一解法について紹介す る．簡単な検証例を示した上で，CADデータか らの解析例を紹介し，最後に「流体−固体−流体」 間での熱のやり取りの例として，二つの円管内外 の流れと熱伝達をとりあげる． 1.1 計算手法 熱流体の流れと固体の熱とを同時に解析するた めの支配方程式として，次に示す非圧縮性流体に 対する連続の式，Navier-Stokes方程式，エネル ギー方程式を用いる．なお，浮力項は無視している． v v v v+ ⋅∇ =−∇ + ∇ ∂ ∂ Re p t 1 ) ( 2 0 = ⋅ ∇ v （1） （2）
白崎 実・岩田正子・姫野龍太郎 ここで，v，p，Tはそれぞれ無次元化された流速 ベクトル，圧力，温度であり，Re，Prは， Rey-nolds数，Prandtl数を表す． 流体と固体とを分離せずにひとつの方程式系と して扱うために，固体部分に常に流速 0 の境界条 件を課し，（1），（2）式の計算の際には流体部分の 寄与のみを考慮する5）．ただし，図 1（A）に示す ように，流体−固体の界面と要素境界面を一致さ せて，この界面をまたぐ要素ができないように格 子を生成しておく．このような非構造格子では， 流体部分の要素と固体部分の要素に個別の物性値 を与えることが可能になる．ここでは（3）式の RePrの値を調整することにより，流体と固体の 熱拡散率の違いを表現することができる．（3）式 では熱に関する他の効果は省略しているが，この 効果を考慮する場合も，流体，固体のそれぞれの 要素において，0 あるいは 1 の値を持つマスキン グパラメータを与えておけば，流体部分あるいは 固体部分でのみ定義される項を導入することがで きる． これらの支配方程式をFractional Step法10）に より時間方向に分解し，有限要素法にもとづいて 空間の離散化を行う．この際，複雑形状への適応 性を考え四面体一次要素を用いる．解析領域はグ ラフ分割ツールMETIS11）により要素をベースと して領域分割され，それぞれの部分領域の計算は 各プロセッサが必要に応じて通信を行いながら並 列に処理される．このメッセージパッシングのた め の 通 信 ラ イ ブ ラ リ と し て はM P I（M e s s a g e Passing Interface） を 使 用 す る ． な お ， 圧 力 Poisson方程式の求解には，共役勾配法（CG法） を用いる． 流体と固体で別々の方程式を立て，境界条件を 通じて両者を連成させる手法は弱連成法と呼ばれ る．この場合には，図 1（B）のような格子が用い られる．これに対して，連成系全体の方程式を一 括して解く方法は強連成法と呼ばれる．一般に弱 連成法では開発済みの既存のコードを利用できる のに対して，強連成法では，新たにコードを開発 する必要がある．このため，強連成法に比べて， 弱連成法の方が連成解析コード開発は容易である といわれる12）．しかし，本手法の場合は，熱流体 解析のコードをほぼそのまま利用して連成解析を 行うことが可能である． 並列環境における弱連成法と本手法との違いに ついて述べておく．弱連成法において複数のプロ セッサを用いた並列化のバリエーションとしては 以下のようなものが考えられる． （A）与えられた複数プロセッサを二つのグルー プに分けて，一方のプロセッサのグループ が「流れ解析」を，もう一方のグループが 「熱解析」を行い，境界条件を交換しながら 解析を進める（図 2）． （B）与えられたすべてのプロセッサを用いて「流 れ解析」→「熱解析」→「流れ解析」→ … というように，境界条件を交換しながら，そ れぞれの解析を交互に行う（図 3 ）． 上記の（A）では，一方の解析を行っている間 は，他方が「待ち状態」になってしまう．この待 ち状態では，必ず一方のグループのプロセッサは 遊んでいることになる．このため，流れ解析，熱 T 1 T ) ( T+ ⋅∇ = ∇ ∂ ∂ RePr t v 2 （3） 図 1 計算に用いる格子
工学的な設計のための流れと熱の数値シミュレーション 解析のそれぞれで仮に 100％ の並列化効率が得ら れたとしても，使用しているプロセッサ全体では， 同様の効率は期待できない．また，（B）の場合， 図 3 からわかるように，プロセッサが待ち状態に なっている時間は発生しない．しかし，一般には 流れ解析と熱解析における計算量（粒度）には差 があるため，与えられたプロセッサ数で，いずれ か一方の解析においては高い並列化効率が得られ たとしても，もう一方の解析では低い並列化効率 しか得られないということが考えられる．すなわ ち，解析全体に対する並列化効率は，計算量が小 さい方の解析に足を引っ張られる可能性がある． 本手法では，形式的に固体部分でも流体解析と同 じ計算を行うことになり，その意味では無駄な計 算を行っていることになる．しかしながら，流体 と固体とを分離せずに同一の支配方程式を扱うこ とになるため，図 4 のように，与えられたプロセ ッサ全てを常に使用することが可能となり，高い 並列化効率が期待できる．本手法では，20 万要素 程度の規模の場合にPentium4を 8 プロセッサ用 いた場合でおよそ 7 倍程度の速度向上比を達成し ている13）．さらに，「流体−固体」間の一対一の やりとりではなく，「複数の流体−複数の固体」 間のやりとりの場合でも，境界条件の交換につい て考慮する必要がない． 1.2 数値計算例 いくつかの数値計算例を紹介する．ここでは， 2 次元計算を除き，Gigabit Ethernetにより接続さ れたPentium4クラスタ16ノード（2.2GHz, RD-RAM 2 GB, Score 4.2.1）を用いて並列計算を行 った． まず，本手法の精度について検証する．温度 T=0 の一様流中にT=1 に加熱された 2 次元円 図 2 弱連成の場合の並列化（A） 図 3 弱連成の場合の並列化（B） 図 4 本手法の場合の並列化 表 1 計算に用いたパラメータ （ 2 次元加熱円柱まわりの流れ） Reynolds数 120 Prandtl数 0.71 時間刻み 0.001 節点数 50,500 要素数 100,000
白崎 実・岩田正子・姫野龍太郎 柱が置かれている．ただし，円柱はその内部まで 格子を生成し，円柱表面とその内部でT=1 であ るとする．表 1 に示した条件で計算を行い，時刻 t=100∼200 における局所Nusselt数の時間平均 を求めた．図 5 は，Eckertらの実験結果14）と比 較したものであるが，よく一致していることがわ かる． 次にCADデータを用いた解析例を示す．キャ ンドルスティックのCADデータからメッシュジ ェネレータを使って図 6 のような格子を生成し た．キャンドルスティックの底面をT=1 に加熱 し，T=0 の一様流中に置いた．そして，表 2 に 示 す 条 件 の も と で 計 算 を 行 っ た ． 図 7 は 時 刻 t =8.75 におけるキャンドルスティック表面の温 度分布および流れと平行なy=0 断面における流 速ベクトルである．ここで，流速ベクトルの濃淡 は流速の大きさを表している．当然のことながら， 物体の下側が高い温度となっているが，底面から の高さが同じであっても，流れによる冷却のため， 下流側の方が高温になっている． 最後に，固体を介して二流体が熱交換を行う場 合の計算例として，図 8 に示すようなモデルにつ いて解析15）を行う．二つの円管の内部を高温の 流体（T=1）が流れており，これらの円管と直 交して，外側を低温の流体（T=0）が流れてい る．二つの円管は，外部の流れに対して直列に L/D=2（D：円管外径, L：中心間距離）で配置 されている．管の外径は 1.0 ，管壁の厚さは 0. 1 である．Reynolds数は，円管の内部流れ，外部 流れともに 120 とした．さらに，本解析では計算 領域全体で，Prandtl数を 0.71 に設定した．計算 条件を表 3 に示す．1,207,136 要素（208,740 節点） 0 2 4 6 8 10 12 -180 -120 -60 0 60 120 180 Eckertet al. Present N us selt Num b θ (degree) θ er 図 5 加熱円柱表面の時間平均Nusselt数 （t=100∼200 の平均）． 表 2 計算に用いたパラメータ （キャンドルスティックまわりの流れ） Reynolds数 100 Prandtl数（流体） 0.71 Prandtl数（固体） 0.14 時間刻み 0.0005 節点数 148,676 要素数 866,574 図 6 キャンドルスティックの有限要素メッシュ． 図 7 キ ャ ン ド ル ス テ ィ ッ ク 表 面 の 温 度 分 布 と y=0断面における流速ベクトル（t=8.75）．
工学的な設計のための流れと熱の数値シミュレーション からなる非構造格子を用いて計算を行った．図 9 に時刻t=120 において，上流側下方から見た円 管の内外壁の温度分布と z=0.5 面の等温度線を示 す．高温流体の流入により，流入口に近いほど， また，管の中心に近いほど，管壁や管内は高温に なっている． 図 10 は，やはり時刻t=120 におけるz= −1.0 面とy=0 断面での温度分布である．二つの円柱 が 直 列 に 配 置 さ れ て い る 場 合 ， 円 柱 間 の 距 離 L/D=2 では，円柱間にはKarman渦は発生しな い16）が，この図から，本計算でも同様の結果と なっていることがわかる．図 11 にそれぞれの円 管外周の最も上流側の位置から測った角度θに対 する時間平均温度分布（t=100∼120）を示す． 上流側の円管と下流側との円管での温度分布の違 いがよくわかる．高温流体の流入口からの距離が 大きくなるにともなって，全体的な温度が低下す ることや，円管の後方（θ = ±180°付近）で温度 が高くなるといった点は上流および下流の円管で 共 通 し て い る ． し か し ， 上 流 側 の 円 管 で は ， θ =0°で最も低い温度であるのに対して，下流側 の管では，θ =0°において極大値を持っているこ とがわかる． 図 12 は，単一円管の場合について，固体部分 のPrandtl数が Pr= 0.71（流体と同じPrandtl数） の場合と，Pr= 0.17 の場合とを比較したものであ る．円管中心を通り流れに対して平行なy =0 断 面における時間平均温度分布（t= 100∼120）をプ ロットしている．高温部分が下流側に偏った分布 となっている．また，流体部分と固体部分での Prandtl数の違い（すなわち熱拡散率の違い）に 表 3 計算に用いたパラメータ （二つの円管内外の流れ） Reynolds数 120 Prandtl数（流体） 0.71 Prandtl数（固体） 0.71 円管の外径 1 管厚 0.1 2円管の中心間距離 2 円管の高さ 2 時間刻み 0.0001 節点数 208,740 要素数 1,207,136 x y z u=1, v=w=0, T=0 u=v=0 , w=1, T=1 1.0 w=0 w=0 v=0 v=0 p=0 p=0 1.0 1.0 7.5 19.5 y=-7.5 0.8 x=-8.0 x=22.0 z=-1.0 z=1.0 y=7.5 図 8 境界条件（二つの円管内外の流れ）． 図 9 管 壁 の 温 度 分 布 と 等 温 度 線 （z=0 . 5面 ， t=120）． 図 10 管壁の温度分布と等温度線（y=0面と z= −1.0面，t=120）．
白崎 実・岩田正子・姫野龍太郎 より，流体と固体の境界において等温度線の屈折 が見られる．また，当然のことながら固体の熱拡 散率が大きい方が，流体によく熱が伝わることも わかる． 2．直交格子法による流体と熱のシミュレーション 直交格子を使う最大のメリットは格子生成にか かる手間が少ないことであるが，物体境界に沿う メッシュを使う場合に比べ計算精度の点では劣る ため，その改良方法について数多くの研究が報告 されている．改良は，大きく二つの方向でなされ ていて，一つは，一般に格子面と一致しない物体 境界を取り扱う方法の研究，もう一つは，格子密 度を部分的に（経済的に）上げる方法の研究であ る． 物体境界の問題については，物体表面が粘性流 体に及ぼす影響を支配方程式の中に取り込むな ど，間接的に表面を表現するPeskin17），
Gold-stein18）によるImmersed Boundary法，Saikiに
よる仮想境界法19），矢部のCIP法20）による方法，
HirtのVOF法21）などがある．あるいはより直接
的に境界形状を表現するいわゆるカットセル法で Quirk22），Berger23），LeVeque24），Forrer25），中
野ら26）によるそれぞれの改良，物体表面と格子 との交点情報を直接スキームに取り込む市川らの 方法27），物体表面の圧力境界条件を正しく評価す る松永らの方法28），また，物体表面のみ境界適合 格子とする朴らの方法29），さらに多方向風上スキ ームにより低コストで形状近似効果を出す桑原の 提案30）などがある． 一方，格子密度を上げることにより精度向上を 目指す方向では，精度向上が必要な領域の格子密 度を部分的に上げる方法に焦点が当てられてい る．小野ら31）は，自動車設計評価に関わる流体 シミュレーションに，必要解像度に応じてブロッ ク 分 け し た マ ル チ レ ベ ル 直 交 格 子 を 用 い た ． Berger32）のAMR法は，構造格子上で，高精度 を要求する領域に自動的に細かい格子のパッチを 当て，小さい時間刻みで計算し，粗い格子と積分 値で受け渡しをする方法であり，内山らが必要精 度を物理量勾配でモニターすることにより衝撃波 の解析を行っている33）．中橋は非構造格子におい て形状・解適合で格子を自動細分化（Refining）， 物体表面層のみ境界適合格子と組み合わせて翼周 りの圧縮性流体の解析を行った34）．Zeeuwはツ リー構造によりセルを管理するツリー型直交格子 で形状・解適合自動格子Refiningを行った35）． ツリー型直交格子においては，流体計算に必要な 隣接セル探索に時間がかかる場合があるが，小川 は，系譜を利用することにより探索時間を短縮し た36）．またビル風の解析等，高解像要請方向が非 図 11 管壁外側面の時間平均温度分布 （t=100∼120の平均）． Te mpera ture x
Tube Wall Tube Wall
0 0.2 0.4 0.6 0.8 1 -1 -0.5 0 0.5 1 Pr=0.71, z=-1.0 Pr=0.71, z=-0.0 Pr=0.71, z=1.0 Pr=0.14, z=-1.0 Pr=0.14, z=0.0 Pr=0.14, z=1.0 Upstream Downstream 図 12 y=0断面における時間平均温度分布 （t=100∼120の平均）．
工学的な設計のための流れと熱の数値シミュレーション 等方な場合に対応して非等方分割を採用した． 理化学研究所では，「ものつくり情報技術統合 化研究プログラム」で直交格子（Voxel）を直接 的に用いて物体を表現し，解析するCAEシステ ムの開発を行っている．この中で表現された直交 格子に基づくデータをV-CADデータと呼び，熱 流体解析のために，有限体積法とVOF法および カットセル法を組み合わせたスキームを開発した 37）．格子密度の向上については，ツリー型直交格 子を用いて，形状適合，解適合でRefining，さ らにCoarseningも取り入れて，熱流体計算を行 う研究を進めている38）． 本節では，その適合ツリー型直交格子における シミュレーションの 2 次元での試みを紹介する． メッシュは，物体表面に近い領域や，物理量勾配 の大きい領域を細かくし，高精度化や安定化に寄 与する．そうでないところは，粗くして計算コス トを下げる．これにより，複雑な形状に対応し， かつ必要な精度を満足する実用的な設計支援シミ ュレーションが可能になると期待している． 2.1 4 分木格子（ 2 次元ツリー型直交格子） まず，メッシュを自動的に生成・変化（リメッ シュ）させていく方針を述べる．V-CADデータ から直交格子と物体表面との交点情報を受け取 り，初めに粗い直交構造格子上に物体を載せ，そ の上で形状適合によりメッシュを再分割する．ツ リー型直交格子では，セル単位で各稜を 2 分割す ることにより新しいセルを生成する．1 回の分割 で 1 つのセルが 2n個のセルに分割される（n=次 元数）ので 2 次元では 4 分木格子となる． 本稿では，オリジナルの直交構造格子セルを 4 分木格子の最も粗なセルとし，ルートセルと呼び， その分割レベルを 0 とする．分割により新しく生 成されるセルの分割レベルは 1 ずつ増加する．細 密レベルは予め設定しておく．流体計算に関与す るのは，末端の分割されていないセルのみであり， これを計算セルと呼ぶ．その分割レベルは 0 から 細密レベルまで分布している． 形状適合リメッシュのルールは，（1）固体流体 境界を含む計算セルとその両隣の計算セルは細密 レベルまで分割，（2）隣り合う計算セルの分割レ ベルの差は 1 以下（Smoothing），とした．細密レ ベル=3 として，このルールに従い円柱を形状適 合リメッシュした様子を図 13 に示す． 次に，解適合リメッシュのルールは，（1），（2） に加えて，（3）物理量勾配の大きい領域の分割レ ベルを上げる（Refining），（4）物理量勾配の小さ い領域の分割レベルを下げる（Coarsening），と した．物理量としては，温度，速さ，速度成分， 渦度，圧力などを用いることが出来る．勾配のチ ェックを一定間隔（図 14 の例では 0.1 無次元時間） で行い，勾配がしきい値以上であればRefining フラグを，しきい値以下であればCoarseningフ ラグを立て，一定間隔（同 1 無次元時間）でフラ グに従いRefining，Coarsening，Smoothingを行 う．このとき，Refining用のしきい値より Coars-ening用のしきい値を低めに設定することにより， 分割レベルが不安定となるセルの発生を抑制する ことができる． なお，Refiningの対象は，分割レベルが細密レ ベルに達していない計算セル．勾配チェックは全 方向の隣接計算セルに対して行い，しきい値以上 の場合が一度でもあれば，自分とその隣接セルに Refiningフラグを立てる．Coarseningの対象は， 1 回だけ分割されているセルである．自分は計算 セルではないので，サブセルの平均値としての自 分の値と隣接セルとの間で勾配チェックを行う． 勾 配 が 全 回 に わ た り し き い 値 以 下 で あ れ ば ， Coarseningフラグが立つことになる． 図 13 細密レベル=3として円柱回りを形状適合 リメッシュした4分木格子．
白崎 実・岩田正子・姫野龍太郎 図 14 は，温度 1 に設定した直径 1 の円柱に向 かって，左空間から温度 0 ，Re=200 の空気がu =1 ，v=0 で流入するときに，細密レベル=2 の 場合で解適合リメッシュをした様子を，時間を追 って示している．図 14-1 は，Refiningのみ行っ た場合，図 14-2 は，RefiningとCoarseningを行 った場合である．Refiningのしきい値は温度勾配 0.011，Coarseningのしきい値は 0.009 である． 時間発展にしたがい，円柱の右側の流体温度が 上昇し，Karman渦の発生につれて温度場も振動 を始める．図 14-1 では，温度勾配大となった領 域で次々と再分割が行われるので，計算セル数は 単調増加し，細密レベルのセルが占める面積は増 大する．一方，図 14-2 では，温度勾配が小さく なった領域でセルの再統合が行われるため，後流 の振動と共に，細密レベルのセルが占めるメッシ ュの領域も振動している． 解適合リメッシュを採用すると，解析者は精度 よい計算を行うためにどの領域のメッシュをどの 程度細かくすればがよいか，予め決定しなくてよ い．流体力学的な勘や経験を持たない開発者が， 小規模な計算機資源で，信頼できる流体解析が出 来る，という期待に応えられると考えている． 2.2 基礎方程式と離散化スキーム ここで，計算に用いた基礎方程式とその離散化 t= 5 t= 20 t= 80 図 14-1 温度勾配によるRefining. t= 5 t= 20 t= 80
工学的な設計のための流れと熱の数値シミュレーション スキームについて簡単に延べておく．基礎方程式 は，非圧縮性のNavier-Stokes方程式（4）と連続の 式（5）およびエネルギー方程式（6）である． Vfはセル内流体体積占有率，Vs=1−Vf，Re はReynolds数，PrはPrandtl数．また，密度ρ， 比熱c，熱伝導率λは，固体の値であり，流体の それぞれの値を代表値として無次元化している． 各式とも，有限体積法に則り，移流項，粘性項 および圧力勾配項はセル内流体部，温度拡散項は セル内で積分し，発散等をグリーンの定理にした がい表面積分で表現して離散化した．格子面と一 致しない物体境界を表現するためには，VOF法 にもとづき，流体体積占有率，セル開口率を用い た．移流項の計算には，一次精度の風上補間を用 いた． 2.3 計算例 図 15 は，A：図 14-1 に示したRefiningのみに よる解適合格子[Mesh (R)]，B：図 14-2 に示した Refining & Coarseningによる解適合格子[Mesh
(R &C ) ]，C： 全 て 細 密 レ ベ ル の 直 交 構 造 格 子
[Voxel Mesh]により，円柱回りの流れを計算した
結果を，抗力係数 (CD)，揚力係数 (CL)の時間変
化で比較したものである．表 4 には，Re=200，
40 におけるCDを実験値と比較した．CDの値は， Refining，Refining & Coarsening，Voxel Mesh
とも，Re=200，40 において実験値のばらつきの
範囲であり，良好な結果が得られたといえる．
速度場，温度場，圧力場については，Refining，
Refining & Coarsening，Voxel Meshで，定性的 な変化はなく，同じような結果が得られた．図 16 にRefining & Coarseningの場合のRe= 200，
t=100 における温度場，図 17 に計算セル数の時
間 変 化 を 示 す ． な お ， 平 均 計 算 セ ル 数 は ， Refining：3,867 ，Refining & Coarsening：3,615 ， Voxel Mesh：16,848 であった． ∂ ∂ = − + − v vv v t ∇⋅( ) ∇ ∇p 1 Re 2 ∇⋅v =0 (V cV) T ( ) t V T V V V Q V Q f s f f s f f s s + ∂ ∂ = − + + + + ρ λ ρ v RePr ∇T ∇
⋅
2 （4） （5） （6） -0.5 0 0.5 1 1.5 2 0 20 40 60 80 100 120 140 160 180 200 time CD CL C:Voxel Mesh A:Mesh B:Mesh (R) (R&C) 図 15 CDとCLの時間変化（Re= 200） Re = 40 Re = 200 1.393 ± 0.004 1.398 ± 0.008 1.375 ± 0.003 1.31 ∼ 1.43 1.676 ± 0.015 1.684 ± 0.012 1.689 ± 0.002 1.68 ∼ 1.83 A: Mesh(R)* B: Mesh(R &C )* C: Voxel Mesh* experiments***: Values are average of t = 100∼200 for Re = 200 , t = 50∼100 for Re = 40.
**: 参考文献39）
表 4 CDin different mesh systems and experiments
白崎 実・岩田正子・姫野龍太郎 おわりに 流体と熱が絡んだ問題を有限要素法と解適合直 交格子法で解く手法と結果を紹介した．固体から 流体への伝熱は固体の表面から境界層を通して行 われるものであり，表面積が正しく評価されなく ては解の信頼性は乏しい．この意味で固体と流体 の境界をメッシュで捕らえる必要があり，従来か らこのように行われてきた．この方向の代表的な 例として有限要素法で解析する場合をあげたので ある．この方法でやることを推奨する．しかしな がら，実際の工業製品では非常に細かい部品やケ ーブル類など，計算メッシュで表現することが非 常に手間を伴い，メッシュ生成に１週間から１ヶ 月という期間を要することや，設計の途中段階で ある程度大雑把な熱収支を見たいこともある．こ のような場合，威力を発揮するのが直交格子にも とづく解析手法である．本稿では，その中でも最 も手軽で，しかも，ある程度の精度の保証のある 解適合直交格子法で解く手法を紹介した．この方 法はまだ２次元で，今後３次元に拡張する必要は あるが，ここで取り上げた２つの方法を状況に応 じて使い分けることで，流体と熱の絡んだ工学的 な問題は解析可能だと考えている． 引 用 文 献 01）例えば, 日本機械学会：伝熱工学資料・改訂 第４版． 02）例えば, 中島 円，柳岡英樹，太田照和：流 路内突起物まわりの流れと熱伝達に関する三 次元数値解析, 日本機械学会論文集 64-628B (1998) 4116-4121. 03）例えば, 中山 恒，朴 商煕：平行平板流路 の壁面に置かれたモジュールからの強制対流 空気への複合熱伝達（熱コンダクタンスの数 値解析予測）, 日本機械学会論文集 61-587B (1995) 2612-2619. 04）小林 孝，大串哲朗，藤井雅雄：自然換気型 電子機器筐体の熱設計手法の研究（発熱密度 の異なる鉛直平行基板群を有する筐体内の自 然対流熱伝達特性）, 日本機械学会論文集 68-669B (2002) 1553-1560. 05）水俣圭子，前田慶太，水上 昭：流体−固体 熱連成問題の有限要素解析，構造工学におけ る数値解析法シンポジウム論文集 15 (1991) 41-44. 06） 増田 糧，中村佳朗：3 次元非構造格子によ る固体を含む熱流体連成計算, 第12回数値流 体力学シンポジウム(1998) 287-288. 07）大西 元，稲岡恭二，松原幸治，中部主敬， 鈴木健二郎：プレートフィン付き管熱交換器 ユニットの流動・伝熱総合特性（単列管に対 する三次元定常数値解析），日本機械学会論 文集64-618B (1998) 534-541. 08）大西 元，稲岡恭二，松原幸治，鈴木健二 郎：プレートフィン付き管熱交換器ユニット の流動・伝熱総合特性（二列管に対する三次 元定常数値解析），日本機械学会論文集 65-634B (1999) 2077-2084. 09）奥田洋司，大城勝史：熱流体／構造系連成問 題の統一解析およびそのデータ並列化に関す る研究，日本計算工学会論文集, Trans. JSCES 1, Paper No. 19990023 (1999).
10）Donea, J., Giuliani, S. & Laval, H.: Finite Element Solution of the Unsteady Navier-Stokes Equations by a Fractional Step Method, Comput. Meths. Appl. Mech. Engrg.
30 (1982) 53-73. 500 1000 1500 2000 2500 3000 3500 4000 4500 5000 0 20 40 60 80 100 Number of Cells Calculated time A:Mesh (R) B:Mesh (R & C) 図 17 計算セル数の時間変化（Re= 200）
工学的な設計のための流れと熱の数値シミュレーション 11）http://www-users.cs.umn.edu/~karypis/metis/ 12）張 群，久田俊明：流体・構造連成有限要素 解析における連成手法に関する検討，日本機 械学会論文集67-662A (2001) 1555-1562. 13）白崎 実，姫野龍太郎：PCクラスタを用い た流れと熱伝達の統一解析, 情報処理学会シ ンポジウムシリーズ2002-4 (2002) 119-124.
14）Eckert, E.R.G. & Soehngen, E.: Distribution of Heat-Transfer Coefficients around Circular Cylinders in Crossflow at Reynolds Numbers from 20 to 500, Trans. ASME 74 (1952) 343-347.
15）白崎 実，姫野龍太郎：管壁を介した二流体 間熱伝達の並列有限要素解析, 熱工学講演会 講演論文集(2002) 97-98.
16）Zdravkovich, M.M.: The Effects of Interference between Circular Cylinders in Cross Flow, J. Fluids and Structures 1 (1987) 239-261.
17）Peskin, C. S.: Numerical Analysis of Blood Flow in the Heart, J. Comput. Phys. 25 (1977) 220-252.
18）Goldstein, D., Handler, R. & Sirovish, L.: Modeling a Non-Slip Flow Boundary with an External Force Field, J. Comput. Phys. 105 (1993) 354-366.
19）Saiki, E. M. & Biringene, S.: Numerical Simulation of a Cylinder in Uniform Flow: Application of a Virtual Boundary Method, J.
Comput. Phys. 123 (1996) 450-465.
20）矢部 孝，肖 鋒：固体・液体・気体の統一 解法とCIP法(2), 数値流体力学 7 (1999)
103-114.
21）Hirt, C. W. & Nichols, B. D.: Volume of Fluid (VOF) Method for the Dynamics of Free Boundaries, J. Comput. Phys. 39 (1981) 201-225.
22）Quirk, J. J.: An Alternative to Unstructured Grids for Computing Gas Dynamic Flows around Arbitrarily Complex Two-Dimensional
Bodies, Computer Fluids 23 (1994) 125-142. 23）Berger, M. J. & LeVeque, R. J.: Stable
Boundary Conditions for Cartesian Grid Calculations, J. Comput. Phys. 120 (1995) 278-304.
24）LeVeque, R. J. et al.: Two-Dimensional Front Tracking Based on High Resolution Wave Propagation Methods, J. Comput. Phys. 123 (1996) 354-368.
25）Forrer, H. & Jeltsch, R.: A Higher-Order Boundary Treatment for Cartesian-Grid Methods, J. Comput. Phys. 140 (1998) 259-277. 26）中野 明，下村信雄，里深信行：デカルト格 子系による任意形状周りの圧縮性粘性流計算, 日本機械学会論文集 61-592B (1995) 4319-4326. 27）市川 治，藤井孝蔵：直交格子を使用した三 次元の任意形状まわりの流体シミュレーショ ン, 日本機械学会論文集 68-669B (2002) 1329-1336.
28）Matsunaga, N., Liu, H. & Himeno, R.: Numerical Analysis of Two-Dimensional Incompressible Viscous Flow in Orthogonal Coordinates, Information 5 (2002) 319-326. 29） 朴 炳湖，黒田成昭：非圧縮性粘性流れの
直交格子解法, ながれ19 (2000) 37-46.
30）Kuwahara, K. et al.: Simulation of High Reynolds Number Flows Using Multidirec-tional Upwind Scheme, AIAA Paper
2002-0133, 40th AIAA Aerospace Sciences Meeting
& Exhibit (2002).
31）Ono, K., Tomita, N., Fujitani, K. & Himeno, R.: An Application of Voxel Modeling Approach to Prediction of Engine Cooling Flow, Society of Automotive Engineers of
Japan, Spring Convention, No. 984 (1998)
165-168.
32）Berger, M. J. & Colella, P.: Local Adaptive Mesh Refinement for Shock Hydrodynamics,
白崎 実・岩田正子・姫野龍太郎
J. Comput. Phys. 82 (1989) 64-84.
33）内山直樹，井上 督：圧縮性流体における渦
構造の数値解析, 第6回数値流体力学シンポ
ジウム講演論文集(1992) 81-89.
34）Nakahashi, K.: An Automatic Grid Generator for the Unstructured Upwind Method, AIAA
9thComputational Fluid Dynamics Conference
89-1985-CP (1989) 533-541.
35）Zeeuw, D. D. & Powell, K. G.: An Adapti-vely Refined Cartesian Mesh Solver for the Euler Equations, J. Comput. Phys. 104 (1993) 56.
36）Ogawa, T.: An Efficient Numerical Algorithm for the Tree-Data Based Flow Solver,
Compu-tational Fluid Dynamics 2000 (2000) 337-342.
37）雷 康斌，岩田正子，野田茂穂，姫野龍太 郎：カットセル直交等間隔格子による非圧縮 粘性流れのシミュレーション（V-CADデー タを直接利用する任意形状流れ場の数値解 析），日本機械学会論文集（B編），掲載予定， 論文No. 02-1151. 38）岩田正子，雷 康斌，姫野龍太郎：直交 4 分 木格子を用いた熱流体数値計算（物理量勾配 に応じたRefiningとCoarseningによるリメ ッシュ），日本機械学会論文集（A編），「メッ シュフリー法とその周辺技術」特集号（2004 年3月号)に投稿中, 論文No. 02-0253. 39）日本機械学会：「機械工学便覧」基礎編 A5 流体工学，丸善(2000) 99. 40）小 野 謙 二 ： 「 数 値 流 体 力 学 ハ ン ド ブ ッ ク 」 （小林敏雄編），丸善 (2003) 10.4 直交格子形 成法，537-545, 12.1.3ボクセル法，586-591.
