Newton-Raphson-Verfahren¶

Linearisierung und iterative Rekursion¶

Aus Virtuelle Arbeit äußerer Kräfte und Assemblierung der globalen Gleichung ergibt sich eine nichtlineare Gleichung für die Knotenverschiebung \(\boldsymbol{u}_{n+1}\) zum Zeitpunkt \(t_{n+1}\), die mit dem Newton-Raphson-Verfahren gelöst wird. Die Knotenverschiebung \(\boldsymbol{u}_n\) bis zum Zeitpunkt \(t_n\) wird als bekannt vorausgesetzt, und das Verschiebungsinkrement \(\Delta\boldsymbol{u}\) wird als unbekannte Variable angesetzt, um

\[ \boldsymbol{u}_{n+1} = \boldsymbol{u}_n + \Delta\boldsymbol{u} \]

zu bestimmen. Im Folgenden wird die Abhängigkeit des äußeren Kraftvektors von der Knotenverschiebung vernachlässigt, und mit \(\boldsymbol{F}(\boldsymbol{u}_{n+1}) = \boldsymbol{F}_{n+1}\) wird

\[ \boldsymbol{Q}(\boldsymbol{u}_n + \Delta\boldsymbol{u}) - \boldsymbol{F}_{n+1} = \boldsymbol{0} \]

gelöst.

Für die aktuelle Lösung \(\Delta\boldsymbol{u}\) wird die Tangentialsteifigkeit

\[ \boldsymbol{K} = \left. \frac{\partial \boldsymbol{Q}}{\partial \boldsymbol{u}} \right|_{\boldsymbol{u}_n + \Delta\boldsymbol{u}} \]

definiert. Damit ergibt die Linearisierung der nichtlinearen Gleichung

\[ \boldsymbol{K}\, d\boldsymbol{u} + \boldsymbol{Q}(\boldsymbol{u}_n + \Delta\boldsymbol{u}) - \boldsymbol{F}_{n+1} = \boldsymbol{0} \]

Sei die Korrektur in der \(i\)-ten Iteration \(d\boldsymbol{u}_i\) und der Residualvektor zu Beginn der Iteration

\[ \boldsymbol{R}_{i-1} = \boldsymbol{F}_{n+1} - \boldsymbol{Q}(\boldsymbol{u}_n + \Delta\boldsymbol{u}) \]

Dann lautet die iterative Rekursion

\[ \boldsymbol{K}_i\, d\boldsymbol{u}_i = \boldsymbol{R}_{i-1}, \qquad \Delta\boldsymbol{u} \leftarrow \Delta\boldsymbol{u} + d\boldsymbol{u}_i \]

Das Residuum \(\boldsymbol{R}_i\) entspricht somit der Kraftunausgeglichenheit gegenüber dem Gleichgewicht.

Aufbau der Tangentialsteifigkeitsmatrix¶

Die Tangentialsteifigkeit \(\boldsymbol{K} = \partial\boldsymbol{Q}/\partial\boldsymbol{u}\) wird aufgebaut, indem der in Diskretisierung der virtuellen Arbeit innerer Kräfte erhaltene Element-Innenkraftvektor nach der Knotenverschiebung partiell abgeleitet, die daraus resultierenden Integranden auf Elementebene über die jeweiligen Elementgebiete integriert und anschließend assembliert werden. Bezeichnen wir den Integranden auf Elementebene mit \(\boldsymbol{K}^e_X\) (Notation in der Referenzkonfiguration, TL-Formulierung) bzw. \(\boldsymbol{K}^e_x\) (Notation in der aktuellen Konfiguration, UL-Formulierung), so lautet die Element-Tangentialsteifigkeit

\[ \boldsymbol{K}^e = \int_{\Omega^e_0} \boldsymbol{K}^e_X\, dV \quad (\text{TL}), \qquad \boldsymbol{K}^e = \int_{\Omega^e} \boldsymbol{K}^e_x\, dv \quad (\text{UL}) \]

Im Folgenden werden die endgültigen Formen der TL/UL-Integranden angegeben. In beiden Fällen zerfallen sie in die Summe aus einem Materialsteifigkeitsterm (Anfangsverschiebungsterm) und einem geometrischen Steifigkeitsterm (Anfangsspannungsterm).

Total-Lagrange-Formulierung¶

In der Total-Lagrange-Formulierung wird eine lineare Beziehung zwischen der Rate der zweiten Piola-Kirchhoff-Spannung \(\dot{\boldsymbol{S}}\) und der Green-Lagrange-Dehnungsrate \(\dot{\boldsymbol{E}}\) angenommen, nämlich \(\dot{\boldsymbol{S}} = \boldsymbol{\mathsf{C}}:\dot{\boldsymbol{E}}\). Dies entspricht Stoffgesetzen für linear-elastische Materialien (St.-Venant-Kirchhoff-Materialien) und hyperelastische Materialien; FrontISTR verwendet für diese Materialien die Total-Lagrange-Formulierung. Der Integrand der Element-Tangentialsteifigkeit wird dann in Tensorform geschrieben als

\[ \delta\boldsymbol{u}^{eT}\, \boldsymbol{K}^e_X\, \dot{\boldsymbol{u}}^e = \dot{\boldsymbol{S}}:\delta\boldsymbol{E} + \boldsymbol{S}:(\delta\boldsymbol{F}^T \dot{\boldsymbol{F}}) \]

Der erste Term auf der rechten Seite ist der Materialsteifigkeitsterm (Anfangsverschiebungsterm), der zweite der geometrische Steifigkeitsterm (Anfangsspannungsterm).

In der FrontISTR-Implementierung wird dieser Integrand in Matrixform unter Verwendung der Voigt-Notation ausgewertet:

\[ \boldsymbol{K}^e_X = (\boldsymbol{B}_L + \boldsymbol{B}_{NL})^T\, \tilde{\boldsymbol{C}}\, (\boldsymbol{B}_L + \boldsymbol{B}_{NL}) + \boldsymbol{F}_9^T\, \boldsymbol{S}_9\, \boldsymbol{F}_9 \]

Die einzelnen Matrizen sind wie folgt definiert. \(\boldsymbol{B}_L, \boldsymbol{B}_{NL}\) sind die in Diskretisierung der virtuellen Arbeit innerer Kräfte eingeführten B-Matrizen, und \(\tilde{\boldsymbol{C}}\) ist die Voigt-Darstellung des konstitutiven Tensors \(\boldsymbol{\mathsf{C}}\), also eine \(6\times 6\)-Materialsteifigkeitsmatrix (Tensor-Notation und mathematische Grundlagen). \(\boldsymbol{S}_9, \boldsymbol{F}_9\) sind die folgenden Umordnungsmatrizen, mit denen der geometrische Steifigkeitsterm als Matrixprodukt dargestellt wird. Für einen \(3\times 3\)-Tensor zweiter Ordnung \(\boldsymbol{A}\) definieren wir zunächst die Notation \([\,\cdot\,]\), die ihn in einen Vektor mit 9 Komponenten umordnet:

\[ [\boldsymbol{A}] = (A_{11}, A_{21}, A_{31}, A_{12}, A_{22}, A_{32}, A_{13}, A_{23}, A_{33})^T \]

Mit dieser Definition stellt \(\boldsymbol{F}_9\) die Variation des Deformationsgradienten in der Form \([\delta\boldsymbol{F}] = \boldsymbol{F}_9\, \delta\boldsymbol{u}^e\) dar und ist eine \(9\times d n_e\)-Matrix. Für den Elementknoten \(\alpha = 1, \ldots, n_e\) lautet der zugehörige \(9\times d\)-Block

\[ [\boldsymbol{F}_9]_\alpha = \begin{bmatrix} (\partial N_\alpha^e/\partial X_1)\, \boldsymbol{I} \\ (\partial N_\alpha^e/\partial X_2)\, \boldsymbol{I} \\ (\partial N_\alpha^e/\partial X_3)\, \boldsymbol{I} \end{bmatrix} \qquad (\boldsymbol{I} \text{ ist die } 3\times 3 \text{ Einheitsmatrix}). \]

Damit ist \(\boldsymbol{F}_9 = [[\boldsymbol{F}_9]_1, \ldots, [\boldsymbol{F}_9]_{n_e}]\), wobei die Blöcke in der Reihenfolge der Elementknoten horizontal angeordnet sind. \(\boldsymbol{S}_9\) wird so gewählt, dass in Kombination mit dieser Matrix der geometrische Steifigkeitsterm als \(\delta\boldsymbol{u}^{eT}\, \boldsymbol{F}_9^T \boldsymbol{S}_9 \boldsymbol{F}_9\, \dot{\boldsymbol{u}}^e\) ausgedrückt wird; es handelt sich um die folgende \(9\times 9\)-Matrix

\[ \boldsymbol{S}_9 = \begin{bmatrix} S_{11} \boldsymbol{I} & S_{12} \boldsymbol{I} & S_{13} \boldsymbol{I} \\ S_{21} \boldsymbol{I} & S_{22} \boldsymbol{I} & S_{23} \boldsymbol{I} \\ S_{31} \boldsymbol{I} & S_{32} \boldsymbol{I} & S_{33} \boldsymbol{I} \end{bmatrix} \]

Dies ist die resultierende Matrix.

Updated-Lagrange-Formulierung¶

In der Updated-Lagrange-Formulierung wird eine lineare Beziehung zwischen der Jaumann-Rate des relativen Kirchhoff-Spannungstensors \(\hat{\boldsymbol{\sigma}}^{\nabla J}\) und dem Deformationsgeschwindigkeitstensor \(\boldsymbol{D}\) angenommen, nämlich \(\hat{\boldsymbol{\sigma}}^{\nabla J} = \boldsymbol{\mathsf{C}}:\boldsymbol{D}\). Dies ist die Form eines hypoelastischen Stoffgesetzes, die linear-elastischen, elastoplastischen und Kriechmaterialien gemeinsam ist; FrontISTR verwendet für diese Materialien die Updated-Lagrange-Formulierung. Der in der aktuellen Konfiguration dargestellte Integrand der Element-Tangentialsteifigkeit wird dann in Tensorform geschrieben als

\[ \delta\boldsymbol{u}^{eT}\, \boldsymbol{K}^e_x\, \dot{\boldsymbol{u}}^e = \boldsymbol{\sigma}^{\nabla T}:\delta\boldsymbol{A}_{(L)} + \boldsymbol{\sigma}:(\delta\boldsymbol{F}_t^T\, \boldsymbol{L}) \]

Dabei ist \(\boldsymbol{\sigma}^{\nabla T}\) die Truesdell-Rate, \(\boldsymbol{A}_{(L)}\) der lineare Anteil der Almansi-Dehnung, \(\boldsymbol{F}_t = \partial\boldsymbol{u}/\partial\boldsymbol{x}\) der Verschiebungsgradient bezüglich der aktuellen Konfiguration und \(\boldsymbol{L}\) der Geschwindigkeitsgradiententensor. Der erste Term auf der rechten Seite ist der Materialsteifigkeitsterm, der zweite der geometrische Steifigkeitsterm.

In der FrontISTR-Implementierung wird dieser Integrand in Matrixform unter Verwendung der Voigt-Notation ausgewertet:

\[ \boldsymbol{K}^e_x = \boldsymbol{b}^T\, (\tilde{\boldsymbol{C}} - \boldsymbol{G})\, \boldsymbol{b} + \boldsymbol{f}_9^T\, \boldsymbol{\sigma}_9\, \boldsymbol{f}_9 \]

Hier ist \(\boldsymbol{b}\) die in der aktuellen Konfiguration aufgebaute B-Matrix (Diskretisierung der virtuellen Arbeit innerer Kräfte). \(\boldsymbol{\sigma}_9, \boldsymbol{f}_9\) erhält man aus den für die TL-Formulierung definierten \(\boldsymbol{S}_9, \boldsymbol{F}_9\), indem die zweite PK-Spannung \(\boldsymbol{S}\) durch die Cauchy-Spannung \(\boldsymbol{\sigma}\) und der Gradient in der Referenzkonfiguration \(\partial N_\alpha^e/\partial X_i\) durch den Gradient in der aktuellen Konfiguration \(\partial N_\alpha^e/\partial x_i\) ersetzt werden.

\(\boldsymbol{G}\) ist eine von der Cauchy-Spannung abhängige Korrekturmatrix, die erforderlich ist, um das hypoelastische Stoffgesetz \(\hat{\boldsymbol{\sigma}}^{\nabla J} = \boldsymbol{\mathsf{C}}:\boldsymbol{D}\) mit dem Tangentialsteifigkeitsrahmen als auf der Truesdell-Rate basierendes Stoffgesetz in Einklang zu bringen. Sie entsteht, indem die Komponenten des Tensors vierter Ordnung \(G_{ijkl} = \delta_{il}\sigma_{kj} + \delta_{jl}\sigma_{ik}\) in einer \(6\times 6\)-Voigt-Darstellung angeordnet werden:

\[ \boldsymbol{G} = \begin{bmatrix} 2\sigma_{11} & 0 & 0 & \sigma_{12} & 0 & \sigma_{31} \\ 0 & 2\sigma_{22} & 0 & \sigma_{12} & \sigma_{23} & 0 \\ 0 & 0 & 2\sigma_{33} & 0 & \sigma_{23} & \sigma_{31} \\ \sigma_{12} & \sigma_{12} & 0 & \tfrac{\sigma_{11}+\sigma_{22}}{2} & \tfrac{\sigma_{12}}{2} & \tfrac{\sigma_{23}}{2} \\ 0 & \sigma_{23} & \sigma_{23} & \tfrac{\sigma_{12}}{2} & \tfrac{\sigma_{22}+\sigma_{33}}{2} & \tfrac{\sigma_{12}}{2} \\ \sigma_{31} & 0 & \sigma_{31} & \tfrac{\sigma_{23}}{2} & \tfrac{\sigma_{12}}{2} & \tfrac{\sigma_{33}+\sigma_{11}}{2} \end{bmatrix} \]

Dies ist die resultierende Matrix.

Assemblierung der globalen Steifigkeitsmatrix¶

Die globale Tangentialsteifigkeit \(\boldsymbol{K}\) wird erhalten, indem jede Elementsteifigkeit \(\boldsymbol{K}^e\) für jedes Knotenpaar in \(d\times d\)-Blöcke \(\boldsymbol{K}^e_{\alpha\beta}\) zerlegt und die in Assemblierung physikalischer Größen an Elementknoten eingeführte Assemblierungsmenge zweiter Ordnung \(\mathcal{E}^2(i_g, i_h)\) verwendet wird:

\[ \boldsymbol{K}_{i_gi_h} = \sum_{(e,\alpha,\beta) \in \mathcal{E}^2(i_g, i_h)} \boldsymbol{K}^e_{\alpha\beta} \]

Die resultierenden Werte werden als Matrix mit Zeile \(i_g\) und Spalte \(i_h\) angeordnet. In der Implementierung wird die Menge \(\mathcal{E}^2\) nicht explizit erzeugt; stattdessen werden die entsprechenden Blöcke innerhalb der Elementschleife direkt addiert. Die Matrix ist quadratisch und besitzt die Dimension Freiheitsgrade pro Knoten \(\times\) Gesamtzahl der Knoten \(n_g\). Da jedoch alle Komponenten außer denen zwischen über Elemente verbundenen Knoten \(0\) sind, wird sie in Form einer dünnbesetzten Matrix gespeichert.

Die Elementsteifigkeitsmatrizen der TL- und UL-Formulierung haben dieselbe Form, abgesehen vom Wechsel der Referenzkonfiguration (Knotenkoordinaten und Quelle für den Aufbau der B-Matrix) und dem Vorhandensein bzw. Fehlen der Matrix \(\boldsymbol{G}\). FrontISTR implementiert daher beide Formulierungen in einer gemeinsamen Subroutine.

Iterationsalgorithmus¶

Zusammenfassend wird zu Beginn der Iteration \(\Delta\boldsymbol{u} = \boldsymbol{0}\) gesetzt und das Anfangsresiduum \(\boldsymbol{R}_0 = \boldsymbol{F}_{n+1} - \boldsymbol{Q}(\boldsymbol{u}_n)\) berechnet. Anschließend wird in der \(i\)-ten Iteration das folgende Verfahren durchgeführt.

Berechnen Sie bei der aktuellen Verschiebung \(\boldsymbol{u}_n + \Delta\boldsymbol{u}\) die Tangentialsteifigkeit \(\boldsymbol{K}_i = \partial\boldsymbol{Q}/\partial\boldsymbol{u}|_{\boldsymbol{u}_n + \Delta\boldsymbol{u}}\) nach dem Verfahren unter Aufbau der Tangentialsteifigkeitsmatrix.
Modifizieren Sie zur Durchsetzung geometrischer Randbedingungen die Tangentialsteifigkeitsmatrix und den Residualvektor für die Freiheitsgrade mit Verschiebungszwangsbedingungen, sodass \(\tilde{\boldsymbol{K}}_i, \tilde{\boldsymbol{R}}_{i-1}\) erhalten werden (Behandlung geometrischer Randbedingungen).
Lösen Sie die lineare Gleichung \(\tilde{\boldsymbol{K}}_i\, d\boldsymbol{u}_i = \tilde{\boldsymbol{R}}_{i-1}\), um die Korrektur \(d\boldsymbol{u}_i\) zu bestimmen. Dieser Schritt macht häufig den größten Teil des Rechenaufwands der Iteration aus.
Aktualisieren Sie das Verschiebungsinkrement gemäß \(\Delta\boldsymbol{u} \leftarrow \Delta\boldsymbol{u} + d\boldsymbol{u}_i\) und berechnen Sie entsprechend den Innenkraftvektor \(\boldsymbol{Q}(\boldsymbol{u}_n + \Delta\boldsymbol{u})\) sowie das Residuum \(\boldsymbol{R}_i = \boldsymbol{F}_{n+1} - \boldsymbol{Q}(\boldsymbol{u}_n + \Delta\boldsymbol{u})\).
Prüfen Sie die Konvergenz und beenden Sie die Iteration, wenn Konvergenz erreicht ist. In den Freiheitsgraden mit geometrischen Randbedingungen erscheinen im Residuum \(\boldsymbol{R}_i\) Komponenten, die den Zwangsreaktionen entsprechen; daher wird der Konvergenzindikator aus \(\tilde{\boldsymbol{R}}_i\) nach Ausschluss dieser Komponenten gebildet. Konkrete Konvergenzindikatoren und Schwellenwerte sind unter Konvergenzkriterien beschrieben. Wird keine Konvergenz erreicht und das Iterationslimit erreicht, gilt die Iteration als fehlgeschlagen.

Wenn die Iteration konvergiert, wird das konvergierte \(\Delta\boldsymbol{u}\) zu \(\boldsymbol{u}_n\) addiert, um die akkumulierte Verschiebung zum Zeitpunkt \(t_{n+1}\), \(\boldsymbol{u}_{n+1} = \boldsymbol{u}_n + \Delta\boldsymbol{u}\), zu erhalten, und anschließend zum nächsten Zeitschritt übergegangen.