Scilabで熱拡散方程式 その2 (無次元化)

Scilabで熱拡散方程式 その1 (陽解法)ではOctaveの精義を参考にして以下に示す一次元の熱拡散方程式を無次元化してScilabで計算しました。

\frac{\partial T'}{\partial t'} = D^2 \frac{\partial^2 T'}{\partial {x'}^2}

ただし 0 ≦ x' ≦ λ, t' ≧ 0

(T': 温度, t': 時間, D2: 熱拡散率, x': 位置)

今回はもう少し具体的に、現実に即した問題を無次元化し、更にもう一度具体的な次元を持った値に直すという手順をやってみます。


問題設定


Scilabで熱拡散方程式 その1 (陽解法)では一次元の熱伝導の問題を解くために熱拡散方程式を無次元化して陽解法でで計算を行いました。
無次元化を行ったということは、計算結果は実際の長さや時間とは違う目盛りで書かれているということなので、無次元化とは逆の処理をして実際のものに直す必要があります。

今回は以下のような具体的な条件でScilabで熱拡散方程式 その1 (陽解法)と全く同じ計算を行います。

長さ1.5mの銅でできた棒の左端を0℃の氷に、右端を100℃のお湯に接しさせたとき、1時間後の温度分布はどのようになるか。ただし、最初の温度分布は全て0℃であったとする。

銅の熱拡散率D2 = 116 × 10-6 m2/sは以下の計算式とパラメータから計算しました。

D^{2} = \frac{k}{\rho C_p}

(D2: 熱拡散率, ρ: 密度, CP: 定圧比熱)

パラメータ
熱伝導度 k398 W/m/K
密度 ρ8.92×103 kg/m3
定圧モル比熱 CP,m24.44 J/mol/K
原子量 M6.3546×10-2 kg/mol
table.1: 熱拡散率D2=116×10-6 m2/s を計算するためのパラメータ。定圧比熱は密度と単位をそろえるためにCP=CP,m/Mとする。各パラメータの出展はWikipedia。


分子のおもちゃ箱 熱拡散長には色々な金属の熱拡散率が表にまとめられています。

熱拡散方程式


Scilabで熱拡散方程式 その1 (陽解法)と全く同じ内容をもう一度書いておきます。今回解くべき偏微分方程式は以下のものになります。

\frac{\partial T'}{\partial t'} = D^2 \frac{\partial^2 T'}{\partial {x'}^2}

ただし 0 ≦ x' ≦ λ, t' ≧ 0

(T': 温度, t': 時間, D2:熱拡散率, x': 位置)

なお文字の上についている'(アポストロフィー)は(後で行う)無次元化をしていないことを意味します(微分を意味するものではありません)。

方程式の無次元化


以下の様に方程式の無次元化することを考えます。

x' = λx
t' = τt
T'(x',t') = T0u(x,t)

ただし

\frac{D^2 \tau}{\lambda^2}=1

の関係に注意が必要です。

いま銅の棒の長さは1.5mなので λ=1.5 m です。
これに熱拡散率 D2 = 116 × 10-6 m2/s をあわせてτが必然的に決まります。

\tau = \frac{\lambda^2}{D^2}

計算結果は τ = 19394.442 s となりました。

最後に温度の無次元化ですが、これは単純に T0 = 100 ℃ と置けばよいでしょう。このように置けば 0 ≦ T' ≦ 100 に対して 0 ≦ u ≦ 1 となります。

結局、無次元化した熱拡散方程式は下記の様になります。

\frac{\partial u(x,t)}{\partial t} = \frac{\partial^2 u(x,t)}{\partial x^2}

0 ≦ x ≦ 1, t ≧ 0

方程式の差分化


無次元化された熱拡散方程式を差分化するに際して、時間と空間の刻み幅をまとめて

r = \frac{\Delta t}{(\Delta x)^2}

と置いたわけですが、これが r < 1/2 では計算が失敗し r = 1/6 のときだけ誤差が小さくなるという特徴があると前回書きました。

ここからは r = 1/6 を選ぶとして話を進めます。rにこれ以外の数字を選ぶにしても、必ず1/2 未満にならないようにチェックをする必要があります。

いま一次元空間 0 ≦ x ≦ 1 の範囲でn点ほど計算するとするならば、Δxの大きさは Δx = 1/(n-1) となります。
rとΔxが決まるとΔtも必然的に決まり Δt = r×(Δx)2 です。

時間はτを使って無次元化されているため、無次元化された時間の刻み幅Δtもτを使って実際の時間の刻み幅に戻すことが出来ます。

\tau \Delta t = \tau r (\Delta x)^2 = \frac{\tau r}{(n-1)^2}

仮に空間の計算点数を n=50 とすると、τΔt= 1.3462753 s となります。
今回の計算では1時間後(60×60=3600秒後)の温度分布を計算するので m ≒ 3600 / (τΔt) ≒ 2674 回の繰り返し計算をします。

計算結果


以下に示すのが 約1時間後(3599.9401秒後)の銅の棒の温度分布です。


001_20140123132748b05.png
Fig.1: 全体が0℃であった銅の棒の左端を0℃に固定、右端を100℃に固定した状態で約1時間(3599.9401秒)保持したときの内部の温度分布。


計算の手順はこのエントリに書いたとおりです。

clear;

// *** 銅の物性値 ***
k = 398; // 熱伝導度 398 W/m/K
rho = 8.92E3; // 密度 8.92×10^3 kg/m^3

cpm = 24.44; // 定積モル比熱 24.44 J/mol/K
m = 63.546; // 原子量 g/mol
mkg = m * 1E-3; // 原子量 kg/mol
cp = cpm / mkg; // 定積質量比熱 J/kg/K

d2 = k / (rho * cp) // 熱拡散率 m^2/s

// *** 無次元化前の時間と空間 ***
lambda = 1.5; // 棒の長さ λ (m)
tdmax = 60 * 60; // t'max (s) まで計算

// *** 無次元化 ***
tau = lambda ^ 2 / d2; // 時間の無次元化係数
tmax = tdmax / tau; // 無次元化時間で tmax まで計算
temp0 = 100; // u=1 に対してT'=100℃

// *** 時間と空間の分割 ***
n = 50; // 空間の分割数
dx = 1 / (n - 1); // 空間の刻み幅

r = 1 / 6;

dt = r * dx ^ 2; // 時間の刻み幅
m = round(tmax / dt); // 時間の分割数

// *** 初期条件と境界条件 ***
// 初期条件
U = zeros(n,m);
// 境界条件
U(1,1) = 0;
U(n,1) = 1;

// *** 行列P ***
s = 1 - 2 * r;
P = toeplitz([s, r, zeros(1, n - 2)]);

// *** 偏微分方程式の計算 ***
for j = 1:(m - 1)
// 境界以外の計算
U(:, j + 1) = P * U(:,j);
// 境界条件
U(1, j + 1) = 0;
U(n, j + 1) = 1;
end

// *** グラフのプロット ***
X = linspace(0,lambda,n);
plot(X,temp0 * U(:,m),'o-b');
zoom_rect([0,0,lambda,temp0]);

xlabel("Position (m)");
ylabel("Temperature (degree C)");


関連エントリ




参考URL




付録


このエントリで使用したScilabのシミュレーション用ファイルを添付します。ファイル名末尾の".txt"を削除して、"_"を"."に変更すれば使えるはずです。(参考:ねがてぃぶろぐの付録)


参考文献/使用機器







フィードバック



にほんブログ村 その他趣味ブログ 電子工作へ

 ↑ 電子工作ブログランキング参加中です。1クリックお願いします。


コメント・トラックバックも歓迎です。 ↓      


 ↓ この記事が面白かった方は「拍手」をお願いします。


tag: Scilab 偏微分方程式 熱拡散方程式 熱伝導 

comment

Secret

FC2カウンター
カテゴリ
ユーザータグ

LTspiceAkaiKKRmachikaneyamaScilabKKRPSoC強磁性OPアンプPICCPA常微分方程式モンテカルロ解析ecaljodeトランジスタ状態密度インターフェースDOS定電流スイッチング回路PDS5022半導体シェルスクリプトレベルシフト乱数HP6632AR6452AI2C可変抵抗分散関係トランジスタ技術ブレッドボード温度解析反強磁性確率論バンドギャップセミナー数値積分熱設計非線形方程式ソルババンド構造絶縁偏微分方程式ISO-I2CLM358フォトカプラ三端子レギュレータカオスLEDシュミットトリガGW近似A/Dコンバータ発振回路PC817C直流動作点解析USBマフィンティン半径数値微分アナログスイッチTL43174HC4053カレントミラーサーボ量子力学単振り子チョッパアンプ補間2ちゃんねる開発環境bzqltyFFT電子負荷LDAイジング模型BSch基本並進ベクトルブラべ格子パラメトリック解析標準ロジックアセンブラ繰り返し六方最密充填構造SMPコバルトewidthFET仮想結晶近似QSGW不規則合金VCAMaximaGGA熱伝導cygwinスレーターポーリング曲線キュリー温度スイッチト・キャパシタ失敗談ランダムウォークgfortran抵抗相対論位相図スピン軌道相互作用VESTA状態方程式TLP621ラプラス方程式TLP552条件分岐NE555LM555TLP521マントル詰め回路MCUテスタFXA-7020ZR三角波過渡解析ガイガー管自動計測QNAPUPSWriter509ダイヤモンドデータロガー格子比熱熱力学awkブラウン運動起電力スーパーセル差し込みグラフ第一原理計算フェルミ面fsolve最大値xcrysden最小値最適化ubuntu平均場近似OpenMP井戸型ポテンシャルシュレディンガー方程式固有値問題2SC1815結晶磁気異方性OPA2277非線型方程式ソルバTeXgnuplot固定スピンモーメントFSMPGAc/a全エネルギーfccフラクタルマンデルブロ集合正規分布縮退初期値interp1multiplotフィルタ面心立方構造ウィグナーザイツ胞L10構造半金属二相共存SICZnOウルツ鉱構造BaO重積分クーロン散乱磁気モーメント電荷密度化学反応CIF岩塩構造CapSenseノコギリ波デバイ模型ハーフメタルキーボードフォノンquantumESPRESSOルチル構造スワップ領域リジッドバンド模型edelt合金等高線線種凡例シンボルトラックボールPC軸ラベルグラフの分割トランス文字列CK1026MAS830L直流解析Excel不規則局所モーメントパラメータ・モデル入出力日本語最小二乗法等価回路モデルヒストグラムGimp円周率TS-110TS-112PIC16F785LMC662三次元specx.fifortUbuntu疎行列不純物問題Realforceジバニャン方程式ヒストグラム確率論マテリアルデザインP-10境界条件連立一次方程式熱拡散方程式AACircuitHiLAPW両対数グラフ片対数グラフ陰解法MBEナイキスト線図負帰還安定性Crank-Nicolson法EAGLE関数フィッティング

最新コメント
リンク

にほんブログ村 その他趣味ブログ 電子工作へ