定常１次元熱伝導方程式の数値計算

定常１次元熱伝導方程式の数値計算

熱伝導率λが温度Tなどによって変化する場合、非定常の１次元熱伝導方程式は次のようになる。

ここで、ρは密度、cは比熱、は単位時間単位体積あたりに発生（消滅）する熱量である。

したがって、熱の発生がない場合、１次元の定常熱伝導方程式は

となり、境界条件が与えられれば、この解を求めることができる。

たとえば、熱伝導率λが一定で、x=0における温度がT₀、x=lにおける温度がT₁の場合、

さらに、単位時間単位面積当たりに通過する熱量（熱流束）は

である。

熱伝導率が一定の場合、温度Tは直線的に変化するので、わざわざ数値計算をする必要はない。

そこで、熱伝導率λが温度によって変化する場合の定常１次元熱伝導方程式の差分法を用いた数値解法について考えてみる。

（２）は

となるので、が

の場合、差分法を用いて

と近似することが可能。

しかし、プログラムがすこし複雑になるので、今回、この計算法は採用しないことにする。

そこで、今回は、

という差分を用い、（２）式を

ここで、とは、それぞれ、との中点における熱伝導率。

それで、

 

そして、熱伝導率λは

として、これを解くプログラムを作ってみた。

ちなみに、この微分方程式の解は

である。

計算領域0≦x≦1をn=１０分割した計算結果は、次の通り。

青い直線は熱伝導率が一定の場合。

計算に使用したプログラムは次の通り。

parameter (n=10)

real t(0:n), gam(0:n)

real a(n),b(n),c(n),d(n)

eps=1.e-6

dx =1./n

! 変数の初期化

t=0.; gam=1.

a=0.; b=0.; c=0.; d=0.

! 境界条件

t(n)=1.

do k=1, 10

! 熱伝導率Γを計算

do i=1, n

gam(i)=1+t(i)

end do

! 連立方程式の係数を計算

do i=1, n-1

dx2=dx*dx

aw=0.5*(gam(i-1)+gam(i))/dx2

ae=0.5*(gam(i)+gam(i+1))/dx2

ap=aw+ae

a(i)=-aw

b(i)=ap

c(i)=-ae

d(i)=0.

end do

a(1)=0.; d(1)=0.5*(gam(0)+gam(1))/dx2*t(0)

c(n-1)=0.; d(n-1)=0.5*(gam(n-1)+gam(n))/dx2*t(n)

! TDMAで連立方程式を解く

call tdma(a,b,c,d,n-1)

! 計算結果の更新と収束判定

err=0.0

do i=1,n-1

err=amax1(err,abs(1-t(i)/d(i)))

t(i)=d(i)

end do

if(err.lt.eps) exit

end do

do i=0,n

x=i*dx

write(*,100) x, t(i), -1+sqrt(1+3.*x)

end do

! ファイルへの出力

open(1, file='Netsu.dat', status='replace')

do i=0,n

x=i*dx

write(1,100) x, t(i),-1+sqrt(1+3.*x)

end do

close(1)

100 format(2(f8.5,1x),f8.5)

end

! TDMA

subroutine tdma(a,b,c,d,n)

real a(n), b(n), c(n),d(n)

! 前進消去

do i=1,n-1

ratio=a(i+1)/b(i)

b(i+1)=b(i+1)-ratio*c(i)

d(i+1)=d(i+1)-ratio*d(i)

end do

! 後退代入

d(n)=d(n)/b(n)

do i=n-1,1,-1

d(i)=(d(i)-c(i)*d(i+1))/b(i)

end do

end

なお、このプログラムでは、とを

と、点における熱伝導率の算術平均を用いて計算している。

熱伝導率のこの計算の是非については議論が必要だろうが、計算結果は極めて良好なのだから、「終わりよければ全てよし」ということで(^^)

本問は、いちおう、Tの非線形の微分方程式なので、計算においては、

と、反復１回前の温度を用いて熱伝導率を計算し、線形化して解いている。

反復計算の収束条件は