隕石衝突後の熱伝導シュミレーション TDC 2006/09/02 23:04:00 
　　└<続き> TDC 2006/09/02 23:06:07 
　　　├単位が合っていない 河川屋 2006/09/18 09:47:59 
　　　└実行時間を大幅に短縮することができますよ... 荒田浩二 2006/09/21 20:18:10 
　　　　├＜続き＞ 荒田浩二 2006/09/21 20:20:06 
　　　　└＜続き＞ 荒田浩二 2006/09/21 20:21:29 

　　隕石衝突後の熱伝導シュミレーション TDC 2006/09/02 23:04:00  ツリーへ

隕石衝突後の熱伝導シュミレーション	返事を書く
TDC 2006/09/02 23:04:00
ＮＨＫの地球大進化という番組で 全海洋蒸発をした際の熱シュミレーションをしていました。 簡易版ができたので発表します。 前半 ! ■■■■巨大隕石衝突後の熱伝導シュミレーション■■■■ !----------初期設定0------------- OPTION ARITHMETIC NATIVE OPTION BASE 0 LET WID=200 LET depth=500 DIM ATM(WID,depth),T(WID,depth) LET K273=273 LET B=10 !標準文字サイズ SET TEXT FONT "ＭＳ ゴシック",B SET WINDOW -1,WID,-depth,depth SET TEXT COLOR 4 ! 白0黒1青2緑3赤4水5黄6 !-------------------------------- !----------初期設定1------------- LET d=1.5 !土の熱伝導率 Ｊ/(sec・ｍ・Ｋ） 温度差が１℃の時、１ｍの厚さを１秒に流れる熱量 !　コンクリートの体積比熱(≒1.92倍)を参考にした仮の値 LET J_BOX=1.9*(4.18*10^6) ! 土の熱容量 Ｊ/(ｍ3・Ｋ) LET σ=5.67*10^(-8) ! (watt/�udeg4) ! 放射の速度は、シュテファン・ボルツマンの法則 ! Ｍ＝σ・Ｔ^4(watt/�u) ! σ:5.67×10^-8(watt/�udeg4) T：絶対温度 !-------------------------------- CALL ini ! 地中の温度を40℃に設定 PLOT TEXT ,AT 0,depth*0.9 : "■熱伝導シュミレーション■ 〜巨大隕石衝突後〜　しばらくおまちください" PLOT TEXT ,AT 0,depth*0.2 : " 巨大隕石が衝突！！" PLOT TEXT ,AT 0,depth*0.05 : "地表面" PRINT "隕石が衝突！！" CALL ini4000 PRINT "地表を4000℃に設定で描画" CALL DRAWING PRINT "本計算開始！！" DO IF TIM<24*365.25 THEN CALL ini4000 ELSE LET T(0,0)=T(0,0)-σ*T(0,0)^4/J_BOX ! 熱放射式 シュテファン・ボルツマン END IF IF INT(TIM/1000)=TIM/1000 THEN PRINT TIM;"時間経過" END IF CALL CALC1 !計算部モジュール CALL TEXT1 !文字表示モジュール LET TIM=TIM+1 LOOP UNTIL T(0,0)<=50+K273 AND TIM>0 PRINT "終了！！"

　　└<続き> TDC 2006/09/02 23:06:07  ツリーへ

Re: 隕石衝突後の熱伝導シュミレーション	返事を書く
TDC 2006/09/02 23:06:07
<続き> !-------------------サブルーチン------------------------- SUB TEXT1 IF TIM<>0 AND INT(TIM/1000)=TIM/1000 THEN CLEAR PLOT TEXT ,AT 0,depth*0.9 : "■熱伝導シュミレーション■ 〜巨大隕石衝突後〜" PLOT TEXT ,AT 0,depth*0.8 : " 隕石衝突後 "&STR$(TIM)&"時間(＝ "&STR$(ROUND(TIM/24/365.25,2))&"年)経過" PLOT TEXT ,AT 0,depth*0.57 : " 土の熱伝導率(仮)： "&STR$(d)&"Ｗ/(ｍ・deg)" PLOT TEXT ,AT 0,depth*0.5 : " 土の熱容量(仮)： "&STR$(J_BOX/1000)&"kＪ/(ｍ3・deg)" IF TIM>24*365.25 THEN SET TEXT COLOR 2 PLOT TEXT ,AT 0,depth*0.2 : "◆地球が冷える・・・熱放射式 Ｗ ＝ σ・Ｔ^4" SET TEXT COLOR 4 ELSEIF TIM<24*365.25 THEN PLOT TEXT ,AT 0,depth*0.2 : " 地球は岩石蒸気で覆われる・・・ 約１年間" END IF IF T(0,0)<100+K273 THEN SET TEXT COLOR 2 PLOT TEXT ,AT 0,depth*0.12 : " そして、雨も降りはじめた・・・ " SET TEXT COLOR 4 END IF PLOT TEXT ,AT 0,depth*0.05 : "地表面 温度 "&STR$(ROUND(T(0,1)-K273,2))&"℃" CALL DRAWING SET TEXT COLOR 4 PLOT TEXT ,AT WID*0.8,depth*0.2 :"■〜400℃" SET TEXT COLOR 6 PLOT TEXT ,AT WID*0.8,depth*0.15 :"■400〜100℃" SET TEXT COLOR 3 PLOT TEXT ,AT WID*0.8,depth*0.1 :"■100〜50℃" SET TEXT COLOR 2 PLOT TEXT ,AT WID*0.8,depth*0.05 :"■50℃以下" SET TEXT COLOR 5 PLOT TEXT ,AT WID*0.45,-depth*0.0 : "深さ："&STR$(depth*0)&"m" PLOT TEXT ,AT WID*0.45,-depth*0.5 : "深さ："&STR$(depth*0.5)&"m" PLOT TEXT ,AT WID*0.45,-depth : "深さ："&STR$(depth)&"m" SET TEXT COLOR 4 ! 白0黒1青2緑3赤4水5黄6 END IF END SUB SUB ini4000 !地表1ｍを4000℃に設定 FOR I=0 TO 1 FOR J=0 TO 1 LET T(I,J)=4000+K273 NEXT J NEXT I END SUB SUB CALC1 FOR J=1 TO depth-1 LET �冲emp=T(0,J-1)-T(0,J) !　境界面の温度差を計算 LET �價=(d*1)*�冲emp*60*60 ! 1m 1時間に移動するエネルギー LET T(0,J)=�價/J_BOX+T(0,J) LET T(0,J-1)=-�價/J_BOX+T(0,J-1) NEXT J END SUB SUB DRAWING ! ///〜400℃赤/400〜100℃黄/100〜50℃緑/50℃以下青/// FOR J=0 TO depth IF T(0,J)<=50+K273 THEN SET LINE COLOR 2 ! 白0黒1青2緑3赤4水5黄6 IF T(0,J)>50+K273 AND T(0,J)<=100+K273 THEN SET LINE COLOR 3 IF T(0,J)>100+K273 AND T(0,J)<=400+K273 THEN SET LINE COLOR 6 IF T(0,J)>400+K273 THEN SET LINE COLOR 4 SET LINE STYLE 1 ! 1実線 2破線 3点線 4一点鎖線 SET LINE WIDTH 1 PLOT LINES:0,-J;WID,-J NEXT J END SUB SUB ini FOR J=0 TO depth LET T(0,J)=40+K273 NEXT J END SUB END

　　　├単位が合っていない 河川屋 2006/09/18 09:47:59  ツリーへ

Re: <続き>	返事を書く
河川屋 2006/09/18 09:47:59
単位が合っていない この部分ですが、 LET T(0,0)=T(0,0)-σ*T(0,0)^4/J_BOX ! 熱放射式 シュテファン・ボルツマン LET T(0,0)=T(0,0)-σ*60*60*T(0,0)^4/J_BOX というふうに、時間(1時間=3600秒)を掛けてやらないと単位が合わないのですが。 それと、この式は、「厚さ１ｍの岩から黒体輻射で放出される分、温度が下がる」 とやっていますが、何か忘れていません？ 初期温度４０００度だから、海水は全部蒸発し大気となっていて、 大気(＝海水)の持っている熱容量は、厚さ１ｍの岩より桁違いに大きいです。 だから、 LET W=2E6 ! 大気(水蒸気)の質量(kg/m2) LET JA=1.5E3 !大気(水蒸気)の比熱(J/kgK) LET T(0,0)=T(0,0)-σ*60*60*T(0,0)^4/W/JA !放射による冷却 あと、地下に伝導する熱より、 ・水蒸気による温室効果 ・太陽からの熱の受領 ・黒体輻射の基準温度を、地表にするか雲頂にするか ・水蒸気→水の潜熱 のほうが効きそう。 よって、簡略版なら、地下への熱伝導は無視し、大気の熱収支だけ考えればよいと思います。

　　　└実行時間を大幅に短縮することができますよ... 荒田浩二 2006/09/21 20:18:10  ツリーへ

Re: <続き>	返事を書く
荒田浩二 2006/09/21 20:18:10
実行時間を大幅に短縮することができますよ。 河川屋さんが専門的なアドバイスをされてますが、 私はそっち方面の知識はないので純粋に実行時間の短縮を図りました。 改良点は次の４点です。 １。print文での大量の出力を回避する print文は出力結果が累積されるにつれ実行時間が遅くなります。 おそらく行送りに時間がかかるのでしょうが、数千行を超えると目に見えて遅くなります。 このプログラムの実行時間を遅くしている最大の要因です。 IF INT(TIM/1000)=TIM/1000 THEN PRINT TIM;"時間経過" END IF この３行を削除するだけでかなりの時間を短縮できます。 ２。ループを分割する 条件文(if文)は単純実行文より実行に時間がかかります。 数万回繰り返されるループの中に、条件に合致すると判定されることが１回しかないif文が何箇所かにあるのならば、ループを分割してしまったほうが良いでしょう。 ソースの記述は冗長になりますが時短を優先させました。 ３。ループ内での演算を少なくする できるだけループ内での演算回数を少なくするために、同じ演算はあらかじめ計算しておき定数化しました。 内部副プログラムの'CALC1'と'DRAWING'のループの中で、いかに演算を少なくするかが鍵です。 ４。グラフィック画面では変更のある箇所だけ描画するようにする これは、以前この掲示板で哲さんに教えていただいたテクニックです。（掲示板No.78） clear文はどうしても画面がちらつきます。 改良したプログラムでは私のＰＣで、雨が降りはじめるまでに１時間４分、終了までに５時間１８分かかりました。 なお、新しく設定した変数は次のとおりです。 interval,year,dep08,dep005,tex1,tex2,WID045,mdep05,mdepth,tex3$,tex4$,dJB,σJB,�價B,C ! ■■■■巨大隕石衝突後の熱伝導シュミレーション■■■■ !----------初期設定0------------- OPTION ARITHMETIC NATIVE OPTION BASE 0 LET WID=200 LET depth=500 DIM T(depth) !! 1次の配列にしました。 ATM(WID,depth)？ LET K273=273 LET B=10 !標準文字サイズ SET TEXT FONT "ＭＳ ゴシック",B SET WINDOW -1,WID,-depth,depth SET TEXT COLOR 4 ! 白0黒1青2緑3赤4水5黄6 SET LINE STYLE 1 ! 1実線 2破線 3点線 4一点鎖線 !! SUB DROWING より移動 SET LINE WIDTH 1 !! SUB DROWING より移動 LET interval=1000 !! 描画間隔 1000時間 LET year=1/24/365.25 !! TIM を 年 に換算 LET dep08=depth*0.8 !! 定数化 LET dep005=depth*0.05 !! 定数化 LET tex1$="" !! LET tex2$="" !! LET WID045=WID*0.45 !! 定数化 LET mdep05=-depth*0.5 !! 定数化 LET mdepth=-depth !! 定数化 LET tex3$="深さ："&STR$(depth*0.5)&"m" !! 定数化 LET tex4$="深さ："&STR$(depth)&"m" !! 定数化 !-------------------------------- !----------初期設定1------------- LET d=1.5 !土の熱伝導率 Ｊ/(sec・ｍ・Ｋ） 温度差が１℃の時、１ｍの厚さを１秒に流れる熱量 !　コンクリートの体積比熱(≒1.92倍)を参考にした仮の値 LET J_BOX=1.9*(4.18*10^6) ! 土の熱容量 Ｊ/(ｍ3・Ｋ) LET σ=5.67*10^(-8) ! (watt/�udeg4) LET dJB=(d*1)*60*60/J_BOX !! 定数化 LET σJB=σ/J_BOX !! 定数化 ! 放射の速度は、シュテファン・ボルツマンの法則 ! Ｍ＝σ・Ｔ^4(watt/�u) ! σ:5.67×10^-8(watt/�udeg4) T：絶対温度 !-------------------------------- ＜続く＞

　　　　├＜続き＞ 荒田浩二 2006/09/21 20:20:06  ツリーへ

Re: 実行時間を大幅に短縮することができますよ...	返事を書く
荒田浩二 2006/09/21 20:20:06
＜続き＞ CALL ini ! 地中の温度を40℃に設定 PLOT TEXT ,AT 0,depth*0.9 : "■熱伝導シュミレーション■ 〜巨大隕石衝突後〜　しばらくおまちください" PLOT TEXT ,AT 0,depth*0.2 : " 巨大隕石が衝突！！" PLOT TEXT ,AT 0,depth*0.05 : "地表面" PRINT "隕石が衝突！！" CALL ini4000 PRINT "地表を4000℃に設定で描画" CALL DRAWING CLEAR !! SUB TEXT1より移動 SET TEXT COLOR 4 PLOT TEXT ,AT 0,depth*0.9 : "■熱伝導シュミレーション■ 〜巨大隕石衝突後〜" PLOT TEXT ,AT 0,depth*0.57 : " 土の熱伝導率(仮)： "&STR$(d)&"Ｗ/(ｍ・deg)" PLOT TEXT ,AT 0,depth*0.5 : " 土の熱容量(仮)： "&STR$(J_BOX/1000)&"kＪ/(ｍ3・deg)" PLOT TEXT ,AT 0,depth*0.2 : " 地球は岩石蒸気で覆われる・・・ 約１年間" SET TEXT COLOR 4 !! SUB TEXT1より移動 PLOT TEXT ,AT WID*0.8,depth*0.2 :"■〜400℃" SET TEXT COLOR 6 PLOT TEXT ,AT WID*0.8,depth*0.15 :"■400〜100℃" SET TEXT COLOR 3 PLOT TEXT ,AT WID*0.8,depth*0.1 :"■100〜50℃" SET TEXT COLOR 2 PLOT TEXT ,AT WID*0.8,depth*0.05 :"■50℃以下" SET TEXT COLOR 5 !! SUB TEXT1より移動 PLOT TEXT ,AT WID045,-depth*0.0 : "深さ："&STR$(depth*0)&"m" PRINT "本計算開始！！" DO WHILE TIM<24*365.25 !! DO文を分割 #1 CALL ini4000 FOR J=1 TO depth-1 LET �價B=(T(J-1)-T(J))*dJB !! LET T(J)=�價B+T(J) LET T(J-1)=-�價B+T(J-1) NEXT J IF TIM<>0 AND INT(TIM/interval)=TIM/interval THEN CALL TEXT1 !文字表示モジュール LET TIM=TIM+1 LOOP SET TEXT COLOR 0 !! SUB TEXT1より移動 PLOT TEXT ,AT 0,depth*0.2 : " 地球は岩石蒸気で覆われる・・・ 約１年間" !! 消去 SET TEXT COLOR 2 PLOT TEXT ,AT 0,depth*0.2 : "◆地球が冷える・・・熱放射式 Ｗ ＝ σ・Ｔ^4" DO !! DO文を分割 #2(TIMを1000時間単位にそろえるための調整) LET T(0)=T(0)-T(0)^4*σJB !! 熱放射式 シュテファン・ボルツマン FOR J=1 TO depth-1 LET �價B=(T(J-1)-T(J))*dJB !! LET T(J)=�價B+T(J) LET T(J-1)=-�價B+T(J-1) NEXT J IF INT(TIM/interval)=TIM/interval THEN EXIT DO LET TIM=TIM+1 LOOP CALL TEXT1 !文字表示モジュール DO !! DO文を分割 #3 CALL CALC1 !計算部モジュール CALL TEXT1 !文字表示モジュール LOOP UNTIL T(0)<100+K273 SET TEXT COLOR 2 !! SUB TEXT1より移動 PLOT TEXT ,AT 0,depth*0.12 : " そして、雨も降りはじめた・・・ " DO !! DO文を分割 #4 CALL CALC1 !計算部モジュール CALL TEXT1 !文字表示モジュール LOOP UNTIL T(0)<=50+K273 PRINT "終了！！" ＜続く＞

　　　　└＜続き＞ 荒田浩二 2006/09/21 20:21:29  ツリーへ

Re: 実行時間を大幅に短縮することができますよ...	返事を書く
荒田浩二 2006/09/21 20:21:29
＜続き＞ !-------------------サブルーチン------------------------- SUB TEXT1 !! 変更のあるテキストのみ印字 SET TEXT COLOR 0 PLOT TEXT ,AT 0,dep08 : tex1$ !! 文字消去 LET tex1$=" 隕石衝突後 "&STR$(TIM)&"時間(＝ "&STR$(ROUND(TIM*year,2))&"年)経過" SET TEXT COLOR 4 PLOT TEXT ,AT 0,dep08 : tex1$ SET TEXT COLOR 0 PLOT TEXT ,AT 0,dep005 : tex2$ !! 文字消去 LET tex2$="地表面 温度 "&STR$(ROUND(T(1)-K273,2))&"℃" SET TEXT COLOR 4 PLOT TEXT ,AT 0,dep005 : tex2$ CALL DRAWING END SUB SUB ini4000 !地表1ｍを4000℃に設定 FOR J=0 TO 1 LET T(J)=4000+K273 NEXT J END SUB SUB CALC1 !! FOR I=1 TO interval !! 描画間隔1000時間の温度計算 LET T(0)=T(0)-T(0)^4*σJB !! 熱放射式 シュテファン・ボルツマン FOR J=1 TO depth-1 LET �價B=(T(J-1)-T(J))*dJB !! LET T(J)=�價B+T(J) LET T(J-1)=-�價B+T(J-1) NEXT J LET TIM=TIM+1 !! NEXT I END SUB SUB DRAWING !! ! ///〜400℃赤/400〜100℃黄/100〜50℃緑/50℃以下青/// FOR J=0 TO depth LET C=T(J)-K273 !! SELECT CASE文と同内容 IF C<=50 THEN SET LINE COLOR 2 ! 白0黒1青2緑3赤4水5黄6 ELSEIF C<=100 THEN SET LINE COLOR 3 ELSEIF C<=400 THEN SET LINE COLOR 6 ELSE SET LINE COLOR 4 END IF PLOT LINES:0,-J;WID,-J NEXT J SET TEXT COLOR 5 !! PLOT TEXT ,AT WID045,mdep05 : tex3$ PLOT TEXT ,AT WID045,mdepth : tex4$ END SUB SUB ini FOR J=0 TO depth LET T(J)=40+K273 NEXT J END SUB END