Pythonで物理シミュ：放射性物質の崩壊

物理モデル

核子$N_0$が半減期$\mu_0$で$N_1$に崩壊するとき、$\lambda_0=\log2/\mu_0$として、
$$\frac{dN_0}{dt} = -\lambda_0 N_0$$
核子$N_1$が半減期$\mu_1$で$N_2$に崩壊するとき、$\lambda_1=\log2/\mu_1$として、
$$\frac{dN_1}{dt} = \lambda_0 N_0 -\lambda_1 N_1$$
以降、同様にして
$$\frac{dN_2}{dt} = \lambda_1 N_1 -\lambda_2 N_2$$
この連立常微分方程式をウラン238系列の値をもとに計算してみる。

import

import numpy as np
import scipy
from scipy.integrate import odeint

import matplotlib.pyplot as plt
%matplotlib inline

コード

NA = 6.02214085774 * 10 ** 23 # アボガドロ数
year = 1.0
day = year / 365.25
hour = day / 24.0
min = hour / 60.0
sec = min / 60.0
msec = 0.001 * sec

L2 = np.log(2.0)

U_238 = L2 / (4.468 * 10 ** 9 * year)
Th_234 = L2 / (24.1 * day)
Pa_234 = L2 / (1.17 * min)
U_234 = L2 / (2.455 * 10 ** 5 * year)
Th_230 = L2 / (7.538 * 10 ** 4 * year)
Ra_226 = L2 / (1600 * year)
Rn_222 = L2 / (3.824 * day)
Po_218 = L2 / (3.101 * min)
Pb_214 = L2 / (26.8 * min)
Bi_214 = L2 / (19.94 * min)
# Po_214 = L2 / (0.1643 * msec)
Pb_210 = L2 / (22.3 * year)
Bi_210 = L2 / (5.013 * day)
Po_210 = L2 / (183.76 * day)
Pb_206 = 1.0

L = (U_238, Th_234, Pa_234, U_234, Th_230, Ra_226, Rn_222, Po_218, Pb_214, Bi_214, Pb_210, Bi_210, Po_210, Pb_206)

def radio(N, t):
    return [
        -L[0]*N[0],
        L[0]*N[0] - L[1]*N[1],
        L[1]*N[1] - L[2]*N[2],
        L[2]*N[2] - L[3]*N[3],
        L[3]*N[3] - L[4]*N[4],
        L[4]*N[4] - L[5]*N[5],
        L[5]*N[5] - L[6]*N[6],
        L[6]*N[6] - L[7]*N[7],
        L[7]*N[7] - L[8]*N[8],
        L[8]*N[8] - L[9]*N[9],
        L[9]*N[9] - L[10]*N[10],
        L[10]*N[10] - L[11]*N[11],
        L[11]*N[11] - L[12]*N[12],
        L[12]*N[12] - L[13]*N[13],
        L[13]*N[13]
    ]

N_0 = [NA, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]

t = np.arange(0, 10*year, day)
v = odeint(radio, N_0, t)

plt.figure(figsize=(16, 10))
plt.ylim(ymin=1.0, ymax=10.0**24)
plt.xscale('linear')
plt.yscale('log')
plt.plot(t, v)
plt.legend([
    r'${}^{238}U$',
    r'${}^{234}Th$',
    r'${}^{234}Pa$',
    r'${}^{234}U$',
    r'${}^{230}Th$',
    r'${}^{226}Ra$',
    r'${}^{222}Rn$',
    r'${}^{218}Po$',
    r'${}^{214}Pb$',
    r'${}^{214}Bi$',
#    r'${}^{214}Po$',
    r'${}^{210}Pb$',
    r'${}^{210}Bi$',
    r'${}^{210}Po$',
    r'${}^{206}Pb$'
],fontsize=18) # 凡例
plt.show()
plt.close()

結果の前に

真面目に値を入れるとエラーが出て特に長期間でのシミューレーションができない。
ポロニウム214の半減期があまりに短すぎることが原因のようなので一旦コメントアウトした。

短期のシミュレーション

これは10年間を1日刻みでのシミュレーションだ。

放射性崩壊の系列では平衡に達すると各物質の存在比は半減期に比例する。
トリウム234とプロトアクチニウム234はこのスケールでは一瞬で平衡に達して、
ウラン238と平行なほぼ水平線になる。
ウラン234は半減期に比べて短いので平衡には達せず右肩上がりになってる。
トリウム230はウラン234にぴったりと追従するがグラフでは判別できない程度に傾きは大きい。
その下はラジウム226と鉛210でウラン234ではなくトリウム230の線に平行になる。
その下の他より傾きが大きそうなのが最終生成物の鉛206の線。
それ以外の核種は半減期が短すぎるのでラジウム226の線に平行に追従していきます。

中期のシミュレーション

これは200万年間を100年刻みでのシミュレーションだ。

ウラン238はこの期間ではびくとも減らない。
ウラン234とトリウム230以外は十分平衡に達している。
トリウム230の線はウラン234に接近していくが、
ウラン234の方がやや半減期が長いので微妙に届かないところで水平になって以降は平衡であるとみなせる。
70万年ぐらいで鉛206が第2位の存在量になっていることがわかる。

長期のシミュレーションの前に

年以下の半減期の核種は無視できるので除きましょう。

NA = 6.02214085774 * 10 ** 23 # アボガドロ数
myear = 1.0
year = 1.0 * myear / 1000000
L2 = np.log(2.0)

U_238 = L2 / (4.468 * 10 ** 9 * year)
U_234 = L2 / (2.455 * 10 ** 5 * year)
Th_230 = L2 / (7.538 * 10 ** 4 * year)
Ra_226 = L2 / (1600 * year)
Pb_210 = L2 / (22.3 * year)
Pb_206 = 1.0

L = (U_238, U_234, Th_230, Ra_226, Pb_210, Pb_206)

def radio(N, t):
    return [
        -L[0]*N[0],
        L[0]*N[0] - L[1]*N[1],
        L[1]*N[1] - L[2]*N[2],
        L[2]*N[2] - L[3]*N[3],
        L[3]*N[3] - L[4]*N[4],
        L[4]*N[4]
    ]

N_0 = [NA, 0.0, 0.0, 0.0, 0.0, 0.0]

t = np.arange(0, 10000.0*myear, 1.0*myear)
v = odeint(radio, N_0, t)

plt.figure(figsize=(16, 10))
plt.ylim(ymin=10.0**14, ymax=10.0**24)
plt.xscale('linear')
plt.yscale('log')
plt.plot(t, v)
plt.legend([
    r'${}^{238}U$',
    r'${}^{234}U$',
    r'${}^{230}Th$',
    r'${}^{226}Ra$',
    r'${}^{210}Pb$',
    r'${}^{206}Pb$'
],fontsize=18) # 凡例
plt.show()
plt.close()

長期のシミュレーション

これは100億年間を100万年刻みでのシミュレーションである。

ウラン238と平行に右下がりになる核種と鉛206がほぼ45億年で存在量1位になることが確認できる。