偏微分方程式の解法熱拡散問題

偏微分方程式の解法
熱拡散問題
図のような金属の棒がある。金属の両端の温度は
T=0 に設定されている
x=L/2
x=0
時刻 t=0 で
温度分布が右図で表されている
ような分布であったとすると
T
T ( x  0, t  0)  0
x=L
T=T0
T ( x  L, t  0)  0
時刻 t=t での温度分布はどうなっているか？
これを数値計算で求める。
偏微分方程式の解法(熱拡散問題）
偏微分方程式の解法（熱拡散方程
式）
熱の熱拡散は次の方程式で求めることが出来る。
T
T
 2
t
x
2
初期条件
x
T (x,0)  2T0
L
x

T (x,0)  2T0  1  
 L
(x  L / 2)
(x  L / 2)
偏微分方程式の解法（差分メッシュ）
T[0] T[1] T[2]
x0
x1
x=0
x軸を N 等分
x2
T[i]
T[N]
xi
x
xN
x=L
x  L / N
位置 xi  ix における温度を T[i] に記憶。
配列要素は N+1 個必要
double
棒の両端の温度は 0 に固定
T[N+1];
T[0]=0
T[N]=0
差分化
T
T
 2
t
x
2
T 1
  T ( x, t  t )  T ( x, t ) 
t t
T
1
 T ( x  x, t )  2T ( x, t )  T ( x  x, t ) 

2
2
x
( x )
2
T (x,t  t)  T (x,t)
t

T (x  x,t)  2T (x,t)  T (x  x,t) 
2
(x)
T[0], T[1], …, T[N] だけでは足りない。
各 i 毎に2個の温度を記憶する必要がある
配列の宣言は double T[N+1][2]; となる。
T[0][0], T[1][0], T[2][0], …, T[N][0]
T[0][1], T[1][1], T[2][1], …, T[N][1]
t
T (xi ,t  t)  T (xi ,t) 
T (xi 1,t)  2T (xi ,t)  T (xi 1,t) 
2
(x)
x0  0
xi 1  xi  x
xi
xi 1  xi  x
t
T[0][now] … T[i-1][now] T[i][now] T[i+1][now] …
t  t
T[0][new] … T[i-1][new] T[i][new] T[i+1][new] …
偏微分方程式の解法（差分メッシュ）
計算の進め方
step=0
(初期状態）
step=0 now=0 new=1
T[i][0] から T[i][1] を計算
step=1
step=2
t 0
t  t
step=3
step=1 now=1 new=0
T[i][1] から T[i][0] を計算
step=5
step=2 now=0 new=1
step=4
t  t
t  2t
T[i][0] から T[i][1] を計算
step=6
t  2t
now=0
now=1
new=1
new=0
t  3t
偏微分方程式の解法（差分メッシュ）
step=0
T(i,0)
step=2
T(i,0)
step=4
T(i,0)
step=1
T(i,1)
step=3
T(i,1)
now= step % 2
new= (step+1) % 2 とすると
step=0
now=0
step=1
now=1
step=2
now=0
…..
が実現出来る。
new=1
new=0
new=1
T (x,t  t)  T (x,t)
t

T (x  x,t)  2T (x,t)  T (x  x,t) 
2
(x)
a = kappa*dt/(dx*dx);
T[i][new] = T[i][now]
+ a*(T[i+1][now]-2*T[i][now]+T[i-1][now]);
#include <stdio.h>
#define N 10
/* x 軸の分割数。定数でなければならない */
int main()
{
double T0, L, dx;
double kappa, dt, tmax, a;
double T[N+1][2];
int i, now, new;
int step, nstep;
int nout, m;
/* 計算条件 */
T0 = 100.0;
/*
L = 1.0;
/*
kappa = 1.0;
/*
dt = 0.002; tmax = 0.5; /*
nout = 10;
/*
dx = L/N;
a = kappa*dt/(dx*dx);
nstep = tmax/dt;
m = nstep/nout;
(続く)
中心の温度 */
棒の長さ */
熱拡散係数 */
時間刻み幅と最終時刻 */
出力回数 */
/* x軸の刻み幅 */
/* 計算回数 */
/* 出力間隔 */
/* 初期温度を設定 */
now = 0; new = 1;
for(i=0; i<N/2; ++i) {
T[i][now] = i*2.0*T0/N;
}
for(i=N/2; i<=N; ++i) {
T[i][now] = (N-i)*2.0*T0/N;
}
T[0][new] = 0.0; T[N][new] = 0.0;
/* 境界条件: 両端は常に T=0 */
/* 各時刻の温度を計算 */
for(step=0; step<nstep; ++step) {
/* m ステップ毎に温度を出力 */
if(step % m == 0) {
for(i=0; i<=N; ++i) {
printf("%f\t%f\n", i*dx, T[i][now]);
}
printf("\n\n");
/* 空行を2行出力 */
}
/* now と new の設定 */
now = step % 2;
new = (step+1) % 2;
/* 各点の新しい温度を計算 */
for(i=1; i<N; ++i) {
/* 両端 i=0,N は計算しない */
T[i][new] = T[i][now]
+ a*(T[i+1][now] - 2.0*T[i][now] + T[i-1][now]);
}
}
return 0;
}
gnuplot に複数のグラフを一度に描かせる
1.0
2.0
3.0
1.0
2.0
2.5
3.0
1.2
1.8
2.2
1.7
1.9
2.3
2.8


1つめのグラフのデータ
空行2行以上
2つめのグラフのデータ

Download Report