Rissausbreitung wirkungsvoll mit Z-set modellieren

Einführung

Die Analyse, wie Risse sich in industriellen Bauteilen ausbreiten, ist entscheidend für deren Sicherheit und Zuverlässigkeit. Das gilt besonders für Branchen mit extremen mechanischen Belastungen und hohen Leistungsanforderungen, etwa die Luft- und Raumfahrt^1–4. Hier verlangt es größte Sorgfalt, die Entstehung und Ausbreitung von Defekten präzise zu modellieren, um schwerwiegende Schäden zu verhindern.

Um diese Herausforderung zu meistern, haben Forscher verschiedene numerische und theoretische Ansätze entwickelt. Die erweiterte Finite-Elemente-Methode (XFEM) bildet Rissdiskontinuitäten anhand lokalisierter Anreicherungsfunktionen ab, ohne das Gitter neu zu erstellen⁵. Dieser effektive Ansatz erfasst lokale Singularitäten und die Verzweigung komplexer Risse präzise. Parallel dazu quantifizieren energiebasierte Kriterien, die auf J- und G-Integralen der linearen elastischen Bruchmechanik (LEFM) beruhen, die Energiefreisetzungsrate sowie die Zähigkeit des Materials unter quasistatischen elastischen Bedingungen⁶. Diese Methoden stützen sich vor allem auf die LEFM-Theorie, die sich für spröde Materialien oder kleine Defekte bewährt. Bei elastoplastischen Materialien stoßen sie jedoch schnell an ihre Grenzen.

Tatsächlich führt Plastizität zu erheblichen Energieverlusten und hinterlässt Restspannungen, was die Analyse und Vorhersage der Rissausbreitung erschwert. Besonders anspruchsvoll ist es, die Wachstumsrichtung eines Risses in dissipativen Umgebungen präzise vorherzusagen.

Die G-Theta-Methode bietet hierfür eine effektive Lösung: Sie kombiniert die Berechnung einer für Brüche in dissipativen Medien typischen Energiemenge mit der genauen Vorhersage des Rissverlaufs.

Das numerische Tool Z-cracks, entwickelt von ONERA⁷ und Transvalor und integriert in die Finite-Elemente-Simulationssuite Z-set, liefert einen generischen und effizienten Ansatz zur Analyse von 3D-Rissen. Das Modul berechnet Spannungsintensitätsfaktoren (SIFs) für statische und gemischte Risskonfigurationen und simuliert die Rissausbreitung mithilfe robuster 3D-adaptiver Remeshing-Techniken sowie einer effizienten Multithread-Nachbearbeitung zur SIF-Extraktion.

In diesem Artikel untersuchen wir die Wirksamkeit der in Z-cracks integrierten G-Theta-Methode. Anhand verschiedener Rissbildungsszenarien analysieren wir ihre Konvergenzeigenschaften, insbesondere im elastoplastischen Bereich.

Rissausbreitung mit der G-Theta-Methode modellieren

Das Z-cracks-Verfahren löst Bruchprobleme numerisch in mehreren Schritten. Der schwierigste Schritt besteht darin, während der Rissausbreitung ein passendes Netz um den Riss zu erstellen. Diese Aufgabe ist komplex, da die Spannungs- und Dehnungsfelder nahe der Rissspitze Singularitäten aufweisen. Zunächst fügt man einen Riss mit vorgegebener Geometrie in ein bestehendes, intaktes Netz ein. Danach generiert man die Schnittmenge zwischen diesem Netz und dem Rissnetz. Das Netz dieser Zwischenstruktur wird so angepasst, dass es eine ausreichend hohe Qualität erreicht, um die Spannungsintensitätsfaktoren (SIFs) an der Rissfront präzise zu berechnen. Für diese Neuvernetzung nutzt Z-cracks externe automatische adaptive Werkzeuge zur Neuvernetzung (z. B. MeshGems), die Finite-Elemente-Berechnungen flexibel unterstützen.

Die Vorhersage der Rissausbreitung stützt sich auf die präzise Bewertung bruchrelevanter Größen, insbesondere der Spannungsintensitätsfaktoren und der Energiefreisetzungsrate G. Eine thermomechanische Finite-Elemente-Analyse ermittelt den Spannungszustand innerhalb der Struktur unter einer definierten Belastung. Aus diesem Spannungsfeld berechnet man die SIFs an mehreren Punkten entlang der Rissfront, um den Ausbreitungswinkel, das Vorrückungsinkrement und die neue Position jedes Knotens an der Front zu bestimmen.

Die Rissfront breitet sich nach bestimmten Gesetzen aus (z. B. Paris-Elber^8,9), die es ermöglichen, die geometrische Entwicklung des Risses im Material abzuschätzen. Dafür muss die Topologie des Finite-Elemente-Modells kontinuierlich aktualisiert werden.

Das Z-Risse-Modul automatisiert den Prozess, indem es Bruchgrößen und den Verzweigungswinkel mithilfe der G-Theta-Methode extrahiert. Diese Methode, entwickelt von Destuynder et al.¹⁰, stützt sich auf Domänenintegrale und die Lagrange-Ableitung der potenziellen Energie in Bezug auf ein virtuelles Geschwindigkeitsfeld der Rissausbreitung, das sogenannte θ-Feld. Wir betrachten einen Referenzbereich Ω, der einen Riss mit der Front Γ ⁰ enthält. M sei ein Punkt auf der Rissfront, S die gekrümmte Abszisse entlang der Front. Die Transformation F ^η verschiebt alle Punkte M im Bereich Ω in einen Bereich Ω^η, ausgelöst durch die Rissausbreitung, und lautet:   

,

wobei η eine Konstante ist (η = 0 entspricht der Referenzkonfiguration) und θ(M) das virtuelle Rissverschiebungsfeld angibt.

Abbildung 1: Darstellung eines Risses innerhalb eines volumetrischen Bereichs Ω.

Man nimmt an, dass das virtuelle Rissausdehnungsfeld θ tangential zur Rissfläche an der Rissfront sowie zu den Körpergrenzen verläuft. Das heißt, die Transformationen F ^η verschieben die Position der Rissfront Γ ⁰, ohne die Körpergrenzen zu verändern.

Die Energiefreisetzungsrate G(θ) für eine virtuelle Rissausdehnung θ ergibt sich aus der Lagrange-Ableitung der potenziellen Energie W:

• Bei thermoelastischem Materialverhalten lässt sich ein domäneninvariantes Integral berechnen. Dieser Ansatz ermöglicht die Berechnung der potenziellen Energie, die bei einer infinitesimalen Verschiebung der Rissfront frei würde – vorausgesetzt, das Material verhält sich rein elastisch und steht unter den durch Thermoelastizität verursachten Restspannungen. Ohne Körperkräfte ergibt sich die Gesamtenergiefreisetzungsrate für diese virtuelle Rissausdehnung wie folgt:

wobei ε^th die thermische Dehnung und ε die Gesamtdehnung darstellt.

• Man kann eine Domäneninvariante auch bei elastoplastischem Materialverhalten Dabei berechnet man die potenzielle Energie, die bei einer infinitesimalen Verschiebung der Rissfront frei wird. Dies geschieht unter der Annahme, dass sich das Material rein elastisch verhält, während es einem festen Restspannungsfeld ausgesetzt ist, das durch thermische und viskoplastische Dehnungen entsteht. Die integrierte Größe, G_p definiert sich durch die folgende Gleichung:

wobei ε^ae die unelastische Dehnung bezeichnet.

Die Energiefreisetzungsrate lässt sich ermitteln, indem man entlang der Rissfront Γ ⁰ integriert und diese mit n_p Kontrollpunkten diskretisiert, die wiederum den Formfunktionen N_j(s) zugeordnet sind. Dazu löst man das folgende Gleichungssystem mit n_p Unbekannten:

Validierung der G-Theta-Formulierung 

Die G-Theta-Formulierung muss validiert werden, um sicherzustellen, dass Risse in einem elastoplastischen Kontext zuverlässig modelliert und evaluiert werden und mit der Materialphysik übereinstimmen. Dafür verglich man die Energiefreisetzungsraten G mit Werten, die für rein elastische Materialien ermittelt wurden. Dieser Vergleich zeigt, wie unelastische Restspannungen – thermische und viskoplastische – die Energiefreisetzung in einem realistischen System beeinflussen, in dem irreversible Verformungen eine entscheidende Rolle spielen.

Um den Ansatz glaubwürdig zu stützen, erfolgte eine Validierung mit analytischen und semi-analytischen Referenzlösungen, sowohl zwei- als auch dreidimensional. Diese Validierungen bestätigten, dass die Methode die erwarteten physikalischen Mechanismen in kontrollierten Konfigurationen korrekt reproduziert. Anhand einer numerischen Konvergenzanalyse wurden die Auswirkungen der Netzverfeinerung und der Auflösungsparameter auf die Stabilität und Genauigkeit der Ergebnisse bewertet.

Durchgeführt wurde zudem eine Sensitivitätsstudie, um den G-Wert (siehe Gleichung (2)) an einem zweidimensionalen gekerbten Probekörper zu berechnen, der der einachsig auf Zug belastet wurde (siehe Abbildung 2(a)). Dabei berücksichtigte man den inneren Integrationsradius r_b und den äußeren Integrationsradius r_a, die den Bereich festlegen, über den G im elastischen und elastoplastischen Fall berechnet wird. Die Analyse bestand darin, einen der Radien konstant zu halten und die Konvergenz von G in Abhängigkeit vom anderen zu untersuchen. Setzte man r_b = 0 und variierte r_a​, zeigten die Ergebnisse im rein elastischen Fall, dass G schnell konvergiert, wenn r_a​ gemäß der üblichen Regel auf das Dreifache der minimalen Elementgröße des Netzes festgelegt wird (siehe Abbildung 2(b)). Im elastoplastischen Fall hingegen trat eine Konvergenz von G erst bei r_a​-Werten über 0,2 mm ein. Dies deutet darauf hin, dass der Außenradius r_a​ mindestens das Dreifache der geschätzten plastischen Zonengröße betragen sollte, um zuverlässige Berechnungen zu gewährleisten.

Abbildung 2: (a) Darstellung einer einseitig gekerbten Probe unter einer Zugbelastung von 50 MPa mit dem zugehörigen Netz. (b) Gegenüberstellung der G-Wert-Entwicklung für den elastischen und elastoplastischen Fall in Abhängigkeit vom äußeren Integrationsradius r_a.

Parallel dazu zeigen die Ergebnisse, dass G vom Innenradius r_b beeinflusst wird, während der Außenradius r_a bei ausreichend großem Abstand der Radien keine Rolle spielt. Wird r_a jedoch kleiner als die plastische Zone, wirken sich Schwankungen von r_b auf die Berechnung von G aus.

Diese Analyse präzisiert die optimale Wahl der Radien, besonders im elastoplastischen Fall. Hier reicht die plastische Zone über drei Netzpunkte hinaus, eine Länge, die im elastischen Fall zuvor r_a bestimmte.

Man untersuchte zudem, wie der Gradient nahe der Rissfront die Werte von G beeinflusst. Der viskoplastische Dehnungsgradient ändert sich über einen Bereich, der nur wenige Elemente der Rissfront umfasst, und wirkt sich auf die Werte von G aus – abhängig vom Parameter h_min. Dieser Parameter gibt den Mindestabstand zur Rissfront an, ab dem man den unelastischen Dehnungsgradienten berücksichtigt. Mit h_min lässt sich der Einfluss des Gradienten in der unmittelbaren Nähe der Rissfront ausschließen, wo starke Schwankungen die Bewertung des mechanischen Feldes stören könnten. Die Sensitivitätsstudien führte man auch auf einem unstrukturierten 2D-Netz und unter Bedingungen großer Verformung durch, um die Konvergenz der G-Werte im elastischen und elastoviskoplastischen Fall zu prüfen.

Die anfängliche Validierung wurde auf dreidimensionale Fälle mit planarer Geometrie (1 mm × 1 mm × 0,5 mm) unter einer Spannung von σ₂₂ = 110 MPa erweitert. In diese Geometrie, die bei z = 0 eine Symmetrieebene besitzt, fügte man einen Randriss ein, wie Abbildung 3(a) zeigt.

Anschließend berechnete man das akkumulierte plastische Dehnungsfeld (epcum) sowohl an der Oberfläche der Geometrie (z = 0,5 mm) als auch in der Mittelebene (z = 0). Die Ergebnisse zeigen leichte Unterschiede zwischen den beiden Bereichen (siehe Abbildung 3(b)). Diese Abweichungen in der plastischen Zone zwischen Mittelebene und Oberfläche entstehen durch Schwankungen der Spannungs- und Dehnungsfelder über die Dicke.

Danach analysierten wir die Konvergenz von G an der Oberfläche und im Kern, indem wir den Außenradius r_a variierten und den Innenradius r_b auf null setzten. Die in Abbildung 3(c) dargestellten Ergebnisse zeigen, dass der Unterschied in G zwischen Kern und Oberfläche sowohl im elastischen als auch im plastischen Bereich vernachlässigbar ist und daher keine Rolle spielt.

(a) - (b)

(c)

Abbildung 3: (a) Die untersuchte Geometrie mit durchgehendem Randriss und das dazugehörige Netz. (b) Das akkumulierte plastische Dehnungsfeld (epcum), das die plastische Zone an der Oberfläche und im Kern der Geometrie darstellt. (c) Ergebnisse der Konvergenzanalyse von G in Abhängigkeit vom Außenradius r_a, getrennt für elastische und plastische Materialien, an der Oberfläche und im Kern.

Bei einem durchgehenden Riss zeigt G im plastischen Fall unterschiedliche Konvergenzen an der Oberfläche und im Kern der Platte, wie Abbildung 4 nachfolgend verdeutlicht. Im Kern schreitet die Konvergenz aufgrund von Schwankungen in der plastischen Zone über die Dicke schneller voran als an der Oberfläche. Im elastischen Fall hingegen konvergiert G an Oberfläche und Kern gleich schnell, obwohl die Grenzwerte variieren – ein Verhalten, das im elastoplastischen Fall nicht auftritt.

Abbildung 4: Die Konvergenzanalyse von G an der Oberfläche und im Kern zeigt die Abhängigkeit vom Außenradius r_a (wenn r_b = 0 ist) für einen durchgehenden Riss, wobei sowohl das plastische als auch das elastische Materialverhalten berücksichtigt werden.

Abbildung 5: Die Konvergenzanalyse von G im Kern (z = 0) zeigt die Abhängigkeit vom Außenradius r_a (unter der Bedingung r_b = 0, h_min  = 0) für zwei wechselwirkende Rissfronten. Untersucht wurden dabei      (a) elastisches und (b) plastisches Verhalten.

Die Analyse zeigt, dass die zweite Rissfront die Konvergenz von G für die erste Rissfront beeinflusst. Tatsächlich sinken im ersten Fall die G-Werte deutlich, sobald der externe Integrationsradius der zweiten Rissfront (r_a = 0,5 mm) näherkommt. Dies stört die Konvergenz der Ergebnisse (siehe blaue Kurven in Abbildung 5).

Im zweiten Fall beeinflusst die Nähe zur geometrischen Grenze (r_a = 0,25 mm) die Ergebnisse (siehe orangefarbene Kurven in Abbildung 5), wodurch die G-Werte jenseits dieser Grenze abfallen. Obwohl beide Konfigurationen mechanisch gleichwertig sind, erklärt die verkleinerte Integrationszone jenseits von r_a  = 0,25 mm diese Abweichung.

Fazit

Diese Studie unterstreicht die Bedeutung der G-Theta-Methode, um die Rissausbreitung in Materialien mit elastoplastischem Verhalten vorherzusagen. Eingebettet in das Z-cracks-Modul der Z-set-suite berechnet dieser Ansatz die beim Risswachstum freigesetzte Energie – selbst in komplexen industriellen Anwendungen.

Dank adaptivem 3D-Remeshing und optimierter Nachbearbeitung bietet Z-cracks eine hocheffiziente Lösung für die dreidimensionale Rissanalyse. Z-set erweist sich damit als leistungsstarke Plattform, die den aktuellen und künftigen Anforderungen der Bruchmechanik gerecht wird.