読者です 読者をやめる 読者になる 読者になる

倭算数理研究所

科学・数学・学習関連の記事を、「倭マン日記」とは別に書いていくのだ!

2.10 原子核の崩壊 : 問題 2.8 原子核の崩壊 a

計算物理学入門 シミュレーション

前回までは単一核種の系を扱ってましたが、今回からは複数の核種の系、すなわちある核種が放射性崩壊してできた核種がさらに放射性崩壊するような系を見ていきます。

放射性崩壊する2核種の系

最初の種類の原子核の個数を { N_1(t) }、崩壊によってできる放射性原子核の個数を { N_2(t) } とし、同様にそれぞれの崩壊定数を { \lambda_1,\, \lambda_2 } とします。 このとき、{ N_1(t),\,N_2(t) } の満たす微分方程式は以下のようになります:

  { \displaystyle
\begin{align*}
    \frac{dN_1}{dt} &= -\lambda_1 N_1 \\[4mm] \frac{dN_2}{dt} &= \lambda_1 N_1 - \lambda_2 N_2
\end{align*}
}

この系は解析的に解けて、初期条件が { N_1(0) = N_{10},\,N_2(0) = N_{20} } のとき

  { \displaystyle
\begin{align*}
    N_1(t) &= N_{10}\,e^{-\lambda_1 t} \\[2mm]
    N_2(t) &= \frac{\lambda_1}{\lambda_2-\lambda_1}N_{10}\,e^{-\lambda_1 t}
                + \left(N_{20} - \frac{\lambda_1}{\lambda_2-\lambda_1}N_{10}\right)\,e^{-\lambda_2 t}
\end{align*}
}

特に、{ N_{20} = 0 } のとき、つまり最初2種目の原子核がないときは

  { \displaystyle
\begin{align*}
    N_1(t) &= N_{10}\,e^{-\lambda_1 t} \\[2mm]
    N_2(t) &= \frac{\lambda_1}{\lambda_2-\lambda_1}N_{10}\left(e^{-\lambda_1 t} - e^{-\lambda_2 t}\right)
\end{align*}
}

となります。 導出はこちらを参照。

問題

a.

  • 76Kr は電子を捕獲して半減期14.8時間で 76Br に崩壊する
  • 76Br は電子を捕獲するか陽電子を放出して半減期16.1時間で 76Se に崩壊する
  • 試料には最初に純粋の 76Kr が1g含まれていたとする

プログラムを修正して 76Kr と 76Se の量の1週間の時間変化を計算せよ。

b.

時間の単位に分を使うと、ある程度の量の 28Mg が崩壊するまでには多数回の反復計算が必要になる。 計算を速めるには、どのような簡単化の過程を置いたらよいか。

c. 211Rn は、下図に示すように、2つに枝分かれした過程で崩壊する。 必要ならばなんらかの近似をし、各同位体元素の量を時間の関数として計算せよ。 試料は初め 1μg の 211Rn からなっていたとする。

今回は問題 a をします。 問題 a に出てくる元素をもう少し詳しく見ると下表の通り:

元素記号 元素名 原子番号 周期・族 Wikipeida
Kr クリプトン (krypton) 36 第4周期18族 wikipedia:クリプトン
Br 臭素 (bromine) 35 第4周期17族 wikipedia:臭素
Se セレン (selenium) 34 第4周期16族 wikipedia:セレン

物理系

問題 a の物理系の簡単なドメイン分析をすると

  • 状態変数
    • 時間 { t } : double
    • Kr, Br, Se の核数 { N_\textrm{Kr},\,N_\textrm{Br},\,N_\textrm{Se} } : それぞれ double
  • シミュレーション・パラメータ
    • Kr, Br の崩壊定数 { \lambda_\textrm{Kr},\,\lambda_\textrm{Br} } : それぞれ double

ってな感じですね。 今回は時間 { t }BigDecimal ではなく double で定義してます。 これは後ほど。 で、この系が従う時間発展の微分方程式

  { \displaystyle
\begin{align*}
    \frac{dN_\textrm{Kr}}{dt} &= -\lambda_\textrm{Kr} N_\textrm{Kr} \\[4mm]
    \frac{dN_\textrm{Br}}{dt} &= \lambda_\textrm{Kr} N_\textrm{Kr} - \lambda_\textrm{Br} N_\textrm{Br} \\[4mm]
    \frac{dN_\textrm{Se}}{dt} &= \lambda_\textrm{Br} N_\textrm{Br}
\end{align*}
}

となります。 Se の核数の微分方程式も載せましたが大丈夫でしょうかね?

時間の単位

2つの核種の崩壊定数が異なるので、半減期も違ってきます。 半減期が異なる場合、短い方の半減期と同程度の大きさの単位で時間を測ると都合がよいとのことなので、シミュレーションを書く前に「2.10 原子核の崩壊 : 問題 2.7 単一核種の崩壊 c」でやったように時間の単位を変更してやりましょう*2。 問題より、各核種の半減期

  { \displaystyle
\begin{align*}
    T_{1/2\, \textrm{Kr}} &= 14.8 \,\textrm{h}, & T_{1/2\, \textrm{Br}} = 16.1 \,\textrm{h}
\end{align*}
}

なので、短い方の Kr の半減期を単位にして時間を測りましょう。 その時間を { \tau } (タウ)とすると

  { \displaystyle
\begin{align*}
    \tau &= \frac{t}{T_{1/2\,\textrm{Kr}}} = \frac{\lambda_\textrm{Kr}}{\log 2}t &
    \left(\because T_{1/2 \, \textrm{Kr}} = \frac{\log 2}{\lambda_\textrm{Kr}}\right)
\end{align*}
}

と定義されます。 { t } から { \tau } への書き換えは、以前行ったのと同じような式の書き換えなので途中は省略して、結果だけを書くと

  { \displaystyle
\begin{align*}
    \frac{dN_\textrm{Kr}}{d\tau} &= -N_\textrm{Kr} \log 2 \\[4mm]
    \frac{dN_\textrm{Br}}{d\tau} &= \left(N_\textrm{Kr} - \lambda N_\textrm{Br}\right)\log 2 \\[4mm]
    \frac{dN_\textrm{Se}}{dt} &= \lambda N_\textrm{Br}\log 2
\end{align*}
}

となります。 ただし

  { \displaystyle
\begin{align*}
    \lambda = \frac{\lambda_\textrm{Br}}{\lambda_\textrm{Kr}} = \frac{T_{1/2\,\textrm{Kr}}}{T_{1/2\,\textrm{Br}}}
\end{align*}
}

です。 ふむ、結構綺麗にまとまりましたね。

シミュレーション・スクリプト

以上を踏まえてシミュレーション・スクリプト。 シミュレーションの反復条件は1週間と指定されてますが、Kr の半減期で測った時間ではキレイな数値にならないので、今までのように時間を BigDecimal にする意味はないでしょう。 したがって、時間も double 型にします。 また、tau と書くのは面倒なので、この tau を改めて t としてスクリプトを書くことにします。

@GrabResolver('https://github.com/waman/groops-core/tree/master/repo')
@Grab('org.waman.groops:groops-core:0.2-beta')

import org.waman.groops.builder.SimulationBuilder
import org.waman.groops.simulation.system.*

double tKr = 14.8d, tBr = 16.1d    // 半減期
double aweek = 24 * 7 / tKr    // 一週間

new SimulationBuilder().simulation(type:'composite'){    // 複合シミュレーション
    simulation{
        system(new RadioactiveDecaySystem(), dt:0.01d, lambda:tKr/tBr)
        init(t:0d, nKr:1d, nBr:0d, nSe:0d)

        count(total:10)    // 10回の反復
    }

    // これらの反復条件、データ出力は10回の反復ごとに行われる(複合シミュレーション)
    continueWhile{ it.t <= aweek }
    console(outputTiming:FINAL)
    png(fileName:'problem2_8a.png', overwrite:true){
        chart(title:'Problem 2.8 a', domainLabel:"Time [$tKr hour]", rangeLabel:'Mass [g]'){
            point(label:'Kr', x:'t', y:'nKr')
            point(label:'Br', x:'t', y:'nBr')
            point(label:'Se', x:'t', y:'nSe')
        }
    }
}.simulate()

class RadioactiveDecaySystem extends AnnotationObservableSystem{

    @State double t, nKr, nBr, nSe
    @Parameter double dt, lambda

    static double LOG2 = Math.log(2d)

    // 今回はオイラー法で。
    @Override
    void evolve(Map<String, Object> params) {
        double dnKr =          nKr * LOG2 * dt
        double dnBr = lambda * nBr * LOG2 * dt

        nKr -= dnKr
        nBr -= dnBr - dnKr
        nSe += dnBr

        t += dt
    }
}
  • problem2_8a.groovy 直
  • 複合シミュレーション (compsite simulation) によって、10回繰り返すごとにデータを出力するようにしています。

シミュレーション結果は以下のようになります:


おー、シミュレーションっぽくなってきたね。 Kr の崩壊は単一核の場合と同じです。 { t = 1 } で初期値の半分 (0.5) になってますね。 一週間経つと、ほとんどの核が最終的な Se に崩壊してしまってます。 

計算物理学入門

計算物理学入門

  • 作者: ハーベイゴールド,ジャントボチニク,Harvey Gould,Jan Tobochnik,鈴木増雄,石川正勝,溜渕継博,宮島佐介
  • 出版社/メーカー: ピアソンエデュケーション
  • 発売日: 2000/12
  • メディア: 単行本
  • 購入: 1人 クリック: 28回
  • この商品を含むブログ (45件) を見る
プログラミングGROOVY

プログラミングGROOVY

*1:『計算物理入門』では Mg ではなく Mn になってますが、まぁ間違いでしょう。

*2:同程度でいいので、正確にこの半減期を単位にする必要はないかもしれませんが。