偏微分方程式にトライ(熱伝導方程式1D)

B1  熱伝導方程式(一次元)

長い金属棒(針金)を伝わる熱伝導を考える(拡散現象も同じ系統である)。長さについては、無限長や半無限長は工学的に重要では無い(数学的にはフーリエ変換を用いるので興味がないわけではない)ので除外し、物理的に長さを指定して考えていく。解plotをみると、初期条件として与えた温度分布が物体内に拡散されて次第になだらかになっていく様子が判るはずだ。

解析的には波動方程式と同様に「フーリエの方法」(変数分離、境界条件の元でODEを解く、解の重ね合わせ、初期条件を連動させて解く、フーリエ級数の利用)が利用される。ほとんどの国産のテキストではこの「フーリエの方法」が解説されている。今回参考にした図書のうちペトロフスキーテキストだけ、「このフーリエの方法に欠点がある」として、より一般化された別角度からのアプローチが追加されていた(最終的にはフーリエの方法に持ち込むのだが)。またコシリヤコフテキストではさらに一般化された例題が多く記載されていた。

工学を意識してより一般化を模索していたら、思っていたより熱伝導方程式は奥が深かった。数学的思考はやはりロシア人が優れているんだろうね。今から2Dを纏められるか心配である。

 

解説

今回は話の流れ上、直ぐに結果を確認したい方には申し訳ないが、解説からスタートしその後に計算例を掲載したい。Mapleで記述した下の解説文書では途中の計算を省かず丁寧に記載したつもりである。このブログの性格上、数学的に深く議論はしていない。けれどこれらの方程式の使用制限を知っていても損はないはずだ。

フーリエは、自ら熱伝導方程式を導き、境界条件をゼロと仮定(両端で断熱)したことで後の計算を簡単にし、最終的にフーリエ級数を利用して解けるようにした。これがフーリエの功績だが、工学的には必ずしも両端が断熱とは限らないから、より一般的なアプローチが必要と考える。

------------------------------------------------------------------------------------

-----------------------------------------------------------------------------------

 

計算例

(1)数式処理ツールMapleによる計算結果を以下に示す、今回は解析解。PDEなのでpdsolveコマンド類を使っている。境界条件や初期条件はリスト[ ]として扱えばunionは要らないようだ。私の使っているバージョンではpiecewise関数がpdsolveコマンド内では利用できなかった、もう少しスマートに終わらせたかったが残念。

-------------------------------------------------------------------------------------

------------------------------------------------------------------------------------

 

(2)次はWolfram CloudMathematicaだ。今回の計算ではNDSolveValueコマンドを利用した(従来のNDSolveではPDEが扱えない)。解析的ではなく数値解析結果になる。

condは初期条件と境界条件をまとめたもの。Piecewise関数が使えるのはMapleより親切だ。具体的な数値代入には/.{ L->1など}を使用する。結果plot(メモリーの都合でt=50までの計算)を見るとMapleと同様に初期温度分布が時間経過に従い減衰し拡散している形が見える。数値解析なのでMapleよりスッキリ見える。解析解にはトライしていない。

------------------------------------------------------------------------------------

---------------------------------------------------------------------------------

 

(3)次は有限要素法(以下FEM)ベースのFlexPDEの計算結果だ、つまりFEM数値解析結果である。余り推薦できない方法であるが、y方向を時間tに見立てて構成方程式を入力している。

プログラムコードの補足を加えたい。SELECT文にNGRID=39を指定しているのは、精度向上のため、通常のグリッド数より多くしたいからだ。DEFINITIONSでは構成式中の定数などを指定する。本問題では計算空間の大きさとしてxは1m、y(つまり時間)は50sec、そして熱伝導率、密度、比熱、最後にTgiven=200+273.15degKを指定。単位系は全てSIである。

BOUDARIESには境界条件を指定する。本問題では大きさはx方向、時間軸はy方向の二次元表現である。FEMはこの二次元平面を要素分割して計算を進めていく。start(0,0)は計算スタート地点の座標で、このスタート地点から二次元平面を反時計回りに計算していく。最後のLine to closeで二次元平面を閉じている。座標間に初期条件境界条件を指示する。valueはディリクレ形境界条件。PLOTS文のelevation関数で、境界条件を指定した部分のど真ん中をy(時間)=0〜50secまでplotした。最後はENDでプログラムを閉じる。

プログラムコードを下に、続けて結果plotを掲載する。

結果plotで気になる点がある、y=0付近で本来なら473.15から滑らかに減衰する所がスタート地点でポッコリしている。今の所解決策が見つけられない。

------------------------------------------------------------------------------------

Code

plot

------------------------------------------------------------------------------------

Mapleの解析解、Mathematicaの数値解析解(差分法?)、FlexPDEのFEM数値解析解のいずれも同様なplotを出力してきた。

 

今回、参考にした図書を挙げておこう。

伝熱の物性係数はJSMEテキストシリーズが重宝した、価格も手頃。

 

Mapleによる解説文書は和達氏のテキストが参考になった。

 

具体例では馬場氏のテキストが役に立った。

 

より一般的(境界条件が非ゼロ)な方法は以下の2冊が参考になる。