2009-12-10[n年前へ]
■エクセルでシミュレーション Vol.8 二次元非定常熱伝導問題を解こう 
表計算ソフト（スプレッドシート）プログラムであるMIcrosoft Excelを使い、（VBAはもちろん関数も使わず）ただセル間の加減乗除演算のみで、二次元非定常熱伝導問題を解きで遊んでみることにしてみました。伝熱方程式を差分化し、二次元の非定常解析を行ってみる、というものです。二次元といっても、比較的薄い円形ベルト（シート）上のものが回転するというモデルで、ベルト（シート）の厚みが（ベルト内の熱拡散に対して）十分に薄いという仮定の下での近似を行った、と考えれば、（そういう条件下の）三次元非定常熱伝導問題を解いているとも言えるシミュレーション計算モデルです。
計算領域は右のようなものになり、2本のローラにより回転駆動させている「下から5点のヒータで温められている茶色のベルト」という計算するモデルです。その概略を大雑把に描いてみたのが、上図になります。
というわけで、下の動画が、「ベルトが回転し動き」「途中にあるヒーターでベルトを温め」「ベルトの平面内では熱伝導が生じ」「周りの空気とベルトの境界から熱が放熱され」・・・といった2.2秒の間に生じることをすべて計算し・結果刻々と表示するシートで非定常伝熱計算を行っているようすです。（熱伝導率等はテキトーな値です）
ちなみに、このシートは、毎年行われている「スプレッドシートを用いた電子写真シミュレーション実習」テキストの「有限差分法による定着プロセスの熱伝導計算（非定常問題）」で使う教材用シート（関連テキスト
）にもとづいています。
動画の上の部分が、10×20の領域に差分化されたメッシュの温度分布を数値で示しており、動画の下部部分には温度分布をグラフで示しています。このグラフを眺めていれば、20℃の室温下で5点のヒータにより温度を上げられていく（回転する）ベルトが、次第に温度分布を持ちながらやがて定常的な温度分布になっていくさまを見てとることができます。
この計算シートは、熱伝導方程式
をCFL条件（Courant-Friedrichs-Lewy Condition)下で陽的に差分化し、エクセルの反復計算を用いることで、解いています。しかも、メッシュ間の計算は、ヤコビ法を用いた計算を行っています。
「メッシュ間の計算は、ヤコビ法を用いた計算を行っている」と、ひとことで書きましたが、これは実は凄いことです。なぜかというと、エクセルは（セル間の循環参照がある場合に用いる）反復（収束）計算時には、表の左上をスタート地点として「Zの法則
定常問題を解くならガウスザイデル法で収束するまで反復計算を行うので構いませんが、非定常問題を解くのにガウスザイデル法を用いてしまうと、正確ではない間違った温度傾斜が生じてしまいます。それは、少し困ります（しかし、現在出版されている”表計算＋差分化 熱伝導”書籍は、多くがガウスザイデル法を用いています）。
というわけで、この計算シートではグラフの下に（動画の）上のメッシュ部分をコピーするだけのセルがあり、上のメッシュ部分は「自分の周囲に相当する下のメッシュを用いてメッシュ間の計算を行う」という仕組みになっています。こうするとなぜ、ガウスザイデル法でなく、ヤコビ法の計算になるかは、「エクセルの反復計算時のセル間の計算順序」を念頭に置いて考えてみると、「なるほど！」と納得できると思います。「道具を理解し使いこなす」ということの意味を、心から実感させられます。
だから、この二次元非定常熱伝導問題を解くエクセル・シートは、とても「単純明快な計算がされているだけ」なのですが、実に巧妙に考え抜かれた上で作られた「素晴らしき一品」です。しかも、セルを少しづつ書き変えていくことで、かなり実用的な計算もできるのです。このシートを見たときには、本当に驚かされました。
「どんな風に書きかえると何がわかるか」とか「さらにこんなこともできる」といったことは、明日にでも書いてみようと思います。（というわけで、この続きとして、温度が時間を追って上昇していく機能も付けてみました）
