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

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

 

熱伝導率λが温度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の非線形の微分方程式なので、計算においては、

　　

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

反復計算の収束条件は