非定常熱伝導解析のエクセルよるグラフの作成

初期温度0度の縦10mm、横20mmの長方形の左右の辺に、いきなり300度の温度を与えた与えた続けた場合の解析です。ansysではステップ加熱を選択しています。このansysでの解析計算のグラフとエクセルでの計算結果のグラフの比較したいのですがエクセルで使う計算式が分かりません。グラフは縦5mmの場所から横方向に0.5mm幅で温度を取りたいです。この時の計算式を教えて欲しいです。今回は非定常なので1秒ごとに解析を行い、１秒から10秒まで解析を行いました。材料特性は熱伝導率λ＝50(W/m・K)、密度ρ＝7300(kg/m^3)、比熱Cp＝520(J/kg・K)です。

質問者からの補足コメント

• 加熱していない上辺と下辺は断熱です。表計算ソフトでチャレンジしてみたんですが計算式するための式が分からず、とまっています。良かったら表計算ソフトで計算に使った式を教えて頂けたらうれしいです。

    補足日時：2016/11/20 14:50

• 加熱していない上辺と下辺は断熱です。表計算ソフトでチャレンジしてみたんですが計算式するための式が分からず、とまっています。良かったら表計算ソフトで計算に使った式を教えて頂けたらうれしいです。

  No.3の回答に寄せられた補足コメントです。 補足日時：2016/11/20 14:54

No.1 回答者： doc_somday 回答日時： 2016/11/19 19:15

これエクセルじゃ無理だしエクセルの微分方程式は力不足。

高価だがmathematicaを薦めます。

この回答へのお礼

回答ありがとうございます。
もし宜しければmathematicaでの計算式を教えて頂けたらうれしいです。

お礼日時：2016/11/21 10:34

No.2 回答者： NemurinekoNya 回答日時： 2016/11/19 21:31

この長方形の片側だけを加熱し、しかも、この長方形が半無限の場合は解析解が存在するけれど、この質問のケースについては解析的な解はないと思う。

あったとしても、その解は無限三角級数（フーリエ級数）の形で与えられるので、エクセルで関数の形で与えて解くことはできないと思う。

温度が300度にどんと上がってから時間が短い場合、半無限の解析解を使って計算することはできるけれど、１秒とか１０秒という長い時間の場合、この式は使えない。

ですから、エクセルを使ってこの問題を解こうとする場合、エクセルのプログラム記述用の言語VBAを使って、この偏微分方程式を数値的に解くことができるプログラムを自力で作り、その結果を使うということになる。

このためには、例えば、
「Excelで◯◯数値解析」
といった本を購入し、勉強をする必要がある。

こうした本には、１次非定常熱伝導方程式を解くプログラムが公開されていると思う。

ですが、
エクセルの表計算の機能を使って、陽解法で、この問題を解くことは可能。

以下のサイトなどが参考になると思います。
http://chemeng.in.coocan.jp/ExcelCe/exhc2.html

ワークシートまで公開されているので、数値を変更するだけで計算できるかもしれない。
ここで使っている解法は陽解法なので、ひょっとしたら、数値解が振動し、とんでもない結果になるかもしれないですね。
この場合、時間間隔を小さくして、数値的に解いた値が振動しないようにしないといけない。

この程度のワークシートならば、自力で簡単にできると思いますので、チャレンジしてみてください。
実は、もっと簡単に解く方法があるんだけどね。

No.3 回答者： NemurinekoNya 回答日時： 2016/11/19 22:32

ちょっと質問しますが、加熱していない両側は断熱ですよね。

ここで、例えば、対流や輻射などによって外部と熱のやりとりがあるわけじゃないですよね。
だとしたら、話がまったく違ってきます。

いま、ちょこっと表計算ソフトを計算してみたところ、計算は可能ですね。
簡単にできるので、チャレンジしてください。

わたしは、「教えてgoo」で回答をしないことにしているのですけれど、乗りかかった船ですから、最後まで面倒をみます。
わからないところがありましたら、遠慮せずに聞いてください。

この回答へのお礼

丁寧な回答ありがとうございます。
補足させて頂けたので良かったら回答お願いします。

お礼日時：2016/11/20 14:53

No.4 回答者： NemurinekoNya 回答日時： 2016/11/20 16:04

例えば、表計算ソフトで次のような表を作る。

ちなみに、B1からL1は計算点のx座標。
B2からL2までは時刻t=0の計算点xにおける温度。
B3以降は、時刻、時刻に応じた計算点の温度。
計算するときには、初期値として、B2、L2は300、C2〜K2は0を入力しておかないといけない。
また、B列のB２〜Bホニャララ、L列のL2〜Lホニャララまでには境界条件として300を事前にセットしておかないといけない。

計算式は、添付した写真に出ているけれど、
=(1-2*$O$7)*C2 +$O$7*(B2+D2)
これをC３からK３までにコピーぺして、
その後、B３からK３までを、B４〜K4までをBホニャララ〜Kホニャララまで計算できるはずです。

それらしい値が出ていますから、たぶん、間違ってはいないと思いますが、計算式は確かめてください。

添付した計算結果は、h=0.002mだから、2mm間隔ということになります。
h=0.0005だとセルの数が膨大になるので、表計算ソフトでの計算は厳しいかと思います。
この計算法ですと、θ=κΔt/h²<1/2にしないと、一般的に数値解が振動するので、θ<1/2におさめるようにΔｔを設定し直す必要があります。

プログラミング言語を何かご存知でしたら、それを使ったほうが賢明かと思います。

先に紹介したものは陽解法と呼ばれるもので、連立方程式を解く陰解法を使うと、θ<1/2という制約条件はなくなり、無条件安定で計算できますが、このためには連立方程式を解くプログラムを作らないといけない。

No.5 回答者： NemurinekoNya 回答日時： 2016/11/20 16:29

画像データが見づらいので、こちらをご覧になってください。

https://goo.gl/vpgTEv

No.6 回答者： NemurinekoNya 回答日時： 2016/11/20 17:26

この問題は、両端が３００℃で対称なので、この対称性を考慮すると、計算量が半分になります。

で、この対称性を考慮し、先のワークシートを少し変更し、1mm間隔、つまり、h=0.001ｍで計算した結果は次のようになるようです。
陽解法ですと、時間間隔Δtを小さくしないと、解が振動して大変なことになりますね(^^ゞ

参考にしてみてください。
それっぽい値が出ていますから、間違ってはいないと思います。
表計算ソフトの基本的な機能を使うだけでも、結構、いろいろなことができますね。
驚いております。

No.7 回答者： NemurinekoNya 回答日時： 2016/11/20 20:19

面倒くさいので、スプレッドシートを作って公開しました。

計算結果だけならば、h=0.001m、1mm間隔で、時間は0.025秒の場合のものを、以下のリンク先で閲覧可能です。
https://docs.google.com/spreadsheets/d/1eyGcSD-- …

以下のアドレスをクリックすると、エクセル形式のファイルがダウンロードされます。
https://docs.google.com/spreadsheets/d/1eyGcSD-- …

わたしはエクセルとは違う表計算ソフトを使っているので互換性が完全ではないようで、フィルを変換した時に、計算結果のグラフの表示が変わってしまっているようですが、エクセルでも計算は可能だと思います。

対称性を考慮して計算させているので、計算領域は加熱面から中心部までになっています。

マクロは一切使っていないので、ウィルス感染の心配は無用かと思いますが、ファイルをダウンロードするときには覚悟してください。

それっぽい値が出ていますから、計算は間違っていないと思いますが・・・。