巨大隕石衝突後の熱シミュレーション （条件／隕石衝突後１年間が、４０００度の場合）
(隕石衝突の後、10年後の結果)
※温室効果ガス（水蒸気など）が考慮されていないのでもっと厳しい結果になる。
※衝突点の反対側の極（北極又は南極）地方は、直射の太陽光の熱供給は考慮し
なくてよい分、このモデルに近くなる。
! 十進ＢＡＳＩＣ
! ■■■■巨大隕石衝突後の熱伝導シミュレーション■■■■
!----------初期設定0-------------
OPTION ARITHMETIC NATIVE
OPTION BASE 0
LET WID=200
LET depth=600
DIM T(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
LET MAX_TE=4000+K273
!--------------------------------
!----------初期設定1-------------
LET d=1.5 !土の熱伝導率 Ｊ/(sec・ｍ・Ｋ） 温度差が１℃の時、１ｍの厚さを１秒に流れる熱量
! コンクリートの体積比熱(≒1.92倍)を参考にした仮の値
LET J_BOX=1.9*(4.18*10^6) ! 土の熱容量 Ｊ/(ｍ3・Ｋ)
LET σ=5.67*10^(-8) ! (watt/㎡deg4)
! 放射の速度は、シュテファン・ボルツマンの法則
! Ｍ＝σ・T^4(watt/㎡) T：絶対温度
! σ:5.67×10^-8(watt/㎡deg4)
!--------------------------------
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)=T(0)-σ*T(0)^4/J_BOX ! 熱放射式 シュテファン・ボルツマン
END IF
IF INT(TIM/1000)=TIM/1000 THEN
PRINT TIM;"時間経過"
! beep 1000,100
END IF
CALL CALC1 !計算部モジュール
CALL TEXT1 !文字表示モジュール
LET TIM=TIM+1
CALL MAXTEMP !最大の温度を検知
LOOP UNTIL MAX_TE<=50+K273 AND TIM>0
PRINT "終了！！"
!-------------------サブルーチン-------------------------
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)<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(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 J=0 TO 1
LET T(J)=4000+K273
NEXT J
END SUB
SUB CALC1
FOR J=1 TO depth-1
LET ⊿temp=T(J-1)-T(J) ! 境界面の温度差を計算
LET ⊿J=(d*1)*⊿temp*60*60 ! 1m 1時間に移動するエネルギー
LET T(J)=⊿J/J_BOX+T(J)
LET T(J-1)=-⊿J/J_BOX+T(J-1)
NEXT J
END SUB
SUB DRAWING
! ///～400℃赤/400～100℃黄/100～50℃緑/50℃以下青///
FOR J=0 TO depth
IF T(J)<=50+K273 THEN SET LINE COLOR 2 ! 白0黒1青2緑3赤4水5黄6
IF T(J)>50+K273 AND T(J)<=100+K273 THEN SET LINE COLOR 3
IF T(J)>100+K273 AND T(J)<=400+K273 THEN SET LINE COLOR 6
IF T(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(J)=40+K273
NEXT J
END SUB
SUB MAXTEMP
LET MAX_TP=50+K273
FOR J=0 TO depth
IF MAX_TP+K273<T(J)+K273 THEN
LET MAX_TP=T(J)
! PRINT MAX_TE
END IF
IF MAX_TP+K273<MAX_TE+K273 THEN LET MAX_TE=MAX_TP
NEXT J
END SUB
END