Autor	Nachricht
TomS	Verfasst am: 19. Feb 2024 10:14    Titel: Anbei noch weitere Ergebnisse zu den Aktivitäten und Energien für alpha- und beta-Strahlung. D.h. z.B.
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)
TomS	Verfasst am: 11. Feb 2024 13:40    Titel: Die Gesamtaktivität für alpha-und beta-Strahlung purzelt auch heraus; gamma-Strahlung nicht, das müsste man separat lösen.
Lorz	Verfasst am: 11. Feb 2024 10:02    Titel: TomS hat Folgendes geschrieben: Lorz hat Folgendes geschrieben: Aber ich verstehen nicht den Doppelindex, an der Zerfallskonstanten. Im Allgemeinen kann ein Nuklid k in verschiedene Tochternuklide zerfallen, also in m1, m2 … – siehe Abbildung. Ist das nicht der Fall, bzw. setzt man näherungsweise an, dass je Nuklid ein Zerfall dominiert, d.h. dass eine Zerfallsreihe vorliegt, dann hast du immer nur zwei Einträge für einen Zerfall: km mit k<m: Nuklid k wird weniger mk mit k<m: Nuklid m wird um den selben Betrag mehr Beide Einträge in der Matrix haben unterschiedliches Vorzeichen, d.h. letztlich kommt je Nuklid k nur eine Zahl vor, die Zerfallskonstante k. Und diese steht genau zweimal da, dämlich bei km mit m=k+1: Nuklid k wird weniger mk mit m=k+1: Nuklid m=k+1 wird um den selben Betrag mehr Ah, OK, verstehe, mehrere Zerfallsmöglichkeiten pro Nuklid, mehrere Zerfallskonstanten. Danke für das Bild, ja, klar, Zerfallsreihen mit mehr als einem Weg kenne ich. Aber nu mach mal ruhig Deinen Urlaub! Ich probiere mich mal an Deine dargebotenen Möglichkeiten, sollte ja recht knapp werden meine Matrix, wenn ich erstmal nur die bis ink. der 3. Generation gehe und dabei nur eine Zerfallskonstante pro Radionuklid annehme. Ansonsten interessiert mich schon die weitere Entwicklung der Gesamtaktivität. Also mein ursprüngliches Ziel war es, zu verstehen, wie es zu einem säkularen Gleichgewicht kommt, bei dem das konstante Plateau recht lang ist.
TomS	Verfasst am: 11. Feb 2024 08:00    Titel: Hier ist das sehr ausführlich dargestellt: https://arminda.whitman.edu/_flysystem/fedora/2022-09/A_method_to_decouple_the_equations_describing_radioactive_decay.pdf
TomS	Verfasst am: 10. Feb 2024 20:54    Titel: Lorz hat Folgendes geschrieben: Aber ich verstehen nicht den Doppelindex, an der Zerfallskonstanten. Im Allgemeinen kann ein Nuklid k in verschiedene Tochternuklide zerfallen, also in m1, m2 … – siehe Abbildung. Ist das nicht der Fall, bzw. setzt man näherungsweise an, dass je Nuklid ein Zerfall dominiert, d.h. dass eine Zerfallsreihe vorliegt, dann hast du immer nur zwei Einträge für einen Zerfall: km mit k<m: Nuklid k wird weniger mk mit k<m: Nuklid m wird um den selben Betrag mehr Beide Einträge in der Matrix haben unterschiedliches Vorzeichen, d.h. letztlich kommt je Nuklid k nur eine Zahl vor, die Zerfallskonstante k. Und diese steht genau zweimal da, dämlich bei km mit m=k+1: Nuklid k wird weniger mk mit m=k+1: Nuklid m=k+1 wird um den selben Betrag mehr
TomS	Verfasst am: 10. Feb 2024 20:50    Titel: Wenn’s dich interessiert, hier der Einzeiler in Python: Code:  qty_sol_space[:,i] = np.dot(la.expm(t_space[i] * dec_mat), qty_nuclides_0) Das löst das DGL-System für beliebige Zeiten bei gegebener Anfangsbedingung. Außenrum noch 50 Zeilen Code zum Befüllen von Lambda und für die Graphik. Ich hätte das für die Zerfallsreihe von Radium fertig und würde es dir nach dem Urlaub schicken.
Lorz	Verfasst am: 10. Feb 2024 20:44    Titel: TomS hat Folgendes geschrieben: Man setzt: Im der k-ten Zeile steht die Änderung des k-ten Nuklids, in der Klammer mit dem Minuszeichen die Zerfälle in die folgenden Tochternuklide k+1, k+2 … , in der Klammer mit dem Pluszeichen die Zerfälle aus den Mutternukliden k-1, k-2 … Das kann man als Gleichung mit einer Zerfallsmatrix Lambda und einem Vektor N auffassen: Auf den Diagonalelementen kk stehen die Summen der Zerfallskonstanten in die Tochternuklide. In einem Dreieck die einzelnen Matrixelemente km mit m>k mit den Zerfallskonstanten von k nach m. Das zweite Dreieck ist Null, keine Zerfälle von m nach k für m>k. In einfachen Fällen d.h. für wenige beteiligte Nuklide bzw. viele Nullen man das mittels Diagonalisierung von Lambda analytisch lösen. Ansonsten füttert man das in einen Computer. OK, damit könnte ich ja tatsächlich beliebig lange Zerfallsreihen behandeln. Aber ich verstehen nicht den Doppelindex, an der Zerfallskonstanten. Nach meiner Idee habe ich für jedes der k Generationen an Nukliden eine Zerfallskonstante . Ich kann diese vielen nicht zuordnen.
Lorz	Verfasst am: 10. Feb 2024 20:39    Titel: 3.Generation hat Folgendes geschrieben: Lorz hat Folgendes geschrieben: Was ich nicht verstehe - warum "explodiert" die Anzahl, wenn ? Da explodiert nichts,weil die Gesamtmenge hier immer N1(0) beträgt Man hat hier eine Grenzwertrechnung So etwa Ah, ich habe nicht Summe bzw Differenz aus beiden Brüche betrachtet. Einzeln gesehen geht der Bruch ja schon gegen unendlich. Ist daher hier nicht das Grenzwertproblem, dass entsteht? 3.Generation hat Folgendes geschrieben: Lorz hat Folgendes geschrieben: Und warum kommt in der Lösung für kein vor? Streng genommen kommt vor Wenn N3 nicht weiter zerfällt ist und damit Kann man denn noch angenehm erweitern auf eine Halbwertszeit 3. Generation kleiner als unendlich?
TomS	Verfasst am: 10. Feb 2024 20:03    Titel: Man setzt: Im der k-ten Zeile steht die Änderung des k-ten Nuklids, in der Klammer mit dem Minuszeichen die Zerfälle in die folgenden Tochternuklide k+1, k+2 … , in der Klammer mit dem Pluszeichen die Zerfälle aus den Mutternukliden k-1, k-2 … Das kann man als Gleichung mit einer Zerfallsmatrix Lambda und einem Vektor N auffassen: Auf den Diagonalelementen kk stehen die Summen der Zerfallskonstanten in die Tochternuklide. In einem Dreieck die einzelnen Matrixelemente km mit m>k mit den Zerfallskonstanten von k nach m. Das zweite Dreieck ist Null, keine Zerfälle von m nach k für m>k. In einfachen Fällen d.h. für wenige beteiligte Nuklide bzw. viele Nullen man das mittels Diagonalisierung von Lambda analytisch lösen. Ansonsten füttert man das in einen Computer.
3.Generation	Verfasst am: 10. Feb 2024 19:44    Titel: Lorz hat Folgendes geschrieben: Was ich nicht verstehe - warum "explodiert" die Anzahl, wenn ? Da explodiert nichts,weil die Gesamtmenge hier immer N1(0) beträgt Man hat hier eine Grenzwertrechnung So etwa Lorz hat Folgendes geschrieben: Und warum kommt in der Lösung für kein vor? Streng genommen kommt vor Wenn N3 nicht weiter zerfällt ist und damit
Lorz	Verfasst am: 10. Feb 2024 18:22    Titel: 3.Generation hat Folgendes geschrieben: Lorz hat Folgendes geschrieben: Ich verstehe allerdings gerade nicht, wieso bereits bei ein auftaucht. Man löst diese Aufgabe am besten mit einer Eigenwertrechnung und der Eigenvektor von ist Es gilt zu jedem Zeitpunkt Lorz hat Folgendes geschrieben: Auch bin ich skeptisch, was die Dimension/Einheiten der Summanden in und betrifft. Also die Zerfallskonstante ist ja nicht einfach dimensionslos (hat doch ) . Oder wie verhält sich dies? Gut erkannt Man könnte die Einheit ausgleichen indem man s bei C1 dazuschreibt OK, verstehe, man hat hier viel mehr Flexibilität. Ich kann mit einer beliebigen Anfangssituation starten. Bei meinem vereinfachten Beispiel (alle Mutternuklide warten mit dem Zerfallen bis , dann gilt Weiter müsste nach der obersten Gleichung gelten Soll denn die Zerfallskonstante sein? Und entsprechend ? Dann wäre Mit erhalten Und somit für Ich prüfe mal für mich die Plausibität von : Also für hebt sich alles weg, , das passt. Nun nehme ich mal den Fall, dass winzig klein ist, dann dominierend gegenüber in der Differenz und kürzt sich im Bruch weg. Wir erhalten bzw Auch das macht Sinn, die Mutternuklide wandeln sich quasi sofort um, und und nur noch bestimmt den Verlauf. Was ich nicht verstehe - warum "explodiert" die Anzahl, wenn ? Und warum kommt in der Lösung für kein vor?
Lorz	Verfasst am: 10. Feb 2024 18:16    Titel: Ich probiere jetzt mal selbst. Also die Gleichung setze ich jetzt erst mal für die zweite Generation an und schaue, ob das was vernünftiges bei rauskommt. Umformungen ...ähhm, ja...
3.Generation	Verfasst am: 10. Feb 2024 18:01    Titel: Lorz hat Folgendes geschrieben: Ich verstehe allerdings gerade nicht, wieso bereits bei ein auftaucht. Man löst diese Aufgabe am besten mit einer Eigenwertrechnung und der Eigenvektor von ist Es gilt zu jedem Zeitpunkt Lorz hat Folgendes geschrieben: Auch bin ich skeptisch, was die Dimension/Einheiten der Summanden in und betrifft. Also die Zerfallskonstante ist ja nicht einfach dimensionslos (hat doch ) . Oder wie verhält sich dies? Gut erkannt Man könnte die Einheit ausgleichen indem man s bei C1 dazuschreibt
Lorz	Verfasst am: 10. Feb 2024 16:55    Titel: @ Gast Danke! Ich verstehe allerdings gerade nicht, wieso bereits bei ein auftaucht. Auch bin ich skeptisch, was die Dimension/Einheiten der Summanden in und betrifft. Also die Zerfallskonstante ist ja nicht einfach dimensionslos (hat doch ) . Oder wie verhält sich dies?
3.Generation	Verfasst am: 10. Feb 2024 16:50    Titel: Lorz hat Folgendes geschrieben: Also wenn mir jemand einfach so die Lösung für Generation 3 sagen kann, bin ich happy :-) Das funktioniert nur wenn Bei wird die Sache wesentlich komplizierter
Lorz	Verfasst am: 10. Feb 2024 15:06    Titel: @TomS Danke für die Antwort. Dachte ich mir schon, dass ich das nicht so einfach ansetzen kann. Bzw, das war mein erster Gedanke. Danach habe ich gedacht - ist ja doch nicht kompliziert. Aber der erste Gedanke war wohl schon richtig. Konkret zu Ergebnissen: Ist die Anzahl noch korrekt? Denn diese Art von Kurve finde ich auch in der Literatur. Von der Berechnung beliebig vieler Generationen sehe ich dann gerne ab. Aber kann mir jemand die Funktion noch für die dritte Generation nennen? Ein bisschen versuche ich mitzuhelfen: Also der zukünftige Bestand (zur Zeit ) der Generation 3 ergibt sich aus dem Zuwachs durch den momentanen Abbau von Generation 2 (gleichbedeutend mit Umwandlung in die Generation 3) innerhalb zur Zeit und dem Abbau des schon vorhandenen Bestands von : Soweit korrekt? Da der Anfangswert für die erste Generation als gegeben vorausgesetzt ist, könnte ich meinem Verständnis nach, die 2. Generation berechnen, denn und sind ja bekannt. Und iterativ sollten dann auch alle anderen Generationen berechenbar sein. Aber ich weiß halt nicht mehr, wie das geht. Also ich kann mich an so ein paar Verfahren (Picard-Lindelöff oder so) erinnern aus der Analysis. Theo-Physik ist leider fast komplett vergessen bei mir. Also wenn mir jemand einfach so die Lösung für Generation 3 sagen kann, bin ich happy :-)
TomS	Verfasst am: 10. Feb 2024 14:38    Titel: Man muss ein gekoppeltes lineares Gleichungssystem ansetzen, das alle Zerfallskanäle berücksichtigt. Wenn man immer nur einen Zerfallskanal je Nuklid hat, wird das ziemlich einfach, da man eine DGL nach der anderen lösen kann.
Lorz	Verfasst am: 10. Feb 2024 14:21    Titel: Gesamtradioaktivität (Mutter-+Töchteraktivitäten) Hallo! Ich interessiere mich für die Berechnung der Gesamtaktivität eines Nuklids. Dazu betrachte ich zunächst die Anzahlen der Generationen. Mit einer gewissen Anfangsanzahl ergibt sich die momentane Anzahl an Mutternukliden (also die 1. Generation) durch Für zweite Generation sollte dann gelten Oder als t-N-Diagramme für zB ein Verhältnis (siehe grüne und rote Kurve) Wenn ich nun noch die dritte Generation betrachten will, so würde ich eigentlich diese Funkton ansetzen Aber wenn ich alle drei zeichne, so sieht das so aus (siehe blaue Kurve) Ich könnte jetzt ja einfach mal die blaue Kurve umklappen. Sehe aber vom Ansatz meines Funktionsterms keinen Grund dafür. Weiß jemand Rat? Edit: Ah, jetzt checke ich das Problem ist Null... Mmmhm, aber was muss ich dann einsetzen? Ist die Berechnung doch komplizierter? Muss ich da mit Differenzialgleichungen ran? Stimmt überhaupt die rote Kurve?

Powered by phpBB © 2001, 2005 phpBB Group