原子動力工学特論

1. 原子動力工学特論 課題　3、4 交通電子機械工学専攻 2003310　　齋藤　泰治

2. 発表内容　1 ・課題.3　（フィン内部の温度分布） ・ガウスの消去法においての留意点 　　　：軸要素（pivot）、正規化（scaling）について ・プログラム ・実行結果

3. 発表内容　2 ・課題.4　（平板の非定常熱伝導解析） ・差分法（FDE)による解法 ・プログラム ・実行結果 ・フーリエ変換（PDE)による解法 ・プログラム ・実行結果 ・FDEとPDEの比較

4. 課題.3 Qｃ T∞ 400℃ Qｎ Qｎ+1 300℃ Δx Qｎ＝Qｎ+1+Qｃ

6. 軸要素（pivot)及び正規化（scaling) 6本の連立方程式として解けばよい。 課題2の問2において3本の連立方程式を解くプログラム（3×3）を作成したので、それを6×6に拡張する。 しかし、課題2でのプログラムは、単純にガウスの消去法を行わしただけであったので、そのままだと不具合を起こす。 というのも、人間が手計算で行う時には起こらない、コンピューターでの数値計算において、留意しなければならない問題が発生するからである。

7. 軸要素（pivot） 各段階においけるannで除算を行ってゆく。 これを軸要素（pivot）といい、 0に近いと次の計算での桁落ちが以降すべてに影響し、計算が不安定になる。

8. 正規化（scaling） 方程式の係数の絶対値が極端に違っていると、計算の過程で桁落ち誤差が発生し易くなる。また、先に述べた軸要素（pivot）の問題も発生しやすくなる。そのため、計算前に各係数の絶対値が同程度になるように、調整を行う。これを正規化（scaling）という。 Siは係数aijの最大値、 ｔiはaij/Siの最大値として、 全ての係数の絶対値をそろえる。 yiを未知数として解けば、解はxiとなる

9. //----------Gauss-Jordan---------- void sweep(double *pa,double *pt,int n,int m,int *piwork,int ILL) { double max,w; int i,j,k,iw,p,q; for(i=0;i<n;++i){ //正規化 max=fabs(*(pa+m*i)); for(j=0;j<n;++j){ if(max<fabs(*(pa+m*i+j))){ max=fabs(*(pa+m*i+j)); } } for(j=0;j<m;++j){ *(pa+m*i+j)=*(pa+m*i+j)/max; } } for(j=0;j<n;++j){ max=fabs(*(pa+j)); for(i=0;i<n;++i){ if(max<fabs(*(pa+m*i+j))){ max=fabs(*(pa+m*i+j)); } } *(pt+j)=max; for(i=0;i<n;++i){ *(pa+m*i+j)=*(pa+m*i+j)/max; } }

10. for(i=0;i<n;++i){ //列の入れ替えをする配列の初期設定for(i=0;i<n;++i){ //列の入れ替えをする配列の初期設定 *(piwork+i)=i; } for(k=0;k<n;++k){ //掃き出し法 max=fabs(*(pa+m*k+k)); //絶対値の最大の要素を求める（pivot選択） p=k; q=k; for(j=k;j<n;++j){ for(i=k;i<n;++i){ if(max=fabs(*(pa+m*i+j))){ max=fabs(*(pa+m*i+j)); p=i; q=j; } } } if(max<=eps){ //最大値が極端に小さい場合は ILL=1; //解不能として終了 printf("Matrix is ILL\n"); return; } for(i=0;i<n;++i){ //k列とq列を入れ替える w=*(pa+m*i+k); *(pa+m*i+k)=*(pa+m*i+q); *(pa+m*i+q)=w; } for(j=k;j<m;++j){ //k行とp行の入れ替え w=*(pa+m*k+j); *(pa+m*k+j)=*(pa+m*p+j); *(pa+m*p+j)=w; }

11. i=*(piwork+k); //入れ替えた列を記憶 *(piwork+k)=*(piwork+q); *(piwork+q)=i; for(j=k+1;j<m;++j){ *(pa+m*k+j)=*(pa+m*k+j)/(*(pa+m*k+k)); //割る } for(i=0;i<n;++i){ //引く if(i!=k){ for(j=k+1;j<m;++j){ *(pa+m*i+j)=*(pa+m*i+j)-*(pa+m*i+k)*(*(pa+m*k+j)); } } } } for(j=n;j<m;++j){ //入れ替えた列を戻す for(i=0;i<n;++i){ iw=*(piwork+i); *(pa+m*iw+n-1)=*(pa+m*i+j); } for(i=0;i<n;++i){ *(pa+m*i+j)=*(pa+m*i+n-1)/(*(pt+i)); //正規化されているので } //作業領域から戻し } //未知数を決める return; }

12. 実行結果 kax=0.235619 hpx=0.084823 X[11]= 1.0 X[12]= 0.5 X[13]= 0.0 X[14]= 0.0 X[15]= 0.0 X[16]= 0.0 X(17)=321.4 X[21]= 0.0 X[22]=-0.8 X[23]= 0.2 X[24]= 0.0 X[25]= 0.0 X[26]= 0.0 X(27)=-365.7 X[31]= 0.0 X[32]= 0.2 X[33]=-0.6 X[34]= 0.2 X[35]= 0.0 X[36]= 0.0 X(37)=-48.6 X[41]= 0.0 X[42]= 0.0 X[43]= 0.2 X[44]=-0.6 X[45]= 0.2 X[46]= 0.0 X(47)=-48.6 X[51]= 0.0 X[52]= 0.0 X[53]= 0.0 X[54]= 0.2 X[55]=-0.6 X[56]= 0.2 X(57)=-48.6 X[61]= 0.0 X[62]= 0.0 X[63]= 0.0 X[64]= 0.0 X[65]= 0.2 X[66]=-3.1 X(67)=-1644.4 Answer Q(1)=1.000000 T(2)=644.0113 T(3)=611.5981 T(4)=593.0802 T(5)=581.7912 T(6)=573.6670

14. 課題.4　平板の非定常熱伝導方程式 a=2×10-4 (m2/s) L=100 (mm) Initial condition: T(x)=100[℃] at t=0.0 Boundary condition T(x)=0[℃] at x=(0,L) 100℃ L 0℃

15. n+1 n： 時間メッシュ n m-1 m m+1 m: 空間メッシュ 差分法（FDE）による解法 熱伝導方程式 これを差分表示し、整理する unkown known

16. #define nn 10000 #define mm 100  double T[nn][mm];  main() { double a,dt,dx; int m,n; a=2.0*(pow(10,-4)); dt=2.0*(pow(10,-3)); dx=1.0*(pow(10,-3)); for(n=1;n<=10000;++n){ //　境界条件 T[n][1]=T[n][100]=0.0; } for(n=1;n<=10000;++n){ //　初期温度（100℃一様） for(m=2;m<=99;++m){ T[n][m]=100.0; } } for(n=1;n<=9999;++n){//　差分計算 for(m=2;m<=99;++m){ T[n+1][m]=T[n][m]*(1-2*a*dt/(pow(dx,2)))+a*dt/(pow(dx,2))*(T[n][m+1]+T[n][m-1]); printf("%lf ",T[n][m]); } printf("\n"); } }

18. 熱伝導方程式（拡散方程式） フーリエ変換（PDE）による解法 これをフーリエ変換して解くと、 従って、これが解析解となる。

19. #include <stdio.h> #include <math.h> #define pi 3.1415926535 main() { double T; double T0,L,a; int t,x,n; T0=100.0; L=100*pow(10,-3); a=2*pow(10,-4); for(t=0;t<=100;++t){ printf("t=%d\n",t); for(x=0;x<=100;++x){ T=0.0; for(n=1;n<=100;n+=2){ T=2/(n*pi)*T0*2*sin(n*pi/L*x*(pow(10,-3)))*exp(-a*pow(n,2)*pow((pi/L),2)*t*0.01)+T; } printf("%lf ",T); } printf("\n"); }  }

20. ×10-3[s]

21. FDE PDE