Presentation is loading. Please wait.

Presentation is loading. Please wait.

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

Similar presentations


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 Qc T∞ 400℃ Qn Qn+1 300℃ Δx Qn=Qn+1+Qc      

5

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

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

8 正規化(scaling) 方程式の係数の絶対値が極端に違っていると、計算の過程で桁落ち誤差が発生し易くなる。また、先に述べた軸要素(pivot)の問題も発生しやすくなる。そのため、計算前に各係数の絶対値が同程度になるように、調整を行う。これを正規化(scaling)という。 Siは係数aijの最大値、 tiは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


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

Similar presentations


Ads by Google