Gleichungslösermethoden bei nichtlinearen Berechnungen

Einführung

Die Finite-Elemente-Methode kommt immer dann zum Einsatz, wenn sich mechanische Problemstellungen analytisch nicht lösen lassen. Häufig werden hierbei auch nichtlineare Effekte wie Ausfall bei Druck (geometrische Nichtlinearität), Plastifizierungen (materielle Nichtlinearität) und Kontakt- beziehungsweise kinematische Freiheitsgrade berücksichtigt. Diese Effekte, insbesondere bei nichtlinearen Materialmodellen, lassen sich über eine iterative Berechnungsmethode berücksichtigen.

FEM-Formulierung

Grundlegende Schritte einer FEM-Formulierung (weitere Informationen finden sich in [1]):

1. Schwache Form des Gleichgewichts
   $${\int }_{\mathrm{B}}\underset{\mathrm{nichtlinearer} \mathrm{Anteil}}{\underbrace{{\mathrm{grad\delta }}_{\mathrm{u}}\cdots \mathrm{\sigma dV}}} = {\int }_{\mathrm{B}}{\mathrm{\delta }}_{\mathrm{u}} · \mathrm{\rho bdV} + {\int }_{\partial \mathrm{B}}{\mathrm{\delta }}_{\mathrm{u}} · {\mathrm{t}}_{\mathrm{o}}\mathrm{dA}$$

   δu	Virtuelle (Test-)Verschiebung
   t0	Anfangslastfaktor
   σdV	Innere Kräfte
   ρbdV	Volumenkräfte
   B	integriertes Gebiet

   Schwache Form
2. Umformung in Voigt Notation mit Tensor 4. Stufe
   $${\int }_{\mathrm{B}}\mathrm{\delta \epsilon }\cdots \mathrm{ℂ}\cdots \mathrm{\epsilon dV} = {\int }_{\mathrm{B}}\underset{\mathrm{Volumenkräfte}}{\underbrace{{\mathrm{\delta u}}^{\mathrm{T}}\mathrm{bdV}}} + {\int }_{\partial \mathrm{B}}\underset{\mathrm{Oberflächenkräfte}}{\underbrace{{\mathrm{\delta u}}^{\mathrm{T}} · {\mathrm{t}}_0\mathrm{dA}}}$$

   C	Steifigkeitsmatrix
   δε	Variation des Verzerrungszustands
   B	integriertes Gebiet
   ε	Dehnung
   u	Verformung

   Schwache Form in Voigt Notation
3. Mit dieser Notation wird im Weiteren die Näherungslösung des FE-Ansatzes bei nichtlinearen Materialien gelöst.
4. Das Verschiebungsfeld wird dazu mit den Ansatzfunktionen multipliziert.
   $$\mathrm{u}(\mathrm{x},\mathrm{t}) = \mathrm{H}\left(\mathrm{x}\right)\hat{\mathrm{u}}\left(\mathrm{t}\right)$$

   u(x,t)	Verschiebung über Zeit (Lastinkrement)
   H	Formfunktion
   û	Knotenverschiebung

   Verschiebungsfeld über Ansatzfunktion

Ableitung der Verschiebung in die schwache Form einsetzen. Über numerische Integration werden Knotenverschiebung und im Post Processing die Spannungen und Verzerrungen über das Stoffgesetz berechnet.

Ablauf einer Newton-Raphson-Iteration

Aufgrund des nichtlinearen Materialverhaltens ändert sich die Materialmatrix C in obiger Gleichung 2 mit jedem Dehnschritt. Das Standard-Rechenverfahren zur Lösung dieses Problems ist die sogenannte Newton-Raphson-Iteration. Mit ihr wird die Funktion in einem Ausgangspunkt linearisiert. In der Iteration wird immer die Steifigkeitsmatrix C des Vorschritts verwendet. In dem linearisierten Iterationsschritt wird eine Tangente an die Nullstelle der Funktion gelegt.

Die zum Ablaufschema im obigen Bild gehörenden Gleichungen:

1. Last in Lastschritte zerlegen.
   $${}^{\mathrm{t} + ∆\mathrm{t}}\mathrm{f}_{\mathrm{ext}} = {}^{\mathrm{t}}\mathrm{f}_{\mathrm{ext}} + ∆\mathrm{f}$$

   fext	Äußere Last
   t	Zeit/Lastschritt
   Δf	Neues Lastinkrement

   Last in Lastschritte
2. Prädiktorschritt
   $${}_0{}^{\mathrm{t}}\mathrm{K}{∆}_0\mathrm{\phi } = {}^{\mathrm{t}+∆\mathrm{t}}\mathrm{f}_{\mathrm{ext}} - {}_0{}^{\mathrm{t}}\mathrm{f}_{\mathrm{int}}$$

   K	Steifigkeitsmatrix des vorherigen Zeitschritts
   t+Δtfext	Äußere Kräfte um eine weitere Laststufe erhöht
   0tfint	Innere Kräfte des vorherigen Zeitschritts
   ϕ	Verzerrung

   Prädiktorschritt
3. Im Punkt 3 Iteration des Ablaufschemas wird die Gesamtverzerrung abzüglich der plastischen Verzerrung berechnet (Korrektorschritt). Das Ziel der iterativen Berechnung ist immer, dass die Summe der Lasten Null ergibt. Numerisch ist dies jedoch nicht möglich. Daher wird eine Abbruchrate ε definiert, bei der die Berechnung als ausreichend genau abgebrochen wird.
   $$\mathrm{R} = \left|{}^{\mathrm{t} + ∆\mathrm{t}}\mathrm{f}_{\mathrm{ext}} - {}_{\mathrm{n}}{}^{\mathrm{t} + ∆\mathrm{t}}\mathrm{f}_{\mathrm{int}}\right| < \mathrm{\epsilon }$$

   R	Abbruchrate
   fext	äußere Kräfte
   fint	innere Kräfte
   ε	Epsilon Abbruchrate
   t	Zeitschritt

   Abbruchrate

Im Programm kann die Abbruchrate unter den Berechnungsparametern gesteuert werden.

Im folgenden Bild wird der Ablauf einer Newton-Raphson-Iteration dargestellt. In der ersten Iteration

wird die Abbruchrate R beziehungsweise ε nicht erreicht. Auch in der zweiten Iteration (rot) wird die Toleranzschranke nicht erreicht. Erst in der dritten Iteration ist der Abstand der Tangentensteifigkeit so klein, dass Konvergenz erreicht wird.

Wie bereits erwähnt, wird die Verformung während der Iteration immer weiter aufsummiert.

Zusammenfassung

Die Newton-Raphson-Iteration hat die Konsistenz- beziehungsweise Konvergenzordnung 2. Die Anzahl der "richtigen" Stellen der Iteration verdoppelt sich bei jedem Schritt. Demzufolge konvergiert eine Newton-Raphson-Iteration quadratisch und die Genauigkeit wird mit jeder Iteration größer, wenn das Verfahren konvergiert. Wenn das Verfahren jedoch nicht konvergiert, geht der Fehler ins Unendliche und die Berechnung wird abgebrochen.

Fehlerursachen sind zum Beispiel eine zu starke Steigung der Last Verformungs-Kurve und zu flache Neigung der Kurve im plastischen Bereich. Wenn die Last-Verformungs-Kurve im obigen Bild im zweiten Iterationsschritt einen zu starken Knick aufwiese, würde die Materialtangente und damit die Steifigkeitsmatrix die Steigung des elastischen Bereichs nicht korrekt abbilden. Die Steigung der Nullstelle wäre in diesem Fall für den plastischen Bereich falsch. Dies ist einer der Gründe, weshalb eine Erhöhung der Lastschritte mit einer verbesserten Konvergenz einhergeht.