Scilabで熱拡散方程式 その1 (陽解法)

Scilabには常微分方程式ソルバや非線形方程式ソルバは存在しますが、偏微分方程式ソルバというものは存在しません。
そこで今回はOctaveの精義をを参考にして、偏微分方程式の最も簡単な例の一つである、一次元の熱伝導の問題を計算しました。

\frac{\partial u}{\partial t} = \frac{\partial^2 u}{\partial x^2}


偏微分方程式


Scilabを利用すると色々な数値計算を行うことが出来ます。
この中で常微分方程式非線形方程式を解くのは比較的簡単でScilabがデフォルトで持っている命令をワンパターンで使うだけです。

これに対して偏微分方程式をワンパターンで解けるような命令はScilabには用意されていません。今回はOctaveの精義を参考に一次元の熱伝導を陽解法(FTCS法:Forward Time Space methods)で解いてみます。

熱拡散方程式


物質の拡散や熱伝導は熱拡散方程式と呼ばれる偏微分方程式で表されます。
今回扱う一次元の物体内部の温度分布は以下のような形になります。

\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

するとxの範囲が 0 ≦ x ≦ 1 となります。

またτは

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

すなわち

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

です。

温度も計算結果の値が大きすぎたり小さすぎたりしないように適切なT0を選んで

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

の様に無次元化します。

すると偏微分方程式は以下のような形になります。

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

0 ≦ x ≦ 1, t ≧ 0

無次元化によって

T' → u
x' → x
t' → t

と3つの変数が変更されました。
以降ではこの無次元化された偏微分方程式について数値計算を行います。

方程式の差分化


更に前述の無次元化された偏微分方程式の微分を差分に変換します。

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

を差分化すると以下の様になります。

\frac{\partial u(x,t)}{\partial t} \sim \frac{u(x,t+\Delta t)-u(x,t)}{\Delta t}

\frac{\partial^2 u(x,t)}{\partial x^2} \sim \frac{u(x+\Delta x,t)-2u(x,t)+u(x-\Delta x, t)}{(\Delta x)^2}

式が長くなるので、以下の様に ui,j を定義します。

u_{i,j} \equiv u((i-1)\Delta x, (j-1)\Delta t)

結局

\frac{u_{i,j+1}-u_{i,j}}{\Delta t} = \frac{u_{i+1,j}-2u_{i,j}+u_{i-1,j}}{(\Delta x)^2}

ここで

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

とおくと

u_{i,j+1} = (1-2r)u_{i,j} + r(u_{i+1,j}+u_{i-1,j})

これには式の見た目が簡単になるということ以上に重要な意味があります。
rが小さいほど計算が高速になるのですが r < 1/2 では計算に失敗し正しい結果が得られないことが知られています。また r = 1/6 のときだけ打切り誤差が特別に小さくなるという性質があることも知られています。従って実際の計算では r = 1/2 または r = 1/6 のどちらかが採用されます。

この事は裏を返せば空間分解能を上げたいときは、同時に時間分解能を上げなければならない(逆も然り)ということでもあります。

Scilabでガウス型波束の散乱でハミルトニアンを行列で表したのと同様に

7861cb402820e6f58fddc8d3abb3e335_90_black.png

となる様な行列Pを考えると

f99827543db336a6148d03e0d830d8ab_90_black.png

ただし

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

s = 1 - 2r

です。

初期条件と境界条件


常微分方程式を解くために初期条件が必要であったのと同様に、偏微分方程式を解くためには初期条件と境界条件が必要になります。
問題設定によって色々な境界条件が考えられ、その境界条件に応じて異なったプログラムを書かなければならない事が偏微分方程式の難しさの元になっているのだと思うのですが、今回は入門書でよく見る以下のような簡単なものを考えます。

初期条件(t=0)
u(x,0) = f(x)

境界条件(x=0, x=1)
u(0,t) = 0
u(1,t) = 1

この初期条件と境界条件は、次のような状態を表していると考えることが出来ます。

まず、金属の棒全体を0℃の状態にしておく(初期条件)。
次に左端を0℃の氷に接しさせ、右端を100℃のお湯に接しさせる(境界条件)。

これを時間発展させれば、金属の棒の内部の温度分布がどうなるかを調べられます。
非常に長い時間を置いた後の最終的な温度分布は、直線的になるであろうことは、計算しなくても直感的に予想が付きます。

計算結果


以下に計算結果を示します。
各線は 200×Δt の時間ごとに結果をプロットしたもので一番最後のものが 1600×Δt の時間のものです。
この段階では温度分布がほとんど直線的になっており、直感的な予想が正しかったことが分かります。


001_20140119171908a1e.png
Fig.1: 計算結果


clear;

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

r = 1 / 6;

dt = r * dx ^ 2; // 時間の刻み幅
m = 1601; // 時間の分割数
dm = 200; // プロットする m の間隔

// *** 初期条件と境界条件 ***
// 初期条件
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

// *** グラフのプロット ***
plot([0:dx:1],U(:,1:dm:m));


関連エントリ




参考URL




付録


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


参考文献/使用機器






フィードバック



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

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


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


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


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

comment

Secret

管理人のみ閲覧できます

このコメントは管理人のみ閲覧できます

Re: 2012-05-14

> タイトル2012-05-14の中ほどの文章が主さんのおっしゃりたい意味合いと異なってるように感じます
> ”従ってからの”表記としては”直流解析では~”だと思うのですが

ご指摘ありがとうございます。下記エントリの内容を修正いたしました。

LTspiceで直流解析(DC解析)
http://gomisai.blog75.fc2.com/blog-entry-511.html

このようなミスの指摘は大歓迎です。
FC2カウンター
カテゴリ
ユーザータグ

LTspiceAkaiKKRmachikaneyamaScilabKKRPSoCOPアンプPICCPA強磁性常微分方程式モンテカルロ解析odeトランジスタ状態密度インターフェーススイッチング回路ecaljPDS5022DOS定電流半導体シェルスクリプト乱数レベルシフトHP6632Aブレッドボード分散関係温度解析R6452Aトランジスタ技術I2C可変抵抗反強磁性セミナー数値積分確率論偏微分方程式バンド構造非線形方程式ソルババンドギャップ絶縁熱設計シュミットトリガLEDA/Dコンバータ三端子レギュレータLM358ISO-I2CGW近似カオスフォトカプラマフィンティン半径TL431数値微分PC817Cアナログスイッチ直流動作点解析発振回路USBサーボカレントミラー74HC4053パラメトリック解析LDAbzqltyチョッパアンプ量子力学FFT2ちゃんねるアセンブラBSch開発環境電子負荷ブラべ格子イジング模型補間基本並進ベクトル標準ロジック単振り子キュリー温度繰り返しMaxima状態方程式失敗談相対論スピン軌道相互作用FETランダムウォーク熱伝導六方最密充填構造コバルトewidthTLP621GGAQSGW不規則合金位相図抵抗SMPcygwinラプラス方程式スレーターポーリング曲線gfortranスイッチト・キャパシタ詰め回路TLP552三角波格子比熱TLP521条件分岐LM555MCUNE555QNAPマントルテスタ過渡解析FXA-7020ZRダイヤモンドデータロガーガイガー管自動計測Writer509UPSシュレディンガー方程式ブラウン運動awk差し込みグラフ熱力学平均場近似仮想結晶近似VCAfsolve井戸型ポテンシャルVESTA起電力スーパーセルOpenMP第一原理計算ubuntu固有値問題L10構造OPA2277interp12SC1815fccウィグナーザイツ胞面心立方構造フィルタジバニャン方程式ヒストグラム確率論マテリアルデザインspecx.f等高線正規分布PGAフェルミ面非線型方程式ソルバ初期値固定スピンモーメントスワップ領域ルチル構造リジッドバンド模型edeltquantumESPRESSO岩塩構造BaOSIC二相共存ZnOウルツ鉱構造フォノンデバイ模型c/aノコギリ波全エネルギーFSMTeXgnuplotmultiplotハーフメタルCapSense半金属合金結晶磁気異方性Ubuntu文字列入出力TS-110TS-112疎行列Excel直流解析ヒストグラム円周率不規則局所モーメントトラックボールPC等価回路モデルパラメータ・モデルキーボードRealforce三次元マンデルブロ集合フラクタル化学反応重積分縮退日本語最小二乗法関数フィッティングGimpMAS830LHiLAPW熱拡散方程式両対数グラフナイキスト線図負帰還安定性陰解法Crank-Nicolson法P-10クーロン散乱境界条件連立一次方程式片対数グラフEAGLEPIC16F785LMC662トランスシンボルCK1026線種凡例MBEAACircuitグラフの分割軸ラベルifort

最新コメント
リンク

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