DACH - Jahrestagung 2004 Salzburg

ZfP in Forschung, Entwicklung und Anwendung

Start > Beiträge > Vorträge > Simulation 1: Print

MONTE-CARLO-SIMULATION DER STRAHLUNGSAUSBREITUNG AN REALEN BAUTEILGEOMETRIEN

M. Zhukovsky, S. Podoliako
Keldysch Institut für Angewandte Mathematik, Moskau, Russland
G.-R. Jaenisch, C. Bellon
Bundesanstalt für Materialforschung und -prüfung (BAM), Berlin
Kontakt: Dr.-Ing. Carsten Bellon

Einleitung

Eine wirtschaftliche und verlässliche Prüfung komplexer Bauteile ist für den sicheren Betrieb von industriellen Anlagen und Baugruppen von entscheidender Bedeutung. Eine der klassischen Prüfmethoden stellt die Durchstrahlungsprüfung dar. Bei der Planung und Optimierung komplexer Prüfaufgaben kommt zunehmend auch die Computersimulation zum Einsatz. Entsprechend muss die Computersimulation alle für die Prüfaufgabe entscheidenden Eigenschaften der Prüfmethode mit ausreichender Genauigkeit widerspiegeln. Bei der Betrachtung der Durchstrahlungsprüfung schließt die Modellierung die Beschreibung der Strahlenquelle, der Wechselwirkung von Strahlung und Materie und des Detektionsprozesses sowie die geometrische Bauteilrepräsentation ein.

Standard-Simulationsprogramme für die Durchstrahlungsprüfung basieren auf analytischen Modellen unter Verwendung des Schwächungsgesetzes und von Aufbau-Faktoren zur Beschreibung der Wechselwirkung von Strahlung mit Materie. Die Einführung von Aufbau-Faktoren zur Beschreibung der Streustrahlung erfasst dabei nur eine einheitliche Kontrastminderung im Durchstrahlungsbild und nicht die in Abhängigkeit von der Prüfteilgeometrie lokal unterschiedlichen Streustrahlungsanteile. Diese Vereinfachung ist für viele Anwendungsbereiche wie z.B. bei der Schweißnahtprüfung ausreichend. Aber nur eine detaillierte Beschreibung der zugrunde liegenden Wechselwirkungsmechanismen ist in der Lage z.B. Mottling oder andere in der Praxis zu beobachtende Effekte wiederzugeben.

Der Monte-Carlo-Code MCNP [1] schließt die Behandlung aller bekannten Wechselwirkungsmechanismen die zur Bildentstehung beitragen ein. Es werden Wechselwirkungen von Photonen (Absorption, inkoherente und koherente Streuung einschließlich Bindungseffekte, Paarbildung) und Elektronen (Bremsstrahlungserzeugung, etc.) beschrieben. Allerdings ist die Beschreibung komplexer Prüfteilgeometrien sehr schwierig. Daher ist es Anliegen des vorliegenden Beitrages einen neuen Ansatz zur direkten Kopplung der Beschreibung des Photonentransports mit einer Oberflächen-Beschreibung komplexer Geometrien, wie sie aus dem CAD bekannt ist, zu diskutieren.

Wechselwirkungsmechanismen

Das Softwarepaket MCNP [1], das hier als Referenz betrachtet wird, behandelt Neutronen-, Photonen- und Elektronen-Wechselwirkungen. Im Folgenden werden nur Wechselwirkungen von Photonen betrachtet (s. Abb. 1). Die physikalische Behandlung des Photonentransports schließt den Photoeffekt sowie die Rayleigh- und Compton-Streuung ein. Durch eine Beschränkung auf Photonenenergien unter 1 MeV tritt die Paarbildung hier nicht auf. Im Vergleich zu MCNP werden weiterhin Sekundäreffekte wie Röntgenfluoreszenz und Elektronentransport vernachlässigt.


Abb 1: Wechselwirkungsmechanismen der Photonen.

Als Photoeffekt wird die Absorption eines einfallenden Photons bezeichnet, dessen Energie E vollständig an die Elektronen der Atomhülle abgegeben wird. Dabei kommt es zur Anregung eines Elektrons mit einer Bindungsenergie e < E oder zur Emission eines Photoelektrons mit einer kinetischen Energie von E-e. Werden bei genügend hoher Energie der anregenden Photonen innere Niveaus ionisiert, so wird der Photoeffekt von der Emission eines Fluoreszenzphotons oder eines Auger-Elektrons begleitet, die als Folge innerer Rekombinationsprozesse entstehen. In dieser Untersuchung wird die Emission eines Fluoreszenzphotons genauso wie die Verfolgung der entstehenden Elektronen vernachlässigt.

Zur Behandlung der Compton-Streuung sind der Streuwinkel q und die neue Energie E' des Photons zu bestimmen. Die Analyse der Compton-Streuung ergibt, dass die Energie des gestreuten Photons immer kleiner als die Energie des Primärphotons E>E'. Der Energieübertrag vom Primärphoton auf das Compton-Elektron hängt nur vom Streuwinkel q ab. Mit wachsendem Streuwinkel wird auch ein Anwachsen des Energieübertrags beobachtet, der sein Maximum für q=180° erreicht. Der differentielle Streuquerschnitt setzt sich aus dem Streuquerschnitt nach Klein und Nishina sKN sowie der nichtrelativistischen inkohärenten Streufunktion I(Z,E,q) [2] zusammen, um Elektronenbindungseffekte zu berücksichtigen:

(1)

mit

(2)

Dabei bezeichnet r0 = e2 / m0 c2 = 2.818 10-15 m den klassischen Elektronenradius mit der Elementarladung e und der Ruheenergie des Elektrons m0, sowie e = E / m0 c2 = E / 511 keV. Die Berücksichtigung der Elektronenbindungseffekte führt zu einer Abnahme des Streuquerschnittes in Vorwärtsrichtung für kleine Energien E und Atome mit großer Ordnungszahl Z. Bild 2 stellt beispielhaft diesen Effekt für Aluminium und Photonenenergien von 40 keV und 100 keV dar.

Abb 2: Vergleich von Klein-Nishina- und Compton-Streuquerschnitten für Aluminium bei Photonenenergien von 40 keV (links) and 100 keV (rechtes).

Bei der kohärenten Streuung, die auch als Rayleigh-Streuung bezeichnet wird, wird das einfallende Photon an einem Elektron in der Atomhülle gestreut, ohne Energie zu verlieren bzw. den Energiezustand des Atoms zu ändern. Im Gegensatz zur Compton-Streuung wird kein Elektron produziert, das weiterverfolgt werden müsste. Um das gestreute Photon weiter verfolgen zu können, ist es lediglich notwendig, den Streuwinkel q zu bestimmen. Für diesen Fall bestimmt sich der Streuquerschnitt aus dem Thompson-Querschnitt sThompson und dem relativistischen Formfaktor C(Z,E,q) [3] zur Berücksichtigung von Bindungseffekten in der Atomhülle

(3)

mit dem Thompson-Querschnitt

(4)

Wie Abb. 3 zu entnehmen ist, führt die Berücksichtigung der Bindungseffekte zu einer Abnahme des Streuquerschnittes für Streuwinkel q >90°, sowie für hohe Energien und Atome mit kleiner Ordnungszahl Z. Das bedeutet, dass im Gegensatz zur Compton-Streuung die Vorwärtsrichtung mit wachsender Energie wahrscheinlicher wird. Abb. 3 stellt den Effekt des Formfaktors für Aluminium für die gleichen Energien aus Abb. 2 dar. Es sei hier darauf hingewiesen, dass für hohe Energien nur ein extrem kleiner Winkelbereich für die kohärente Streuung relevant ist. Vom Standpunkt der Transporttheorie aus betrachtet bedeutet das, dass im Prinzip kein Streuakt stattgefunden hat.

Abb 3: Vergleich von Thomson- und Rayleigh-Streuquerschnitten für Aluminium bei Photonenenergien von 40 keV (links) und 100 keV (rechts).

Implementation der Monte-Carlo-Simulation

Bei der Erstellung des Modells für den Strahlungstransport wurden verschiedene Annahmen gemacht. Erstens wurde eine geradlinige Bewegung der Photonen zwischen den Wechselwirkungen ohne Energieverlust angenommen. Weiterhin wird genau eine Wechselwirkung pro Einheitslänge, gegeben durch den linearen Schwächungskoeffizienten µ, angenommen. Die Dauer einer Wechselwirkung wird im Vergleich zur Zeit zwischen zwei Wechselwirkungen vernachlässigt. Für den Photonentransport werden drei unabhängige zufällige Ereignisse angenommen: Absorption t, Compton-Streuung sC und Rayleigh-Streuung sR. Die Wahrscheinlichkeit einer Wechselwirkung für jedes der drei Wechselwirkungsereignisse ergibt sich aus dem Verhältnis der Wechselwirkungsquerschnitte eines Ereignisses und der Gesamtschwächung µ:

Der Ort der nächsten Wechselwirkung kann aus dem Schwächungsgesetz abgeleitet werden. Die Wahrscheinlichkeit dafür, dass ein Teilchen einer Wechselwirkung zwischen l und l+dl entlang seiner Bahn unterliegt, ergibt sich zu

(5)


Abb 4: Monte-Carlo-Simulation der Bahn eines einzelnen Photons. Gestrichelte Linien verlängern die Bahnabschnitte zum Detektor, wo ein gewichteter Anteil des Photons zum Zweck der Erhöhung der Statistik registriert wird.

Es folgt für die Wechselwirkungslänge l durch Ziehen einer Zufallszahl x aus [1-P(l)]

(6)
mit Bei der Annahme stückweiser homogener Materialparameter ergibt sich die Wechselwirkungslänge zu

(7)
Das Wiederholungsschema für alle Realisierungen besteht im wesentlichen aus drei Elementen: Bestimmung der Anfangsbedingungen, Bestimmung des Ortes der nächsten Wechselwirkung aus der Wechselwirkungslänge l, Bestimmung des Wechselwirkungsereignisses ei=(t ,sC ,sR ) . Verlässt ein Photon das Durchstrahlungsobjekt, so werden keine weiteren Wechselwirkungen berücksichtigt und die Realisierung ist beendet. Bei Absorption des Photons wird die Realisierung ebenfalls beendet. Im Fall eines Streuereignisses erfolgt die Bestimmung der neuen Bahnrichtung (E,W) → (E¢,W¢) und des möglichen Energieverlusts. Weiterhin erfolgt bei Streuereignissen ein erneuter Durchlauf des Wiederholungsschemas für dieselbe Realisierung ab dem zweiten Schritt.

Zur Varianzverminderung wird bei Streuereignissen ein besonderes Vorgehen eingeführt. Zuerst wird die Teilchenbahn zum Detektor weiterverfolgt und ein gewichteter Anteil des Photons registriert. Anschließend wird die beschriebene Behandlung von Streuereignissen angewandt. Hiermit wird die Statistik der nicht wechselgewirkten sowie der gestreuten Photonen verbessert. Die Verbesserung der Streustrahlungsstatistik wirkt sich insbesondere auf die ersten Streuereignisse einer Realisierung aus, die von höherer Bedeutung sind. Abb. 4 veranschaulicht beispielhaft den vorgestellten Algorithmus.

Bauteil-Repräsentation

Neben dem Modell des Strahlungstransports benötigt die Durchstrahlungssimulation eine virtuelle Bauteil-Repräsentation. Im Gegensatz zu den Standard-Monte-Carlo-Paketen, z.B. MCNP, wurde hier eine direkte Schnittstelle zum CAD realisiert. Geometriedaten werden im STL-Format (Dateiformat für die Stereolithographie, einem Verfahren zur Erzeugung von Werkstückmodellen aus einem Photopolymer) gelesen, einem Dateiformat das von allen 3D-CAD-Programmen unterstützt wird.

Das Geometriemodell der Prüfteile besteht aus ebenen Dreiecken die eine geschlossene Oberfläche bilden und das homogene Material des Bauteils von der Umgebung abgrenzen (s. Abb. 5).


Abb 5: Bauteil-Repräsentation in einem CAD-Programm (links) und facettiertes Oberfläche (Drahtmodell, rechts).

Simulationsrechnungen

Erste Berechnungen wurden zum Vergleich mit dem Standard-Paket MCNP durchgeführt. Dazu wurde die Durchstrahlung einer begrenzten ebenen Platte gewählt. Wie in MCNP ist es auch hier möglich verschiedene physikalische Modelle auszuwählen. Entsprechend wird im Folgenden zwischen der "einfacher" und "voller" Physik unterschieden. Die Komplexität des Modells zu vereinfachen bringt Vorteile bei der Überprüfung des Modells aber auch bei der Anwendung da sich der Rechenaufwand und damit die Rechenzeit reduziert.

Im Modus "einfache Physik" wird die Rayleigh-Streuung vernachlässigt und die Compton-Streuung wird mit der Klein-Nishina-Formel beschrieben, d.h. Bindungseffekte werden vernachlässigt. Im Modus "volle Physik" wird die Klein-Nishina-Formel durch die Streufunktion zur Berücksichtigung von Bindungseffekten der Elektronen erweitert und mit der Thomson-Formel einschließlich der Formfaktoren zur Berücksichtigung der Bindungseffekte wird die Raileigh-Streuung berücksichtigt. Abb. 6 zeigt die Streustrahlenverteilung bei der Durchstrahlung einer 1 cm-Stahlplatte für "einfache" und "volle" Physik. Die Helligkeit in den Bildern ist Maß für die Intensität der Streustrahlung.


Abb 6: Synthetische Streuverteilungen mit (links) und ohne (rechts) Rayleigh-Anteil für monochromatische (100 keV) Durchstrahlung einer 1 cm-Stahlplatte der Größe 5 x 5 cm.

Zwischen "einfacher" und "voller" Physik ist im vorliegenden Beispiel ein deutlicher Unterschied zu erkennen (s. Abb. 6). Der Anteil der Rayleigh-Streuung ist von der Strahlungsenergie abhängig und nimmt mit der Photonenenergie ab. Abb. 7 stellt das Verhältnis der Streustrahlenintensitäten bei "einfacher" und "voller" Physik in Abhängigkeit von der Strahlungsenergie grafisch dar. Diese Streuintensitäten wurden für die Durchstrahlung von 1 cm Aluminium berechnet. Für diesen Fall konnte damit gezeigt werden, dass für Photonenenergien über 200 keV der Modus "einfache Physik" eine gute Näherung darstellt und entsprechende Berechnungen ohne Reduzierung der Genauigkeit beschleunigen kann.


Abb 7: Vergleich der Streuintensitäten für "einfache" und "volle" Physik bei der Durchstrahlung einer 1 cm-Aluminiumplatte und monochromatischer Strahlung.

Zum Vergleich des vorgestellten Monte-Carlo-Programms (KIAM) mit MCNP wird noch einmal die schon bei Abb. 6 zugrunde liegende Durchstrahlung einer 1 cm-Stahlplatte betrachtet. In Abb. 8 werden Intensitätsprofile gezeigt, die mit dem KIAM-Algorithmus bzw. mit MCNP Version 4 berechnet wurden. In der Ansicht der Durchstrahlungsanordnung (Abb. 8a) ist die Profillinie auf der Detektorfläche weiß markiert. Weiterhin ist das Primärstrahlenbündel durch einen transparenten Kegel wiedergegeben.

Die Profile der Primärintensitäten (Abb. 8b) weisen im mittleren Teil, im von der Stahlplatte überdeckten Detektorbereich, eine relative Intensität nahe Null auf, die nach außen an der Kante der Stahlplatte sprunghaft auf einen Wert von 1, die gesamte Primärintensität erreicht den Detektor, ansteigt. Weiter nach außen schließt sich ein sprunghafter Totalabfall der Intensitätswerte an, wo das Primärstrahlenbündel verlassen wird. Die Abweichungen zwischen den beiden Primärintensitätsprofilen sind rein statistischer Natur.

Abb 8: Vergleich des vorgestellten KIAM-Algorithmus mit dem Standard-Monte-Carlo-Code MCNP: Durchstrahlungssimulation für eine 5x5x1 cm-Stahlplatte; a) Durchstrahlungsanordnung mit transparentem Strahlenbündel und Profillinie in der Detektorfläche; b) Primärintensitätsprofile; c) Streustrahlungsprofile

Auch die Profile der Streuintensitäten weisen einen vergleichbaren Verlauf auf. Im Bereich der Stahlplatte wurde ein Abschnitt konstanter mittlerer Streuintensität registriert. Mit zunehmender Entfernung von der Stahlplatte nimmt die Streuintensität ab. Zu beachten ist das deutliche Überschwingen des Intensitätsverlaufs im Kantenbereich, das im Verhältnis von freier Weglänge zu Materialdicke (hier kleiner Eins) begründet ist. Unterschiede zwischen den Profilen beider Algorithmen sind in der Höhe der absoluten Streuintensitätswerte und im Rauschen erkennbar. Der Unterschied der Profilhöhen wird Gegenstand weiterer Untersuchungen sein. Das geringere Rauschen bei gleicher Anzahl von Realisierungen zeigt die deutliche Wirkung der varianzverbessernden Techniken im KIAM-Algorithmus.

Neben der Varianz ist der KIAM-Algorithmus auch in Bezug auf die Rechenzeit erheblich effizienter. Im betrachteten Beispiel wurde für 100 Millionen Realisierungen bei 300 x 300 Detektorelementen mit MCNP eine Rechenzeit von 298,8 h (12,5 h Rechenzeit eines 25-Prozessor-Clusters) benötigt. Die vergleichbare Rechnung mit dem KIAM-Algorithmus ergab dagegen eine Rechenzeit von 22,2 min (5,5 min Rechenzeit eines 4-Prozessor-Clusters).

Ein detaillierter Vergleich des vorgestellten Monte-Carlo-Programms mit MCNP hinsichtlich Genauigkeit und Rechenzeit bei komplexer Geometrie ist in Vorbereitung.

Zusammenfassung

Mit dem KIAM-Algorithmus wurde eine Monte-Carlo-Simulation der Strahlungsausbreitung unter Berücksichtigung der in der Radiologie relevanten Wechselwirkungsmechanismen bei effektiver Strahlverfolgung an facettierten Bauteilrepräsentationen vorgestellt. Durch eine einfache Formulierung virtueller Durchstrahlungsanordnungen und erheblich reduzierten Rechenzeiten im Vergleich zu Standard-Monte-Carlo-Codes wurde eine deutliche Verbesserung der praktischen Anwendbarkeit der Monte-Carlo-Simulation für die industrielle Radiologie erreicht.

Literatur

  1. Briesmeister, J. F. (ed.), MCNP-A General Monte Carlo N-Particle Transport Code, LANL Report LA-13709-M, Los Alamos, 2000
  2. Hubbell, J. H., Veigele, W. J., Briggs, E. A., Brown, R. T., Cromer, D. T., and Howerton, R. J.: Atomic Form Factors, Incoherent Scattering Functions, and Photon Scattering Cross Sections. J. Phys. Chem. Ref. Data 4, 471-538 (1975); erratum in 6, 615-616 (1977)
  3. Hubbell, J. H. and Overbo: Relativistic Atomic Form Factors and Photon Coherent Scattering Cross Sections. J. Phys. Chem. Ref. Data 8, 69-105, (1979)

STARTHerausgeber: DGfZPProgrammierung: NDT.net