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

Presentation on theme: "原子動力工学特論 課題　3、4 交通電子機械工学専攻 2003310　　齋藤　泰治."— Presentation transcript:

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

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

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

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

5

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; max=fabs(*(pa+j)); for(i=0;i<n;++i){ *(pt+j)=max;

10 for(k=0;k<n;++k){ //掃き出し法
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){ *(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); *(pa+m*i+j)=*(pa+m*i+n-1)/(*(pt+i)); //正規化されているので } //作業領域から戻し } //未知数を決める return;

12 実行結果 kax= hpx= 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)= Answer Q(1)= T(2)= T(3)= T(4)= T(5)= T(6)=

13

14 課題.4 平板の非定常熱伝導方程式 a=2×10-4 (m2/s) L=100 (mm) 100℃ Initial condition:
課題.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 差分法（FDE）による解法 熱伝導方程式 n： 時間メッシュ m: 空間メッシュ m-1 m m+1 n n+1 これを差分表示し、整理する
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){　　　　　　　//　差分計算 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");

17

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

19 #include <stdio.h>
#include <math.h> #define pi 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