朝眠いながらもなんとか起きて学校へ。ツナギを置いておく場所がほしいなあ。
1・2限は加工実習。今日も3時間立っているだけなので暇でしょうがない。旋盤で軸の大きいほうを加工しているのだけれど、旋盤を動かすのは一人だし、径とか記録するのも一人でいいし、やったことと言えばマイクロメータでたまに径とかを測っただけ。なのでひたすらストレッチとかしていた。暇だ。
昼はおにぎりを食べて、プログラムを印刷してきた。まだ途中なので授業中に紙の上でプログラムを書こうと思って。最近学食を全く使わない。それはお金がないから。
3限は燃焼工学。いつも暇でしょうがない授業なので、プログラムでも書いていようと思っていた。内容はNS方程式とかエネルギ式とかをいじっていた気がする。後はエンジンの乱流がどうこうって言ってたなあ。途中で当てられたとき、「○○研の彼」（○○はタカフの研究室）とか言われて俺じゃないなと思ってたら俺だったらしい。どうして先生はそんな間違いをしたのやら。
空き時間にもプログラムをする。半径方向熱伝導なので一次精度風上差分+拡散方程式の差分なのだけれど、ここでのCFL条件は拡散の式が利いてくる気がするので、時間ステップと空間ステップと熱伝導率に注意しなければならないなと思う。
5限は航空構造力学。今日はなんかこの前とかと雰囲気が変わって、曲面の方程式をやるとか言っていて、話が数学的になっていた。曲率半径とかそんなこと言っていたなと思うけれど、最後のほうがちょっと聞いていなかったら分からなくなってので諦めた。どこかでfollowしないとなあ。
家に帰り、ご飯を食べて少ししたら寝てしまった。そして起きてまたプログラムを始める。考えつつ打っていたらプログラム自体は完成して、実行できるようになったのだけれど、実行してみるとめっちゃ時間がかかる。目安がないと計算時間が1分なのか1時間なのか1日なのか分からないので％表示をしてみた。それからコードもさらしておこう。CFL条件が満たされているのか激しく不安。
と思っていたけれどやっぱこんなに遅いのはプログラムのミスっぽい。またデバッグの繰り返しかあ。色々おかしいところがあって今も直してて大変だわ。
実行画面。パソコンの性能が如実に速度に出るんだと思う。CeleronD340(2.93GHz)の限界。
c ノズル熱伝達+熱伝導解析： c 暗黙の型宣言を使わない宣言 c implicit none c 刻み時間、刻みセル設定 real*8 dt,dr parameter(dt=0.01,dr=0.0001) c kiza_numはy方向刻み数、kiza_xはx方向刻み数 integer kiza_num,kiza_x,kiza_t parameter(kiza_num=150,kiza_x=2370,kiza_t=6000) c kiza_num=15.0/dr c kiza_t=60.0/dt c 繰り返しの数の宣言 integer i,j,k real*8 l c temp(x,y)はx,y位置の温度、一次温度の宣言 real*8 temp(1:kiza_x,1:kiza_num) real*8 temp_ex(kiza_num) real*8 temp_temp(3) c 熱伝達率hgの設定： real*8 hg_int(kiza_x) c Baltzの式の配列: real*8 a(kiza_x) real*8 c(kiza_x) real*8 r(kiza_x) c ガス温度の配列: real*8 temp_gas(kiza_x) c ノズル内側半径の配列: real*8 radius(kiza_x) c q()は熱流束の配列 real*8 q(kiza_x) c 熱伝導計算の半径r real*8 rad_in(1:kiza_num) c hgの定数部分 real*8 :: t0=1773.15 real*8 :: u=0.0000181 real*8 :: cp=1.007 real*8 :: pr=0.70 c c-c複合材の初期温度 real*8 :: T_shoki=20.0 c c-c複合材の比熱と密度と円周率設定 real*8 :: hinetsu=1800.0 real*8 :: mitsudo=1800.0 real*8 :: pi=3.14159 c 熱伝導の係数関数 real*8 :: k_heat01=50.0 real*8 k_heat02,k1 c Baltzの計算部分 c hg4のxはマッハ数 hg4(x)=((1.0+0.4/1.4*(x**2.0))**(-0.12)) c hg3のyが壁面温度、xがマッハ数 hg3(x,y)=((0.5*y/t0*(1.0+0.4/1.4*(x**2.0))+0.5)**(-0.68)) c hg2のxが断面積 hg2(x)=((0.428319/x)**0.8) c hg1のxは等価直径 hg1(x)=(0.026/x**0.2)*((u**0.2)*cp/(pr**0.6)) c 熱伝導率の係数k2 k2(x)=k_heat02*dt/x/dr c 指標1 write(*,*)'壱：ここまで終わったよ！' c 断面積、マッハ数、等価直径、ガス静温、内側半径の読み込み open(10,file='areaonly.d',status='old') open(11,file='machonly.d',status='old') open(12,file='r.d',status='old') open(13,file='temponly.d',status='old') open(14,file='heat_radius.d',status='old') do i=1,kiza_x read(10,*)a(i) read(11,*)c(i) read(12,*)r(i) read(13,*)temp_gas(i) read(14,*)radius(i) r(i)=r(i)*0.001 a(i)=a(i)*0.000001 radius(i)=radius(i)*0.001 end do close(10) close(11) close(12) close(13) close(14) c tempの初期化 do i=1,kiza_x do j=1,kiza_num temp(i,j)=273.15+T_shoki end do end do c 指標2 write(*,*)'弐：ここまで終わったよ！' c 温度伝導率の係数の計算 k_heat02=k_heat01/(mitsudo*hinetsu) k1=k_heat02*dt/(dr**2.0) c 時間刻み。60秒/刻み時間 do i=1,kiza_t c 指標3 c do k=10,kiza_t,100 c if(i==k)then l=i/(kiza_t*1.0)*100.00 write(*,*)i,'/',kiza_t,'回 ',l,'％終了' c end if c end do c x方向刻み do j=1,kiza_x c 熱伝達率の計算 hg_int(j)=hg1(r(i))*hg2(a(i))*hg3(temp(j,1),c(i))*hg4(c(i)) c 熱流束qの計算 q(j)=hg_int(j)*(temp_gas(j)-temp(j,1)) c 壁部分の初期温度計算 temp(j,1)=temp(j,1)+q(j)/(hinetsu*mitsudo*2.0*pi*radius(j))*dt c 前ステップの温度をtemp_exに代入 do k=1,kiza_num temp_ex(k)=temp(j,k) end do c 熱伝導計算 do k=2,kiza_num rad_in(k)=radius(j)-(k-1)*dr temp_temp(1)=k1*temp_ex(k+1) temp_temp(2)=(1+k2(rad_in(k))-2.0*k1)*temp_ex(k) temp_temp(3)=(k1-k2(rad_in(k)))*temp_ex(k-1) temp(j,k)=temp_temp(1)+temp_temp(2)+temp_temp(3) c 最後だけ条件分岐 if(k==kiza_num)then temp_temp(1)=k1*temp_ex(k) temp_temp(2)=(1+k2(rad_in(k))-2.0*k1)*temp_ex(k) temp_temp(3)=(k1-k2(rad_in(k)))*temp_ex(k-1) temp(j,k)=temp_temp(1)+temp_temp(2)+temp_temp(3) end if end do end do end do open(20,file='heat_results.d',status='unknown') do j=1,kiza_x write(20,*)temp(j,1),j end do close(20) end
コメント