新しく発言する EXIT インデックスへ
隕石衝突後の熱伝導シュミレーション

  隕石衝突後の熱伝導シュミレーション 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
NHKの地球大進化という番組で
全海洋蒸発をした際の熱シュミレーションをしていました。
簡易版ができたので発表します。

前半
! ■■■■巨大隕石衝突後の熱伝導シュミレーション■■■■

!----------初期設定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 "MS ゴシック",B
SET WINDOW -1,WID,-depth,depth
SET TEXT COLOR 4 ! 白0黒1青2緑3赤4水5黄6
!--------------------------------

!----------初期設定1-------------
LET d=1.5 !土の熱伝導率 J/(sec・m・K) 温度差が1℃の時、1mの厚さを1秒に流れる熱量
! コンクリートの体積比熱(≒1.92倍)を参考にした仮の値
LET J_BOX=1.9*(4.18*10^6) ! 土の熱容量 J/(m3・K)
LET σ=5.67*10^(-8) ! (watt/udeg4)
! 放射の速度は、シュテファン・ボルツマンの法則
! M=σ・T^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)&"W/(m・deg)"
PLOT TEXT ,AT 0,depth*0.5 : " 土の熱容量(仮): "&STR$(J_BOX/1000)&"kJ/(m3・deg)"


IF TIM>24*365.25 THEN
SET TEXT COLOR 2
PLOT TEXT ,AT 0,depth*0.2 : "◆地球が冷える・・・熱放射式 W = σ・T^4"
SET TEXT COLOR 4
ELSEIF TIM<24*365.25 THEN
PLOT TEXT ,AT 0,depth*0.2 : " 地球は岩石蒸気で覆われる・・・ 約1年間"
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 !地表1mを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秒)を掛けてやらないと単位が合わないのですが。

それと、この式は、「厚さ1mの岩から黒体輻射で放出される分、温度が下がる」
とやっていますが、何か忘れていません?
初期温度4000度だから、海水は全部蒸発し大気となっていて、
大気(=海水)の持っている熱容量は、厚さ1mの岩より桁違いに大きいです。
だから、
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
実行時間を大幅に短縮することができますよ。

河川屋さんが専門的なアドバイスをされてますが、
私はそっち方面の知識はないので純粋に実行時間の短縮を図りました。

改良点は次の4点です。

1。print文での大量の出力を回避する
print文は出力結果が累積されるにつれ実行時間が遅くなります。
おそらく行送りに時間がかかるのでしょうが、数千行を超えると目に見えて遅くなります。
このプログラムの実行時間を遅くしている最大の要因です。
IF INT(TIM/1000)=TIM/1000 THEN
PRINT TIM;"時間経過"
END IF
この3行を削除するだけでかなりの時間を短縮できます。

2。ループを分割する
条件文(if文)は単純実行文より実行に時間がかかります。
数万回繰り返されるループの中に、条件に合致すると判定されることが1回しかないif文が何箇所かにあるのならば、ループを分割してしまったほうが良いでしょう。
ソースの記述は冗長になりますが時短を優先させました。

3。ループ内での演算を少なくする
できるだけループ内での演算回数を少なくするために、同じ演算はあらかじめ計算しておき定数化しました。
内部副プログラムの'CALC1'と'DRAWING'のループの中で、いかに演算を少なくするかが鍵です。

4。グラフィック画面では変更のある箇所だけ描画するようにする
これは、以前この掲示板で哲さんに教えていただいたテクニックです。(掲示板No.78)
clear文はどうしても画面がちらつきます。



改良したプログラムでは私のPCで、雨が降りはじめるまでに1時間4分、終了までに5時間18分かかりました。
なお、新しく設定した変数は次のとおりです。

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 "MS ゴシック",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 !土の熱伝導率 J/(sec・m・K) 温度差が1℃の時、1mの厚さを1秒に流れる熱量
! コンクリートの体積比熱(≒1.92倍)を参考にした仮の値
LET J_BOX=1.9*(4.18*10^6) ! 土の熱容量 J/(m3・K)
LET σ=5.67*10^(-8) ! (watt/udeg4)
LET dJB=(d*1)*60*60/J_BOX !! 定数化
LET σJB=σ/J_BOX !! 定数化
! 放射の速度は、シュテファン・ボルツマンの法則
! M=σ・T^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)&"W/(m・deg)"
PLOT TEXT ,AT 0,depth*0.5 : " 土の熱容量(仮): "&STR$(J_BOX/1000)&"kJ/(m3・deg)"
PLOT TEXT ,AT 0,depth*0.2 : " 地球は岩石蒸気で覆われる・・・ 約1年間"

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 : " 地球は岩石蒸気で覆われる・・・ 約1年間" !! 消去
SET TEXT COLOR 2
PLOT TEXT ,AT 0,depth*0.2 : "◆地球が冷える・・・熱放射式 W = σ・T^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 !地表1mを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


インデックスへ EXIT
新規発言を反映させるにはブラウザの更新ボタンを押してください。