Excelによる熱流体

日本航空宇宙学会誌, Vol.53, No.621, pp.299-305, 2005
Excelによる熱流体シミュレーションと可視化＊1
Thermo-Fluid Dynamics Simulation and Visualization using Spreadsheet of Excel
富村 寿夫 ＊２ ･ 平澤 茂樹 ＊３ ･ 岩井 裕＊４ ･ 大村 高弘＊５
小林 健一＊６ ･ 羽田 光明＊７ ･ 吉田 英生＊４
Toshio TOMIMURA, Shigeki HIRASAWA, Hiroshi IWAI, Takahiro OOMURA
Kenichi KOBAYASHI, Mitsuaki HANEDA, and Hideo YOSHIDA
Key Words: Temperature Distribution, Velocity Distribution, Numerical Analysis, Visualization, Excel
科学技術計算では Basic，Fortran あるいは C などが
よく使われている．そして，そこでは当然のことなが
ら，解析のための設計図とも云えるフローチャートに
基づくプログラミングを行い，デバッグを何度も繰り
返した末に目的とする解を得ることができる．しかし，
苦労して得られたそのような解も単に数値の並びと
してしか表示されないので，多くの場合は別途用意し
たグラフィックソフトの助けを借りて図示し，やっと
対象とする現象の全体像が把握できるというのが一
般のようである．
一方，最近，Excel の表計算機能を利用した熱流体
シミュレーションも行われるようになってきており，
特 に ， 東 京 大学の森下悦生先生がお書きになった
「Excel で学ぶ流体力学」1)では数多くの解析例が示さ
れている．Excel による数値解析ではワークシート上
に物理モデルのイメージ描くことが可能であり，差分
化式が書き込まれた内点セルや境界点セルの色分け
などをすることでプログラミングやデバッギングを
ビジュアルに行うことができる．また計算結果は付属
のグラフウィザードにより直ちに図示して確認する
ことができる．このような特徴から，Excel 解析には，
これまで当然とされてきた面倒なプログラミングを
必要としない，視覚支援型の新しい数値解析手法の一
つとして大きな可能性が秘められているともいえる．
著者らは，これまでに「ニュートン法を用いた 1 つ
のセル内での 4 次方程式の解法」2)，「ワークシート
を時系列に対応させた 2 次元非定常数値解析法」3)，
「アイコン化セルによる数値解析法」4)などを新たに
開発し，教育用そして研究用の簡易解析ツールとして
利用している．また，「解析結果のアニメ化手法」な
らびに数多くの「身近な伝熱問題や最先端分野におけ
る伝熱問題への適用例」も開発した．
本記事では，丸善（株）より出版した「エクセルと
マウスでできる熱流体のシミュレーション」5)に掲載
した内容を要約して紹介する．
２．アイコン化セルによる数値解析法
基礎方程式や境界条件などの差分式があらかじめ
書き込まれ，かつその機能がビジュアルに表示された
セルが用意されていると，数値解析を非常に効率的に
行うことができる．ここではそのようなセルを「アイ
コン化セル」と名付け，そのセルを用いた解析コード
の構築方法について，2 次元直交座標系の定常熱伝導
問題を具体例として示す．
図 1 に物理モデルの具体例を示す．熱伝導率λが
20W/mK の縦 0.15m，横 0.12m の矩形領域が，下方か
ら熱流束 qwb= 1×104W/m2 で一様加熱され，単位時間，
y
0.12m
60°C
0.15m
１．はじめに
100°C
λ = 20W/mK
Q = 5×105W/m3
20°C
＊1
©2005 日本航空宇宙学会
平成 17 年 5 月 24 日原稿受理原稿受理.
＊2
九州大学 先導物質化学研究所
＊3
神戸大学 機械工学科
＊4
京都大学 工学研究科 機械工学専攻
＊5
ニチアス（株） 浜松研究所
＊6
明治大学 機械工学科
＊7
（株）レーザックス レーザテクニカルセンター
x
0
qwb = 1×104W/m2
図1
物理モデルと座標系
単位体積当たり Q=5×105W/m3 で一様発熱していると
する．また，左壁，右壁および上壁は，それぞれ，100℃，
20℃および 60℃に保たれているとする．このモデルに
対する 2 次元直交座標系の定常熱伝導方程式は
き込まれている．以下に，数値解析プロセスを示す．
図 2 のブック上で，「挿入」→「ワークシート」と
進むと， 新たに「Sheet1」が作成される．このシート
上で，「ツール」→「オプション」→「計算方法」と
進み，自動，反復計算をチェックし，例えば，最大反
∂ 2T ∂ 2T Q
+ 2 + =0
(1)
2
復回数を 1000，変化の最大値 を 0.001 とする．以下
λ
∂x
∂y
のプロセスによりプログラムをビジュアルに作成し，
で与えられる．以下では，アイコン化セルを用いた解
図 3 に示す結果を得ることにするが，ここでは簡単の
析方法に焦点を絞って説明するために，x 方向，y 方向
ため，アイコン化セルが配置されているシート「2 次
いずれも等間隔メッシュ∆x，∆y で差分化した場合につ
元・定常・直角座標系・等メッシュ」を(a)，新たに追
いて示す．内点(i，j)に関する式(1)の差分式は
加したシート「Sheet1」を(b)と表記する．
1
⎧
2
2
2
2 Q ⎫ (2)
Ti , j =
⎨(Ti +1, j + Ti −1, j )∆y + (Ti , j +1 + Ti , j −1 )∆x + ∆x ∆y
⎬
λ⎭
「Sheet1」上への物理モデルの輪郭の描画 例えば，
2(∆x 2 + ∆y 2 ) ⎩
「セルの書式設定」から「罫線」の「外枠」を選択し，
で与えられる．一様熱流束で加熱される下壁面では
対象とする物理モデルの輪郭を「Sheet1」上にあらか
4Ti , 1 − Ti , 2 + 2qwb ∆y λ
Ti , 0 =
(3)
じめ描いておくと，以下の操作が楽になる．
3
特異点セルへの式の書き込み 1. (a)のセル D2 をコピ
となる．なお，特異点となる凸状，凹状の各コーナー
ーし(b)の D10 に貼り付ける．2. (a)のセル K2 をコピー
の温度に関しては，例えば，
し(b)の P10 に貼り付ける．3. (a)のセル K15 をコピー
T1 , 0 + T0 , 1
(4)
x = 0 , y = 0 : T0 , 0 =
し(b)の P25 に貼り付ける．4. (a)のセル D15 をコピー
2
し(b)の D25 に貼り付ける．
のように，壁面に沿う隣接点の算術平均値とした．
左壁セルへの境界値の入力
図 2 に，アイコン化セルを内点，境界点，特異点な
5.
(a)のセル B7 をコピーし(b)の D11 に貼り付け D11
どに分類して配置した，ブック内のシート例を示す．
に 100 を入力する．6. セル D11 を左クリックし右下
セル内に差分化した式が書き込まれており，例えば，
コーナーにカーソルを移動する．7. フィルハンドル
セル I9，G15，D15 に，それぞれ，式(2)，(3)，(4)が書
(■)を左クリックし D24 ま
でドラッグする．
上壁セルへの境界値の入力
8. (a)のセル F2 をコピーし
(b)の E10 に貼り付け E10 に
60 を入力する．9. セル G11
を左クリックし右下コーナ
ーにカーソルを移動する．
10. フィルハンドル(■)を
左クリックし O10 までドラ
ッグする．
右壁セルへの境界値の入力
11. (a)のセル M7 をコピー
し(b)の P11 に貼り付け P11
図 2 「アイコン化セル」のシート例
に 20 を入力する．12. セル
P11 を左クリックし右下コ
ーナーにカーソルを移動す
る．13. フィルハンドル(■)
を左クリックし P24 までド
ラッグする．
下壁セルへの式の書き込み
14. (a)のセル G15 をコピー
し(b)の E25 に貼り付ける．
15. セル E25 を左クリック
し記号を確認する．16. セ
ル D2， D3， D4 に記号 qwb，
∆y， λを書き込み，対応す
る値 10000， 0.01， 20 を
E2， E3， E4 に入力する．
図 3 「アイコン化セル」による解析結果の例
全体のリンクが張られる前
のエラー (#NAME?) 修正
17. セル E2 を左クリック→「挿入」→「名前」→「定
義」→名前欄に qwb が表示されているので「OK」とす
る．18. セル E3 を左クリック→「挿入」→「名前」→
「定義」→名前欄に∆y が表示されているので「OK」
とする．19. セル E4 を左クリック→「挿入」→「名前」
→「定義」→名前欄にλが表示されているので「OK」
とする．
下壁セルへの式の書き込み (続き) 20. セル E25 を
左クリックし右下コーナーにカーソルを移動する．21.
フィルハンドル(■)を左クリックし O25 までドラッグ
する．
内点セルへの式の書き込み 22. (a)のセル I9 をコピー
し(b)の E11 に貼り付ける．23. セル E11 を左クリック
し記号を確認する． 24. セル F2， F3 に記号∆x， Q
を書き込み，対応する値 0.01， 500000 を G2， G3
に入力する．
全体のリンクが張られる前のエラー (#NAME?) 修正
25. セル G2 を左クリック→「挿入」→「名前」→「定
義」→名前欄に∆x が表示されているので「OK」とす
る．26. セル G3 を左クリック→「挿入」→「名前」
→「定義」→名前欄に Q が表示されているので「OK」
とする．
内点セルへの式の書き込み (続き) 27. セル E11 を
左クリックし右下コーナーにカーソルを移動する．28.
フィルハンドル(■)を左クリックし E24 までドラッグ
する．29. セル E11 から E24 まで選択し，右下コーナ
ーにカーソルを移動する．30. フィルハンドル(■)を
左クリックし O24 までドラッグする．
この段階で反復計算は全て終了している．その結果
は，グラフウィザードを使って直ちに確認することが
でき，図 3 のグラフがその出力例である．
合の内点(i, j)に関する式(5)の差分式は
Tik, +j 1 =
+
(
⎧ cx k +1
Ti +1, j + Tik−1+,1 j + Tik+1, j + Tik−1,
⎨
1 + cx + c y ⎩ 2
1
cy
2
(T
k +1
k +1
k
k
i , j +1 + Ti , j −1 + Ti , j +1 + Ti , j −1
⎛ ∂ 2T ∂ 2T
∂T
= a⎜⎜ 2 + 2
∂t
∂y
⎝ ∂x
⎞
⎟
⎟
⎠
⎫
− c y Tik, j ⎬
⎭
)
y
y0
T (x, y, t)
Q
c ρ λ
x
0
図4
t (k)
Time
Worksheet
x0
物理モデルと座標系
Crank-Nicolson
Molecular structure of Ti, jk+1
t (k)
T(k+1)
T(k)
(5)
ここでは，Crank-Nicolson の陰解法を適用し，空間微
分に関しては中心差分近似により離散化した．この場
x
)
(6)
で与えられる．ここで，cx=a∆t/∆x2，cy=a∆t/∆y2 である．
複数のワークシートを時系列に対応させた 2 次元非
定常解析プロセスの概要を図 5 に示す．ここでのプロ
セスは，以下のようになる．
[1] まず，図 6 と図 7 に例示するような，計算パラメ
ータに関するワークシート P と初期条件に関するワー
クシート T(0)を用意する．
[2] 次に，編集機能に含まれる「シートの移動または
コピー」によりワークシート T(0)の複製 T(1) (図 8 参
照) を作成する．そして，セル D5 に Th の値(この例で
は 50)を入力したあと D5 を左クリックし，右下コーナ
ーにカーソルを移動，そしてフィルハンドル(■)を左
クリックし N5 までドラッグする．
[3] 式(6)で記述される内点の温度に関し，ワークシー
ト T(1)上の各セルに時刻 t=∆t における差分式を書き込
む．例えば，まずセル E6 (図 8 参照) に，ワークシー
ト P と T(0)上に配置した必要な値を参照し，式(6)を作
成する．その際，cx，cy，Ti-1,jk+1，Ti,jk などは，その値
が入っているセルを左クリックするだけで自動的に
代入される．
[4] 上記[3]の式(6)が書き込まれたセル E6 を選択し
E14 までドラッグする．さらに，E6 から E14 まで選
３．ワークシートを時系列に対応させた 2 次元
非定常数値解析法
1 次元の非定常数値解析は，例えば，時系列を各行
に対応させることにより，1 枚のワークシートを用い
て簡単に実行することができる．しかし，2 次元の場
合となると，1 枚のワークシートによる非定常解析は
原理的に可能ではあるものの実用的であるとは言い
難い．ここでは，複数のワークシートを時系列に対応
させた 2 次元の非定常数値解析を，直角座標系での熱
伝導問題を例として説明する．
図 4 に物理モデルを示す．初期温度 T0，比熱 c，密
度ρ，熱伝導率λ (温度伝導率 a) の矩形領域 (x0×y0)の
一辺が，時刻 t > 0 で温度 Th に加熱されるとする．こ
のモデルに対する 2 次元直角座標系の非定常熱伝導方
程式は次式で与えられる．
)+ (1 − c
j
x (i)
T(2)
T(1)
T(0)
P
図5
y (j)
y (j)
x (i)
Initial condition
Parameters
解析プロセスの概要
図6
図8
計算パラメータ「P」の例
時刻 t=∆t における解析結果「T(1)」の例
択し M14 までドラッグすると，全ての内点に式(6)が
書き込まれる．この段階で時刻 t=∆t における反復計算
は全て終了している．その結果は，グラフウィザード
を使って直ちに確認することができ，図 8 の下部に示
したグラフがその出力例である．
[5] 上記[2]～[4]で作成したワークシート T(1)を，「シ
ートの移動またはコピー」機能により複製し，ワーク
シート T(2)を作成する．そして，内点の各差分式の中
に含まれる T(0)を[編集]メニューの[置換] に含まれ
る[すべて置換]ボタンにより T(1)に変更し，時刻 t=
2∆t における計算を実行する．なお，ワークシートの
コピーは，[Ctrl]キーを押しながらタブを左クリック
し右側に移動することで簡単に行うことができる．
[6] 必要な時刻まで，上記[5]の操作を反復実行する．
以上から，一旦，所定の時刻までのワークシート群
すなわちブックが完成すると，その後は計算したいパ
ラメータの値を変更するだけで自動的に全時間ステ
ップの計算が実行され，計算結果の数値だけをコピー
してコンパクトに保存することもできる．なお，時間
や空間の分割数が多くなる場合には，それぞれに不等
間隔メッシュを適用することにより，分割数の低減を
図ることができる．
４．解析結果のアニメ化手法
Excel による温度等の計算結果のグラフを，動画で
プレゼンテーションするための，アニメーションの作
成方法を紹介する．
図7
初期条件「T(0)」の例
基本的な方法はパラパラ漫画と同じ原理である．各
グラフを 1 枚ずつの画像ファイルにして連続的に表示
する．
非定常温度シミュレーション（後述の「５.２ アー
ク溶接時の温度分布」など）のように，時間ごとの計
算結果が各シートに同じ形式，寸法でグラフ化されて
いるとして，具体的な作成方法について説明する．図
9 にアニメーションの作成手順を示す．図 9 にある４
つ方法を以下に説明する．
（１）PowerPoint の自動プレゼンテーションを利用す
る方法
まず，少し時間がかかってたいへんだが，特別なソ
フトを使用しなくても，PowerPoint の自動プレゼンテ
ーションを利用し，手作業でアニメーションを表示す
る簡単な方法を次に説明する．
①Excel シートからグラフをコピ－・・・，同じサイ
ズのグラフが各シートに１枚ずつある．
シ ー ト を 選 択 ⇒ グラフを選択 ⇒ 編集 ⇒
コピー
②PowerPoint にグラフを貼り付け
挿入 ⇒ 新しいスライド ⇒ 編集 ⇒ 貼り付
け
・・・シート（グラフ）の枚数分①～②を繰り返す．
③PowerPoint の設定
スライドショー ⇒ スライドショーの設定・・・「自
動プレゼンテーション」，「保存済みのタイミング」
に設定する．
スライドショー ⇒ 画面切り替えの設定・・・「自
動的に切り替え」に設定する（表示間隔（秒単位））
④PowerPoint のファイルを保存
ファイル ⇒ 名前を付けて保存
⑤アニメーションを表示
・・・PowerPoint を起動してファイルを開く．
（２）画像ファイル（JPEG）をアニメーションにする
方法
これ以降の方法については概要のみ説明する．詳細
や実例は文献 5）を参照して頂きたい．
各シートから取り出したグラフを PowerPoint のスラ
イドに 1 枚ずつ貼り付け，画像ファイル（JPEG）とし
て保存する．JPEG の画像ファイルからは，HTML 形
式のファイルを作成して，ブラウザでアニメーション
を表示する．
（３）画像ファイル（GIF）をアニメーションにする
方法
各シートから取り出したグラフを PowerPoint のスラ
イドに 1 枚ずつ貼り付け，画像ファイル（GIF）とし
て保存する．GIF の画像ファイルからは，フリーソフ
トでアニメーション GIF 形式のファイルを作成して，
ブラウザ，あるいは PowerPoint でアニメーションを表
示する．
を表示する．
（４）Excel VBA により画像ファイルを取得する方法
各シートからの画像ファイルの取得を手作業では
なく，Excel のマクロを利用する方法もある．
濃度拡散計算
相似則
12
ホットプレートの温度 発熱量のＰＩＤ制御
制御
昇との関係を計算する．
計算結果の例を図 11，12 に示す．基板上に３個の
発熱する電子デバイスを取り付けた場合，発熱する電
子デバイスやその近傍の配線基板の温度が上昇し（図
11），発熱部品を冷却した空気は温度上昇して流れ去
る（図 12）．
手作業による方法
各グラフをＥｘｃｅｌシートから
１枚ずつ取り出す
各グラフを１枚ずつ
ＰｏｗｅｒＰｉｏｎｔに貼り付ける
ＪＰＥＧファイルで保存
５．伝熱問題への適用例
ＨＴＭＬ形式の
ファイルを作成
以下に身近な伝熱問題や最先端分野における伝熱
問題への適用例のいくつかの概要を説明する．各適用
例の Excel シートや他の多数の適用例（表１参照）は
文献 5）に添付の CD-ROM を参照して頂きたい．
５.１ 基板上の発熱部品の空冷
パソコンなどの電子機器内部には半導体部品など
の発熱する電子デバイスが配線基板に多数搭載され，
狭い空間内に納められている（図 10）．電子デバイス
が発熱することによって温度上昇し損傷することを
防止するため，配線基板の表面に沿った方向に空気を
流し冷却する．高発熱の電子デバイスにはフィンをつ
けて放熱面積を増やす．ここでは，冷却空気を流す配
線基板に搭載する電子デバイスの配置条件と温度上
No.
1
2
3
4
5
6
7
8
9
10
11
ＥｘｃｅｌＶＢＡによる方法
各グラフをマクロで自動的に
画像ファイルにして保存
ＰｏｗｅｒＰｏｉｎｔの
自動プレゼンテーションで
アニメーション表示
図9
アニメーションＧＩＦ形式の
ファイルを作成（その１）（その２）
Ｉｎｔｅｒｎｅｔ
Ｅｘｐｌｏｒｅｒまたは
ＰｏｗｅｒＰｏｉｎｔ
に挿入して
アニメーション表示
Ｉｎｔｅｒｎｅｔ
Ｅｘｐｌｏｒｅｒで
アニメーション表示
アニメーションの作成手順
発熱部品
冷却風
表１ 本記事で紹介した以外の適用例
内容
含んでいる伝熱現象
ふとんの伝熱現象
定常熱伝導
ゆで卵の伝熱現象
非定常熱伝導
高温炉内の伝熱現象
放射伝熱
含水した断熱材の伝熱 非定常熱伝導
配管の断熱問題
非定常熱伝導
摩擦攪拌接合時の温度 強制対流伝熱
分布
電子機器筐体内部の冷 強制対流伝熱
却空気流れと温度分布
流れ関数を使った自然 自然対流伝熱
対流と温度計算
行列を用いたふく射計 ふく射伝熱
算
ニュートン法を用いた ふく射伝熱
ふく射計算
海岸の汚染物質流入の 熱 伝 導 と 濃 度 拡 散 の
ＧＩＦファイルで保存
配線基板
図10
基板上の発熱部品空冷の計算モデル
温度上昇
S1
S1
S2
S2
S3
S3
S4
S4
S5
S5
S6
S6
120.0 -150.0
S7
90.0 -120.0
60.0 -90.0
S8
30.0 -60.0
S9
0.0 -30.0
200.0 -250.0
150.0
S7 -200.0
100.0 -150.0
50.0S8-100.0
0.0 -50.0
S9
1
2
3
図11
4
5
6
7
8
9
S10
10
発熱部品，配線基
板温度分布計算結果
1
2
3
図12
4
5
6
7
8
9
空気温度分布
計算結果
S10
10
使用している基板の定常熱伝導（上面にて対流熱伝
達ある場合）の計算式を次に示す
1800
1600
1400
1200
1000
800
600
400
200
0
時刻１
S9
9
S11
S1
S3
S5
S9
S11
S5
9
S7
1
5
S7
1
5
時刻２
溶接トーチがある位置が局所
的に溶ける温度になる
1800
1600
1400
1200
1000
800
600
400
200
0
1
S9
9
S11
S9
9
S7
5
S5
S1
1
5
S3
1800
1600
1400
1200
1000
800
600
400
200
0
時刻３
時刻４
図 14 鉄板の温度分布の時間変化
5.00
図 15
鉄板の温度分布
0.60
1.40
-0.20
-1.00
-2.60
-1.80
-3.40
-4.20
-5.00
-6.60
-5.80
0.00
1500
1350
1200
1050
900
750
600
450
300
150
0
-7.40
ここで，T：温度，Q：発熱密度，ａ：温度伝導率，ｃ：
比熱，ρ：密度である．
次に，上記の非定常熱伝導の式を解く方法とは別の
方法で，アーク溶接時の温度分布を計算する方法を説
明する．その方法は，アークの発熱分布を点発熱の集
まりと考え，移動点発熱源の準定常温度分布における
理論解の総和として温度分布を計算するものである．
半無限板の表面に，一定速度 v で移動する点熱源が
ある場合の準定常温度分布の理論解は次式となる．
T＝(q/2πλr)*exp{-v*(x+r)/2a}
(9)
r＝(x2+y2+z2)0.5
(10)
ここで，T：初期温度からの温度上昇，q：点発熱量，
r：発熱源からの距離，x：発熱源から移動方向の距離，
λ：熱伝導率，ａ：温度伝導率である．この移動点発
熱源の理論解の総和として全体の温度分布を計算す
る．
図 15 に溶接トーチが一定発熱で，移動する場合の
温度分布の計算結果を示す．
アーク溶接
S11
(8)
1800
1600
1400
1200
1000
800
600
400
200
0
-8.20
⎞ Q
⎟+
⎟ cρ
⎠
図 13
-9.80
-9.00
2
溶融液
S5
⎛∂ T ∂ T
∂T
= a⎜ 2 + 2
⎜ ∂x
∂t
∂y
⎝
2
鉄板
S7
ここで，T：基板の温度，Tair：空気の温度，Q：発熱
密度，λ：基板の熱伝導率，α：対流熱伝達率，b：
基板の厚さである．空気と基板の対流熱伝達率は，空
気流のレイノルズ数により層流か乱流かを区別して
計算式を変えるようになっている．
この応用例の特徴は次である．温度計算は Excel の
表計算機能を使うが，さらにマクロ機能を使って計算
条件の設定を行っている．すなわち，発熱部品の位置
や寸法を座標数値データで与え，マクロ機能ボタン
「部品条件設定」を押すことで，対応する基板のセル
に発熱部品の発熱密度，対流熱伝達率が記入される．
このマクロ機能を使うことにより，発熱部品の寸法を
一定としたまま，発熱部品の配置を変える計算条件の
設定が容易にできる．さらに，マクロ機能ボタン「基
板温度セルの式設定」を使うことで，各セルの周囲条
件（中央，左側境界，コーナ，切り欠きなど）に応じ
てセルに式を入力していくことができる．
５.２ アーク溶接時の温度分布
２枚の鉄板の突き合せ部をアーク溶接するときの
温度分布を計算する（図 13）．アーク溶接とは移動す
る溶接トーチと鉄板との間にアーク放電を起こし，そ
のアークの発熱で鉄板を局所的に溶かして接合させ
るものである．熱は鉄板内を伝導する．
計算結果の例を図 14 に示す．アークの発熱量を時
間変化させて溶接トーチを移動させた場合の温度分
布を計算した．溶接トーチがある位置が局所的に溶け
る温度になり，発熱量に依存して溶融領域の幅が時間
変化する．
非定常熱伝導の式を次に示す．
溶接トーチ
アーク
-10.60
(7)
-11.40
α(Tair - T )
=0
bλ
S1
λ
+
S3
Q
S1
∂y 2
+
S3
∂ 2T
-12.20
∂x 2
+
-13.00
∂ 2T
移動
-5.00
1350.00 -1500.00
1200.00 -1350.00
1050.00 -1200.00
900.00 -1050.00
750.00 -900.00
600.00 -750.00
450.00 -600.00
300.00 -450.00
150.00 -300.00
0.00 -150.00
６．おわりに
4) 富村寿夫，ほか２名, Excel の表計算機能を利用し
た数値計算プロセスの可視化(3), 可視化情報，
24-Suppl.1 (2004), pp.47-50.
5) 岩井裕，ほか６名，エクセルとマウスでできる熱流
体のシミュレーション，丸善（2005）
製氷室
水→氷
図 16
計算モデル
20.0
中心温度[℃]
15.0
水が氷になる時
の凝固潜熱のた
め温度低下速度
が遅くなる
10.0
5.0
0.0
-5.0
-10.0
0
図 17
10
20
時間[min]
30
40
中心位置の温度の時間変化
2.0
2.0
1.8
1.6
1.8
1.6
1.4
1.4
1.2
1.0
0.8
0.6
0.4
0.2
y [nm]
y [nm]
５.３ 冷蔵庫における氷の製成
冷蔵庫の製氷室に入れた水が，冷えて氷になるまで
の温度変化を計算する（図 16）．水は氷になるときに
凝固潜熱を放出する．また，水と氷は物性値が異なる．
計算結果の例を図 17 に示す．製氷室の温度-10℃，
初期の水の温度 20℃について，水（もしくは氷）の中
心温度の時間変化を示す．２０分が経過したころには，
水が氷になる時の凝固潜熱が発生するため，温度低下
速度が遅くなっていることがわかる．
エンタルピーH を用いた非定常熱伝導の計算式を次
に示す．
ρ(∂H/∂t) =λ(∂2T/∂x2) +λ(∂2T/∂y2)
(11)
ここで，T：温度，λ：熱伝導率，ρ：密度である．
氷点 0℃基準のエンタルピーH と温度 T との関係は凝
固潜熱をＬとした場合に次となる．
H＝Ｌ＋c(水)＊T
（T≧0℃）
(12)
H＝c(氷)＊T
（T＜0℃）
(13)
５.４ 気体分子拡散計算
密閉された空間内において，最初に右半分と左半分
にある気体の分子が，時間の経過と共に拡散し混合し
ていく分子の位置の変化を分子動力学法により計算
する．
図 18 は４角形空間内において最初の気体分子の位
置（左の気体分子を薄い色,右を濃い色で示す），図
19 は時間が経過した後の気体分子の位置を示す．拡散
によって分子が移動する様子がわかる．
計算方法を次に示す．分子間には距離に依存した相
互作用がある．分子間ポテンシャルをφとして次式の
Lennard-Jones(12-6)ポテンシャルが良く使われる．
φ=4ε{(σ／r)12－(σ／r)6}
(14)
ここで，ε：分子間の結合エネルギ，σ：分子直径，
r：分子間距離である．すべての分子間ポテンシャルを
加えあわせて，各分子の運動方程式によって分子の位
置変化を計算する．
1.2
1.0
0.8
0.6
分子拡散
0.4
0.2
学生や技術者などのほとんどの人が Excel を使える
0.0
0.0
時代であり，Excel の表計算機能を利用した熱流体シ
0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0
0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0
ミュレーションは，その計算結果をすぐにビジュアル
x [nm]
x [nm]
に表示でき，新しい数値解析手法の一つとして大きな
図 18 初期状態の分子位置
図 19 反復計算後の分子位置
可能性がある．本紹介記事により，読者の皆様が自分
流の効率的な数値解析法を開発される上での「きっか
け」となることを切に願う．
参考文献
1) 森下悦生：Excel で学ぶ流体力学，丸善 (2000).
2) 富村寿夫，ほか２名, Excel の表計算機能を利用し
た数値計算プロセスの可視化(1), 可視化情報，
23-Suppl.3 (2003), CD-ROM 5-3.
3) 富村寿夫，ほか２名, : Excel の表計算機能を利用し
た数値計算プロセスの可視化(2), 可視化情報，
23-Suppl.2 (2003), pp.161-164.