エクセルで2次元熱伝導方程式を解く

Tweet

Pocket

どうもYukiです。
エクセルで1次元熱伝導方程式を解くについて書いてきましたが、今回はエクセルを使って実際に2次元熱伝導方程式を解いていきたいと思います。

計算条件

まず、熱伝導方程式を解くために計算条件を設定します。
定常計算、発熱あり 、2次元の熱伝導方程式、境界条件は左側1000[K]、上側800[K] 、下側400[K]で縛り、右側はフリーにします。

熱伝導率k = 80.3 [W/m・K]、ΔX = 1 [m]、発熱は中央のみ（1000[W]）初期温度は0[K]、計算点数は5×5点で中心差分法を使用します。

0 = k(∂²T/∂x² +∂²T/∂y² ) + q′_v…　2次元の熱伝導方程式

計算手順

条件をプロットするとこのようになります。

プロットするとこのような表になると思います。

言葉で詳しく説明すると、境界条件（一番端っこの所）より1 , 5点目に1000 , 800 , 400を入力します。
右側はfree（吹き抜け）設定、中央は発熱するので発熱設定です。

片側差分法

Free条件にするために片側差分法というものを使います。
前回の記事に中心差分法について書きました。
この差分法は求めたい場所の前後を使って計算します。
片側差分法は求めたい場所の前後どちらかを使って求めます。
差分式は

dƒ/dx    = (ƒ(i) – ƒ(i-1)) / 2Δx

d²ƒ/dx² = (ƒ(i) – 2ƒ(i-1) + ƒ(i-2)) / Δx²

dƒ/dx    = (ƒ(i) – ƒ(i+1)) / 2Δx

d²ƒ/dx² = (ƒ(i) – 2ƒ(i+1) + ƒ(i+2)) / Δx²

です。

表でにすると

ここの値を使ってFreeを求めます。

補足

二次元なのでFreeの部分も2次元（縦方向）にも差分を取らないといけないが、収束しなかったので1次元の片側差分法を使いました。

この問題なら外挿（四点目）の値をFreeの部分（5点目）にそのまま挿入する方法のあります。

計算手順2

次に2～4点目に中心差分をした2次元の熱伝導方程式を入れます。

となるはずです。

もう少し詳しく説明すると差分を入れた、2次元の熱伝導方程式は

T(i,j) = 0.5(T(i-1,j) + T(i+1,j) +

T(i,j-1) + T(i,j+1))

になります。

表にするとこんな風に求めることが出来ます。

これで2次元熱伝導方程式が解けるはずです。

グラフの作り方

2次元熱伝導方程式のグラフの作り方、正直面倒くさいです。

等高線図を使えば何とか、それらしい物が出来ます。

こんな感じで作れます。
これを作るのに2,3時間位かかりました。

そこで、グラフではありませんが、条件式書式によってセルに色を塗り、温度分布を表現しようと思います。

まず、色を付けたいセルを選択し、ホームの条件式書式の中の新しいルールをクリックします。

すると下の図のような画面が出ると思います。

書式スタイルを三色スケールにし中間値の種類を数値にし、値を任意に入れます。（下図では450）

入れ終わったら、OKを押します。

すると下の図になります。

このようにして、グラフに近いものが出来ました。
数字が邪魔な場合はセルを小さくする等工夫すると消えます。

これが一番簡単に温度分布を表す方法だと思います。

まとめ

今回は一つの境界条件がFreeなので片側差分法と中心差分法を併用しました。

Freeが入っているので条件設定では、発散する場合があります（1回発散しました（笑））

気を付けましょう

Fortranのプログラムコードも載せておきます。

プログラムコード

PROGRAM MAIN
INCLUDE ‘implicit.for’
C———————————————— Global Variable —–
REAL*4 X_T(ISSS:IEEE,JSSS:JEEE)        !— Position X-component
REAL*4 Y_T(ISSS:IEEE,JSSS:JEEE)        !— Position Y-component
REAL*4 T(ISSS:IEEE,JSSS:JEEE)          !— Thermal
REAL*4 T0 (ISSS:IEEE,JSSS:JEEE) !— Temperature for Meomory
REAL*4 T1                       !— Upeer Temperature
REAL*4 T2                       !— Left Temperature
REAL*4 T3                     !— Bottom Temperature
C———————————————————————-

C
INTEGER I0, I1, J0, J1 ,N0 ,N1   !— Mesh Size
INTEGER S0, S1                   !— Mesh Size

C
REAL*4 RES_P                     !— Residual P
REAL*4 RES_PP                    !— Residual PP
REAL*4 RES_PPP                   !— Residual PPP
REAL*4 RES_T                     !— Residual T
REAL*4 RES_E                     !— Residual E
C
REAL*4 DELY_T                     !— deltaY T
C
INTEGER NMAX_P, NMAX_T, NMAX_Y  ,ALL ,NMAX_E

C————————————————- Local Variable —–
REAL*4  DX , DY
INTEGER I,J,S ,BB
INTEGER N
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
110 FORMAT(A)
WRITE(*,110) ‘**************************************************’
WRITE(*,110) ‘************                          ************’
WRITE(*,110) ‘************     Main Program         ************’
WRITE(*,110) ‘************      Start !!!     ************’
WRITE(*,110) ‘************                          ************’
WRITE(*,110) ‘**************************************************’
WRITE(*,110) ‘ ‘
C
111 FORMAT(A,E16.8E3,A)
116 FORMAT(A,I5,10X,I4)
C
WRITE(*,110) ‘———– Reading Parameter of Mesh ————‘
OPEN  (11,FILE=’./para_mesh.txt’,STATUS=’unknown’)
READ(11,*) I0
READ(11,*) I1
READ(11,*) J0
READ(11,*) J1
READ(11,*) S0
READ(11,*) S1
CLOSE (11)
WRITE(*,110) ‘———————————— has done. —‘
WRITE(*,110) ‘ ‘
WRITE(*,110) ‘>>> Grid Parameter — start ——- end ———-‘
WRITE(*,116) ‘    I-Direction :   ‘,I0,I1
WRITE(*,116) ‘    J-Direction :   ‘,J0,J1
WRITE(*,110) ‘                                                 ‘
c      GOTO 999

NMAX_T = 10000
C—————————–Position X,Y——————————-
DO J = J0, J1
DO I = I0, I1
X_T(I,J) = 5/ 5 * REAL(I)
Y_T(I,J) = 5/ 5 * REAL(J)
END DO
END DO

C—————————–initial——————————-
T1= 800
T2= 1000
T3= 400
DO I = I0+1, I1-1
DO J = J0+1, J1-1
T(I,J) = 500
END DO
END DO
C
DO J=J0,J1
DO I=I0,I1
T(I,J0) = T3
T(I0,J) = T2
T(I,J1)= T1
END DO
END DO
C

C—– Computation of Thermal Eq. ————————————
DO BB = 1 , NMAX_T
CALL MEM_T   (I0,I1,J0,J1,T,T0)
C
CAll CAL_T   (I0,I1,J0,J1,T,T0,X_T,Y_T)
C
write(*,*)’T(2,2)’,T(2,2)
CAll BOU_T   (I0,I1,J0,J1,T,T1,T2,T3)
C
CAll CAL_R_T (I0,I1,J0,J1,T,T0,RES_T)

CALL OUTPUT_T1 (I0,I1,J0,J1,X_T,Y_T,T)
C
END DO

WRITE(*,110) ‘**************************************************’
WRITE(*,110) ‘************                          ************’
WRITE(*,110) ‘************      Main Program        ************’
WRITE(*,110) ‘************       End !!!      ************’
WRITE(*,110) ‘************                          ************’
WRITE(*,110) ‘**************************************************’
WRITE(*,110) ‘ ‘
C
999 STOP
END PROGRAM
C *********************************************************************
SUBROUTINE MEM_T (I0,I1,J0,J1,T,T0)
INCLUDE ‘implicit.for’
C———————————————— Global Variable —–
REAL*4 T(ISSS:IEEE,JSSS:JEEE)          !— Thermal
REAL*4 T0 (ISSS:IEEE,JSSS:JEEE) !— Temperature for Meomory
INTEGER I0, I1, J0, J1           !— Mesh Size
C————————————————- Local Variable —–
INTEGER I,J
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
DO J = J0, J1
DO I = I0, I1
T0(I,J) = T(I,J)
END DO
END DO
C
999 RETURN
END
C
C
C
C *********************************************************************
SUBROUTINE CAL_T (I0,I1,J0,J1,T,T0,X_T,Y_T)
INCLUDE ‘implicit.for’
C———————————————— Global Variable —–
REAL*4 T(ISSS:IEEE,JSSS:JEEE)          !— Thermal
REAL*4 T0 (ISSS:IEEE,JSSS:JEEE) !— Temperature for Meomory
REAL*4 X_T(ISSS:IEEE,JSSS:JEEE)  !— Position X-component
REAL*4 Y_T(ISSS:IEEE,JSSS:JEEE)  !— Position Y-component
INTEGER I0, I1, J0, J1           !— Mesh Size

C————————————————- Local Variable —–
INTEGER I,J,S
REAL*4 DX, DY
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
DX = X_T(1,0) – X_T(0,0)
DY = Y_T(0,0) – Y_T(0,1)
C
DO J = J0+1, J1-1
DO I = I0+1, I1-1
T(I,J) = 0.5 / (DX**2+DY**2)
&         *(T0(I-1,J)+T0(I+1,J)+T0(I,J-1)+T0(I,J+1))
C
END DO
END DO
C
DO J = J0, J1
T(I1,J) = T0(I1-1,J)
END DO

999 RETURN
END
C
C
C
C *********************************************************************
SUBROUTINE BOU_T(I0,I1,J0,J1,T,T1,T2,T3)
INCLUDE ‘implicit.for’
C———————————————— Global Variable —–
REAL*4 T(ISSS:IEEE,JSSS:JEEE)          !— Thermal
REAL*4 T1                              !— Upeer Temperature
REAL*4 T2                              !— Left Temperature
REAL*4 T3                            !— Bottom Temperature
INTEGER I0, I1, J0, J1           !— Mesh Size
C————————————————- Local Variable —–
INTEGER I,J
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

C
DO J=J0,J1
DO I=I0,I1
T(I,J0) = T3
T(I0,J) = T2
T(I,J1)= T1
END DO
END DO
C
999 RETURN
END
C
C
C
C *********************************************************************
SUBROUTINE CAL_R_T (I0,I1,J0,J1,T,T0,RES_T)
INCLUDE ‘implicit.for’
C———————————————— Global Variable —–
REAL*4 T(ISSS:IEEE,JSSS:JEEE)          !— Thermal
REAL*4 T0 (ISSS:IEEE,JSSS:JEEE) !— Temperature for Meomory
INTEGER I0, I1, J0, J1           !— Mesh Size
REAL*4 RES_T                     !— Residual T
C————————————————- Local Variable —–
INTEGER I , J
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
RES_T=-1.0E+10

DO J=J0+1,J1-1
DO I=I0+1,I1-1
RES_T = AMAX1( RES_T , ABS(T(I,J)-T0(I,J)) )
END DO
END DO
C
999 RETURN
END
C
CC**********************************************************************
SUBROUTINE OUTPUT_T1 (I0,I1,J0,J1,X_T,Y_T,T)
INCLUDE ‘implicit.for’
C———————————————— Global Variable —–
REAL*4 X_T(ISSS:IEEE,JSSS:JEEE)        !— Position X-component
REAL*4 Y_T(ISSS:IEEE,JSSS:JEEE)        !— Position Y-component
REAL*4 T(ISSS:IEEE,JSSS:JEEE)          !— Thermal
INTEGER I0, I1, J0, J1           !— Mesh Size
C————————————————- Local Variable —–
INTEGER I,J
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
112 FORMAT(2(E16.8E3,1X))
113 FORMAT(3(E16.8E3,1X))
114 FORMAT(4(E16.8E3,1X))
C
OPEN (11,FILE=’./data/temp.dat’,STATUS=’unknown’)
DO J=J0,J1
DO I=I0,I1
WRITE(11,113) X_T(I,J), Y_T(I,J) , T(I,J)
END DO
END DO
CLOSE(11)
C
RETURN
END
C
C