工学解析の中で重要な分野でもあるのに、その専門の参考書類が構造系の書籍に比べて極端に少ないのが、有限要素法による熱伝導解析の書籍である。たいがいは、有限要素法構造解析本の一章を飾るに過ぎないことが多い。数少ない熱伝導単独の書籍を開いても内容が豊富とは言い難い。

そんなことで、今回は熱伝導解析分野に関係される読者の一助となればと思ってワンポイントテーマを選んでみた。

今、過渡的な時刻で温度分布を考える、すなわち非定常熱伝導解析の基本式を書き下せば下の通りである。下式は１次元の基本式であるが、シンボリックな式と捉えればいい。むしろ、ここでは左辺に注目するので、２、３次元の問題も話の展開としては同じである。

ここに、α は熱拡散率（熱伝導率／熱容量）であり、熱伝導問題を物理的に考察する際に重要なファクターとなる。ここの話でも後で少し顔を出す。　さて、有限要素法というのは、言うなれば偏微分方程式の解法の１つにすぎないが、式⑴の右辺の解法がそれに当たり、空間の離散化による一連の手続がそうである。左辺の時間項は通常、差分法がとられることが多い。

ちょうど、振動問題で現れる時間の２階微分方程式が“ニューマークの β 法”という便利な公式的表現で定式化できるように、１階微分方程式もまた下式で表現できる。

この式は、時刻 n での既知の温度 T_n の値から、時刻 n＋1 での未知の温度 T_n＋1 を求めようとする式である。T の上にあるドットは時間微分を意味する。また、Δ_t は時間刻みである。

式(2)は仮に温度 T を変位に、温度の時間勾配を速度とみなせば分かりやすい。すなわち、時刻 n と時刻 n＋1 での速度値を使って時刻 n＋1 での変位を求めるというアナロジーとなっている。

式中のパラメータβの値により、それぞれにネーミングのついた有名な時間積分法となる。

• β＝０ ：オイラーの前進差分法、いわゆる陽解法
• β＝1/2：クランク・ニコルソン法（第34話参照）。陰解法の１つ
• β＝１ ：オイラーの完全陰解法

この３法は特に有名な解法で、１階の時間微分式の積分法を説明するどの書籍にも出ているが、β＝2/3もあることを読者はご存知だろうか。このときは、ガレルキン法と呼ばれている。偏微分方程式の数値解法を説明する教科書に出てくる、あのロシア人の名がここにも登場するが、以前はこの方法も結構使用されていて、初期時刻から解析時間の短い範囲では、クランク・ニコルソン法などに比べて精度がいいという報告もある。

さて、各積分法の中でも、陽解法だけが時間刻みΔt の値によっては計算が発散してしまう、数学的には条件付の安定性ということで、汎用コードでは避けられる傾向がある。ただ、陽解法も、よく言われているアルゴリズムの簡単さからくる、プログラムのシンプル性という開発者側のメリットだけでなく、利用者側のメリットも少なからずあるので、他の陰解法とうまく使い分ければいいと筆者は思う。

一方の陰解法の方は、Δt の値に関係しない無条件安定性が数学的に保証されているので、開発者側にはよく利用されているが、利用に当たっては、やはり注意点がある。いくらΔt を気にしなくてもいいといっても、それはあくまでも数学的に安定しているということであって、物理的には問題が残る場合がある。最終的には収束するとしても、時刻歴過程では温度の値が振動したり、物理的にはあり得ない値を取ったりすることがある。これについては、後でまた言及するとして、先に、陰解法の中でもライバル関係の位置にあるクランク・ニコルソン法と完全陰解法について、ちょっと話そう。

この二種の優劣をつけるのは難しい。好みと経験上からくる知識で使い分けるしかないのではないだろうか。筆者自身の経験では、どちらに軍杯を上げるかの判断基準は持てなかった。式(2)を見れば、２時刻間を線形補間するクランク・ニコルソン法の方が完全陰解法よりも直感的によさそうに見える。また、時間刻みが小さい場合は、クランク・ニコルソン法がいいと言及している書籍もあったりする。

しかし、完全陰解法の方がいいという意見も一方にあり、そのもっともらしい理由が次の通りである。β＝１の完全陰解法を解釈すれば、時刻 n では温度が T_n であり、それ以外の時刻 n＋1 までは T_n＋1 となることを意味する(図39-1参照)。すなわち、時刻 n＋0 で急激に温度変化しそれ以後は一定値を維持しているわけである。このカーブが指数関数の描くカーブに似ている。温度場を解析的に求めた場合の応答式には時刻を変数とする指数関数が出てくる。この類似が完全陰解法をよしとする所以である。

次に、陰解法でも時間刻みとメッシュ分割とのバランスが必要という大事な話をしよう。ときおり、温度の時刻歴応答が計算的には安定的に得られても、結果が物理的にあり得ない場面に遭遇することがある。これを具体例で紹介しよう。

いま、図39-2のような疎メッシュの解析モデルを考える。この計算モデルは一定の熱流束を入射する左端側面を除いて周りは断熱されている。Δt＝0.05とし、その他の諸元は図に示してある。

ここで、左端から右へ１メッシュ進んだ位置（L＋）に注目してみよう。外部から熱を与えて１秒後のこの位置の温度は表39-1に示すように何とマイナス温度となってしまう。すなわち、外から熱を与えているにもかかわらず、温度が下がるという奇妙な現象になっているのだ。もっとも、位置 L＋より右側の節点位置でもマイナス温度となっている場所もあるが、こちらの方はそもそも絶対値が小さくて数値誤差とみなすことができる。

次に、図39-2のメッシュを長手方向に細分割して、元のモデルの４倍のメッシュ数としてみよう。このときの結果も表39-1に記載してあるが、空間離散化をより密にしたこちらの方は、位置 L＋からマイナス温度が解消し、かつ、そもそも値そのものも小さかったことを教えてくれる。

ところで、本問題は物理的には半無限固体の表面から熱流束を与えたときの温度分布を求める問題と同等になり、本質的に式⑴が支配方程式である。ただ、この微分方程式の解析解には、フーリエ級数の由来となった平板における解と違って、誤差関数という高等関数が入ってくる。ここらのことは伝熱工学の教科書にはよく記載されていることで、結果が図表化されていることも多い。これらの図表を見るとき、冒頭でいった熱拡散率というパラメータが重要な役割を演じていることが分かる。

上の例題は熱流束を与えた問題なので、厳密にはその解析解を言うべきだが、この解式はちょっと複雑なので、ここでは類似した別の問題にすりかえて話を進める。熱流束の代わりに固体表面を瞬間的にある温度に設定し、その後も保持させる問題とする。この問題の場合も、固体内の温度分布は熱流束の場合と定性的には同じ分布グラフを描く。

さて、代替の解析モデルからは覚えておいたら役に立つのではと思う次の式がある。

γ は“Penetration Depth”と呼ばれている。和名では温度貫入深さまたは温度浸透深さと言えばいいのか。そして、γ が意味するところは、時刻 tにおける固体内部の温度が表面温度の１％となる位置である（図39-3）。

この γ の値を知ることが、メッシュ分割の密度決定で役に立つことにもなる。上の例題での物性値からは、熱拡散率 α＝0.02であり、時刻 t＝1.0ともども γ の式に代入すると、γ＝0.5となる。一方、７分割モデルの１メッシュのサイズは約0.5である。ということは、時刻 t＝1.0における温度分布が占める範囲を１メッシュでモデリングしたことを意味している（上で断ったように、厳密には熱流束入射と温度指定の違いはあり、温度貫入深さの数値に違いはある）。これが、表中にある奇妙な数値となった原因である。あまりに疎メッシュの解析だったのである。

γ は熱拡散率 α の平方根に比例するので、もちろん、α が大きくなる、すなわち、熱容量が小さい場合は大きくなっていく。上の例では比較的、熱容量が大きかったサンプルなので、今度は、この数値を下げてみよう。すなわち、熱拡散率 α＝１としてみよう。このとき、γ＝3.6となって温度貫入深さはモデル全長に近くなり、その結果は表39-2で示すとおりとなる。この場合は、メッシュの粗密が結果に大きな影響を与えることはない。

もちろん、ここで紹介している話は１次元解析なので、FEM ユーザが実際に対象とする２、３次元解析には直接、役に立つ話ではないが、γ のことは、メッシュ分割の際の参考にはなるはずである。

最後に。有限要素法（FEM）が大衆化してしまった現在、安易な気持ちで解析ソフトを利用する FEM ユーザが見受けられる傾向にある。力学的あるいは物理的考察などの予備知識のための勉強をしなくとも、市販のパッケージソフトを買ってくれば、適当にメッシュ分割さえすれば答えを出してくれる、といったお手軽過ぎる考えに陥っている人も見受けられる。

たしかに、それで事無きに終わるケースのあることも事実だが、一方で、いくつも落とし穴があって、予備知識なしでは危険きわまりないことがあるのも有限要素法である。特にここの話にあったようなタイムスケールと空間スケールが絡むような工学問題では要注意である。

2006年６月 記