TomS |
Verfasst am: 19. Feb 2024 09:25 Titel: |
|
Ich habe jetzt mal die vollständige Zerfallskette
mit dem Differentialgleichungssystem
für den Vektor N nach zwei Methoden berechnet, 1. numerische Lösung, 2. Matrix-Exponential
Beide Methoden liefern bis auf numerische Fehler < 1e-3 übereinstimmende Ergebnisse. Anbei zunächst die Anzahl der jeweiligen Nuklide N, die Aktivitäten
sowie die freiwerdenden Energien
Als Anfangsbedingung setzte ich 1 mol Ra226 an. Hier der Python Code
Code: | # """ decay chain; example: Radium-226 attetion: works for linear decay chains only Author: Tom Date: 2024-02-15 """ from enum import Enum, auto import numpy as np import scipy.constants as sp_sc import scipy.linalg as sp_la import scipy.integrate as sp_in import matplotlib.pyplot as plt plt.rc('font', size=8) import mplcyberpunk plt.style.use("cyberpunk") # the solution ------------------------------ # using the matrix exponential: calculate the quantity, activity and energy for all nuclides as a function of time def calc_dec_chain(Nuclides, dec_coeffs, dec_energies, qty_nuclides_0, t_start, t_end, t_space): def init_dec_mat(Nuclides, dec_coeffs): """ restricted to linear decay chains """ dec_mat = np.zeros((len(Nuclides),len(Nuclides))) for m in range (len(Nuclides)): for n in range (len(Nuclides)): if (m == n): dec_mat[m,n] = -dec_coeffs[n] if (n > 0): dec_mat[m,n-1] = dec_coeffs[n-1] return dec_mat dec_mat = init_dec_mat(Nuclides, dec_coeffs) qty_sol_space = np.zeros((len(Nuclides),t_nsamples)) act_sol_space = np.zeros((len(Nuclides),t_nsamples)) eny_sol_space = np.zeros((len(Nuclides),t_nsamples)) for i in range(t_nsamples): qty_sol_space[:,i] = np.dot(sp_la.expm(t_space[i] * dec_mat), qty_nuclides_0) act_sol_space[:,i] = dec_coeffs * qty_sol_space[:,i] eny_sol_space[:,i] = dec_energies * dec_coeffs * qty_sol_space[:,i] return qty_sol_space, act_sol_space, eny_sol_space # using a 1st order system of differential equations: calculate the quantity, activity and energy for all nuclides as a function of time def solve_dec_chain(Nuclides, dec_coeffs, dec_energies, qty_nuclides_0, t_start, t_end, t_space): def def_decay_chain(t, q, c): """ restricted to linear decay chains """ dot_q = np.zeros_like(q) dot_q[0] = -c[0] * q[0] dot_q[-1] = c[-2] * q[-2] for n in range(1, len(q)-1): dot_q[n] = c[n-1] * q[n-1] - c[n] * q[n] return dot_q t_span = (t_start, t_end) sol = sp_in.solve_ivp(def_decay_chain, t_span, qty_nuclides_0, args=(dec_coeffs,), t_eval=t_space, method='BDF') return sol.y, None, None # extract alpha and beta channels, calculate sums for alpha and beta def sum_alpha_beta_channels(sol_space): alpha = np.zeros(t_nsamples) beta = np.zeros(t_nsamples) for i in range(t_nsamples): alpha[i] += (sol_space[Nuclides.Ra_226.value, i] + sol_space[Nuclides.Rn_222.value, i] + sol_space[Nuclides.Po_218.value, i] + sol_space[Nuclides.Po_214.value, i] + sol_space[Nuclides.Po_210.value, i]) beta[i] += (sol_space[Nuclides.Pb_214.value, i] + sol_space[ Nuclides.Bi_214.value, i] + sol_space[ Nuclides.Pb_210.value, i] + sol_space[ Nuclides.Bi_210.value, i]) return alpha, beta # plot functions ------------------------------ def subplot_quantities(t_space_years, qty_sol_space): for n in Nuclides: plt.plot(t_space_years, qty_sol_space[n.value], lw=0.75, ls='-', label=n.name) return def subplot_activities(t_space_years, qty_sol_space): for n in Nuclides: if (n.value == len(Nuclides) - 1): break plt.plot(t_space_years, act_sol_space[n.value], lw=0.75, ls='-', label=n.name) return def subplot_energies(t_space_years, qty_sol_space): for n in Nuclides: if (n.value == len(Nuclides) - 1): break plt.plot(t_space_years, eny_sol_space[n.value], lw=0.75, ls='-', label=n.name) return def subplot_channels(t_space_years, radiation): plt.plot(t_space_years, radiation[0], lw=0.75, ls='-', label='alpha') plt.plot(t_space_years, radiation[1], lw=0.75, ls='-', label='beta') return def plot_frame(t_space, qty_sol_space, x_label, y_label, title, subplot_func): plt.figure() plt.xscale('log') plt.yscale('log') plt.xlabel(x_label) plt.ylabel(y_label) plt.title(title) subplot_func(t_space/sp_sc.year, qty_sol_space) plt.legend() plt.show() return # the nuclides in the Ra-226 decay chain, their decay half-times [seconds], their decay energies [MeV], and their quantity ------------------------------ Nuclides = Enum('Nuclides', ['Ra_226', 'Rn_222', 'Po_218', 'Pb_214', 'Bi_214', 'Po_214', 'Pb_210', 'Bi_210', 'Po_210', 'Pb_206'], start=0) dec_times = np.array([1600 * sp_sc.year, 3.82 * sp_sc.day, 3.098 * sp_sc.minute, 27.06 * sp_sc.minute, 19.9 * sp_sc.minute, 0.0001643, 22.20 * sp_sc.year, 5.012 * sp_sc.day, 138.376 * sp_sc.day, np.inf]) dec_coeffs = np.log(2) / dec_times dec_energies = np.array([4.7, 5.5, 6.0, 0.3, 1.0, 7.7, 0.006, 0.3, 5.3, 0.0]) qty_nuclides_0 = np.zeros_like(dec_times) qty_nuclides_0[0] = sp_sc.Avogadro # driver code ------------------------------ # set time span t_start = 1.0 * sp_sc.hour t_end = 10000.0 * sp_sc.year t_nsamples = 1000 t_space = np.logspace(np.log10(t_start), np.log10(t_end), t_nsamples) # solve equations qty_sol_space, act_sol_space, eny_sol_space = calc_dec_chain(Nuclides, dec_coeffs, dec_energies, qty_nuclides_0, t_start, t_end, t_space) qty_sol_space_2, _, _ = solve_dec_chain(Nuclides, dec_coeffs, dec_energies, qty_nuclides_0, t_start, t_end, t_space) # sum for alpha and beta channels act_alpha, act_beta = sum_alpha_beta_channels(act_sol_space) eny_alpha, eny_beta = sum_alpha_beta_channels(eny_sol_space) # plot plot_frame(t_space, qty_sol_space, 'time [years]', 'quantities [N]','quantity', subplot_quantities) plot_frame(t_space, qty_sol_space_2, 'time [years]', 'quantities [N]','quantity', subplot_quantities) plot_frame(t_space, qty_sol_space_2/qty_sol_space-1.0, 'time [years]', 'relative error','relative error', subplot_quantities) plot_frame(t_space, qty_sol_space, 'time [years]', 'activities [Bq]','activity', subplot_activities) plot_frame(t_space, qty_sol_space, 'time [years]', 'energies [MeV]','energy', subplot_energies) plot_frame(t_space, (act_alpha, act_beta), 'time [years]', 'activities [Bq]','activity', subplot_channels) plot_frame(t_space, (eny_alpha, eny_beta), 'time [years]', 'energies [MeV]','energy', subplot_channels) |
|
|