第67話　陽解法の旨み

本エッセイ第39話で、熱伝導解析と陽解法の組み合わせについて、ほんの少しだけ触れた。一般に陽解法は、その欠点が強調されがちで、陰解法に比べると使用頻度が少ないと思われるが、うまく利用すれば、なかなか捨て難い手法でもある。特に熱伝導解析分野では意外な展開もあって面白い。

読者の中には、数値解析分野に不案内な方もおられるかもしれないので、改めて“陽解法”なるものを紹介しておくと、やはり同第39話、式(2)でβ＝0の場合を見れば分かる通り、ある時刻の物理量(ここの場合は温度)を知りたければ、その直前時刻の既知物理量の情報さえあればよく、しかも、両者が陽な形式で表現される場合の数値解法を言う。

さて、同第39話、式(1)にある1次元の熱伝導方程式を一般化して、さらに離散化した方程式が下の式である。この式が熱伝導問題を表現する偏微分方程式に有限要素法(FEM)を適用した場合の熱平衡方程式である。

式(1)における1階時間微分項を差分項に置き換える、すなわち、

を適用すると、下式となる。

式(3)を見ると、陽解法の形式になっていることが見て取れる。この式をそのまま解いてもいいのだが、それでは、左辺の係数マトリックス [Ｃ] がどんと控えているので、連立方程式の解式が必要になって厄介なことになってしまう。

そこで、陽解法の数値解析では、一般に一つのテクニックを使用するのが通常である。振動解析の場合の集中質量マトリックスのように、左辺項にある係数マトリックス [Ｃ] を対角項だけに非零数値を残す集中マトリックス化するのである。これを集中熱容量マトリックスと言う。すなわち、要素の熱容量を近似的に節点に集中させるのである。そうすると、式(3)に残る行列演算は右辺の熱伝導マトリックス [Ｋ] の乗算項だけとなり、数値解析対象の方程式が非常にシンプルとなり、プログラミングもぐっと楽になる。これが熱伝導解析に陽解法を採用した場合の姿である。ここで、陽解法の長短を述べておくと次の通りである。

［長所］

アルゴリズムが非常にシンプルで、プログラム開発が容易である。
アウトプットが陰解法に比べても正確である。
非線形問題も線形問題と区別なく統一的に扱える。
メッシュサイズが比較的揃っているモデル、熱拡散率が小さい問題などでは、陰解法に比べて驚くほど計算速度が速い。

［短所］

陰解法のように解が無条件安定性を持たないので、解式過程で発散することが多い。
厳しい時間刻みが要求されるので、場合によっては、莫大な計算時間となってしまう。
差分法と違って、メッシュサイズが不揃いなFEMモデルでは計算の安定条件を示す明確な指標がない。

短所のダメージが長所を上回るものだから、あらゆる層を対象とする市販の汎用型FEMソルバーでは陽解法が避けられている理由だと推測する。この感覚は、ひょっとしたら、実際に陽解法ソルバーを開発したり、利用した人でないと分からないかも知れない。筆者にも経験があるが、使用していて腹が立つぐらい、よく発散してくれるものだ。

だから、よほどの理由がない限り、陽解法は捨てられる運命なのだが、実はよほどの理由もあるのである。長所の3番目に記した非線形問題の存在である。一般に非線形問題のアルゴリズムの複雑さは線形問題の比ではない。通常の問題は連立方程式の解式が必要となるので、陰解法を採用していると、なおさら複雑で、このFEMソルバーを開発する人は、かなりスキルを持った人でもない限りできない相談である。利用する側でも制御パラメータの扱いに注意を要して、慎重に利用すべき手法である。一方、陽解法を採用するとなると、そんなにスキルを必要としない。ここが、厳しい条件付き安定性という短所を覚悟しながらも陽解法に魅惑される所以である。

熱伝導解析分野での非線形問題といえば、第一に熱伝導率、熱容量等の熱物性値が温度依存するという問題がある。それから、輻射問題のように物理現象そのものが非線形現象であり、本質的に非線形問題である熱伝導問題がある。

そして、三番目の非線形問題として、相変化の問題がある。液相から固相へ変化する凝固現象、その逆の融解現象を取り扱う問題である。相変化問題では、相変化温度になると、突如として、潜熱の吸収、放出が始まるという強い非線形問題である。相変化問題を陰解法でプログラム開発するとなると、ちょっと厄介な業を必要とする。

ここに登場するのが、温度―エンタルピーの関係を熱平衡式(式(1))に導入する“エンタルピー法”である。エンタルピーというのは、熱力学の物理に出てくるあの“エンタルピー”である。“エントロピー”の方と間違わないでいただきたい、念のため。

エンタルピーというのは、一種のエネルギー概念であり、絶対量でなく相対的量の変化を調べる物理量である。地面の位置が決まって初めて、位置エネルギーの量が決まるように、ある基準を決めてから、その基準からの相対量を評価していくものである。なんだか小難しい説明をしたかも知れないが、熱伝導解析では、下の定義式さえ認識してもらえればいい。

式(1)の熱容量マトリックスを対角項に集中させ(節点集中化)、かつ式(2)と式(4)を考慮すると、下の式のように、熱平衡方程式中の時間微分項にあった温度がエンタルピーに入れ替わり、ぐっと簡明な式になることが容易に分かるであろう。

式(5)がエンタルピー法と陽解法を組み合わせた熱平衡方程式であり、単純な漸化式の形になっていることが分かる。 {Ｈ} は節点毎のエンタルピーを集めたベクトルである。あと準備するものといえば、温度とエンタルピーの関係を表す“エンタルピーカーブ”だけである。

この、(陽解法＋エンタルピー法)の組み合わせ解法は、鋳物の凝固解析で使用されたのが始めのようである¹ 。そこでは、差分法が採用されることが多いようだが、もちろんFEMで利用することに問題はない。

ところで、読者には聞き慣れないと推測されるエンタルピーカーブだが、種明かしすれば、実は単純なものである。図67-1を見ていただきたい。

図67-1 エンタルピーカーブ

この図にある直線は熱容量値が温度に依存しない線形問題でのエンタルピーカーブを表している。直線の勾配＝熱容量値ということは、覚えておいて損はない。また、熱容量値が温度に依存する非線形問題では、エンタルピーカーブはまさにカーブとなることは直ぐ理解できるであろう。

エンタルピー法を利用した陽解法では、初期温度から開始して、エンタルピーカーブのグラフ上を温度→エンタルピー→温度のサイクルを繰り返し求めながら、式(5)を進んでいけばいいのである。なお、先にも言ったように、エンタルピーは相対量なので、利用者は初期温度でのエンタルピー量を任意に定めておけばいい。

さて、相変化が存在する場合の解析であるが、この場合、固相、液相で熱容量値が違ってくるので、相変化点を境にして、エンタルピーカーブの勾配が変化するのは当然として、もう一つ、重要な物理量を導入する必要がある。それは、潜熱のことである。潜熱の存在により、相変化点でエンタルピーカーブに段差が出現することになる。筆者は、この面に詳しくはないのだが、純金属の場合と合金の場合では、潜熱とエンタルピーの段差量の関係が違うとのことである。下に両者のエンタルピーカーブを描いていておく。図中の式にあるLとは潜熱のことである。