Ga naar inhoud

Newton-Raphson-methode

Linearisatie en iteratieve recursie

Virtuele arbeid van externe krachten en assemblage van de globale vergelijkingen geeft een niet-lineaire vergelijking voor de knooppuntverplaatsing op tijdstip \(t_{n+1}\), \(\boldsymbol{u}_{n+1}\), die met de Newton-Raphson-methode wordt opgelost. De knooppuntverplaatsing tot tijdstip \(t_n\), \(\boldsymbol{u}_n\), wordt als bekend verondersteld en het verplaatsingsincrement \(\Delta\boldsymbol{u}\) wordt genomen als de te bepalen onbekende variabele

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

Hierna wordt de afhankelijkheid van de externe-krachtvector van de knooppuntverplaatsing verwaarloosd en wordt, met \(\boldsymbol{F}(\boldsymbol{u}_{n+1}) = \boldsymbol{F}_{n+1}\),

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

opgelost.

Definieer bij de huidige oplossing \(\Delta\boldsymbol{u}\) de tangentiële stijfheid

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

Met behulp hiervan geeft linearisatie van de niet-lineaire vergelijking

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

Laat de correctie in de \(i\)-de iteratie \(d\boldsymbol{u}_i\) zijn en de residuvector aan het begin van de iteratie

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

Dan is de iteratieve recursie

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

Het residu \(\boldsymbol{R}_i\) is dus een grootheid die overeenkomt met de onbalans van krachten ten opzichte van het evenwicht.

Opbouw van de tangentiële stijfheidsmatrix

De tangentiële stijfheid \(\boldsymbol{K} = \partial\boldsymbol{Q}/\partial\boldsymbol{u}\) wordt opgebouwd door de in Discretisatie van interne virtuele arbeid verkregen interne-krachtvector van het element partieel te differentiëren naar de knooppuntverplaatsing, de resulterende integranden op elementniveau over elk elementdomein te integreren en ze te assembleren. Als de integrand op elementniveau wordt aangeduid met \(\boldsymbol{K}^e_X\) (notatie in de referentieconfiguratie, TL-formulering) of \(\boldsymbol{K}^e_x\) (notatie in de huidige configuratie, UL-formulering), is de tangentiële elementstijfheid

\[ \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}) \]

Hieronder staan de uiteindelijke vormen van de TL/UL-integranden. In beide gevallen worden ze ontbonden in de som van een materiaalstijfheidsterm (initiële-verplaatsingsterm) en een geometrische-stijfheidsterm (initiële-spanningsterm).

Total Lagrange-formulering

In de Total Lagrange-formulering wordt een lineair verband aangenomen tussen de snelheid van de tweede Piola-Kirchhoff-spanning (2) \(\dot{\boldsymbol{S}}\) en de Green-Lagrange-reksnelheid \(\dot{\boldsymbol{E}}\), namelijk \(\dot{\boldsymbol{S}} = \boldsymbol{\mathsf{C}}:\dot{\boldsymbol{E}}\). Dit komt overeen met constitutieve wetten voor lineair-elastische materialen (St. Venant-Kirchhoff-materialen) en hyperelastische materialen, en FrontISTR gebruikt voor deze materialen de Total Lagrange-formulering. De integrand van de tangentiële elementstijfheid wordt dan in tensorvorm geschreven 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}}) \]

Term 1 aan de rechterzijde is de materiaalstijfheidsterm (initiële-verplaatsingsterm) en term 2 is de geometrische-stijfheidsterm (initiële-spanningsterm).

In de FrontISTR-implementatie wordt deze integrand in matrixvorm geëvalueerd met behulp van Voigt-notatie:

\[ \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 \]

De matrices zijn als volgt. \(\boldsymbol{B}_L, \boldsymbol{B}_{NL}\) zijn de B-matrices die zijn geïntroduceerd in Discretisatie van interne virtuele arbeid, en \(\tilde{\boldsymbol{C}}\) is de Voigt-weergave van de constitutieve tensor \(\boldsymbol{\mathsf{C}}\), dus een \(6\times 6\) materiaalstijfheidsmatrix (Tensornotatie en wiskundige grondslagen). \(\boldsymbol{S}_9, \boldsymbol{F}_9\) zijn de volgende herschikkingsmatrices die worden gebruikt om de geometrische-stijfheidsterm als matrixproduct uit te drukken. Definieer eerst voor een \(3\times 3\) tensor van orde 2 \(\boldsymbol{A}\) de notatie \([\,\cdot\,]\), die deze herschikt tot een vector met 9 componenten, als

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

Met deze definitie drukt \(\boldsymbol{F}_9\) de variatie van de vervormingsgradiënt uit in de vorm \([\delta\boldsymbol{F}] = \boldsymbol{F}_9\, \delta\boldsymbol{u}^e\) en is het een \(9\times d n_e\)-matrix. Voor elementknooppunt \(\alpha = 1, \ldots, n_e\) is het overeenkomstige \(9\times d\)-blok

\[ [\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{ is de } 3\times 3 \text{ eenheidsmatrix}) \]

en wordt gegeven door \(\boldsymbol{F}_9 = [[\boldsymbol{F}_9]_1, \ldots, [\boldsymbol{F}_9]_{n_e}]\), waarbij de blokken horizontaal in de volgorde van de elementknooppunten zijn gerangschikt. \(\boldsymbol{S}_9\) wordt zo gekozen dat, in combinatie met deze matrix, de geometrische-stijfheidsterm wordt uitgedrukt als \(\delta\boldsymbol{u}^{eT}\, \boldsymbol{F}_9^T \boldsymbol{S}_9 \boldsymbol{F}_9\, \dot{\boldsymbol{u}}^e\); dit is de volgende \(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} \]

Dit is de resulterende matrix.

Updated Lagrange-formulering

In de Updated Lagrange-formulering wordt een lineair verband aangenomen tussen de Jaumann-snelheid van de relatieve Kirchhoff-spanningstensor \(\hat{\boldsymbol{\sigma}}^{\nabla J}\) en de vervormingssnelheidstensor \(\boldsymbol{D}\), namelijk \(\hat{\boldsymbol{\sigma}}^{\nabla J} = \boldsymbol{\mathsf{C}}:\boldsymbol{D}\). Dit is de vorm van een hypo-elastische constitutieve wet die gebruikelijk is voor lineair-elastische, elastoplastische en kruipmaterialen, en FrontISTR gebruikt voor deze materialen de Updated Lagrange-formulering. De integrand van de tangentiële elementstijfheid, uitgedrukt in de huidige configuratie, wordt dan in tensorvorm geschreven 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}) \]

waarbij \(\boldsymbol{\sigma}^{\nabla T}\) de Truesdell-snelheid is, \(\boldsymbol{A}_{(L)}\) het lineaire deel van de Almansi-rek is, \(\boldsymbol{F}_t = \partial\boldsymbol{u}/\partial\boldsymbol{x}\) de verplaatsingsgradiënt ten opzichte van de huidige configuratie is en \(\boldsymbol{L}\) de snelheidsgradiënttensor is. Term 1 aan de rechterzijde is de materiaalstijfheidsterm en term 2 de geometrische-stijfheidsterm.

In de FrontISTR-implementatie wordt deze integrand in matrixvorm geëvalueerd met behulp van Voigt-notatie:

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

Hierbij is \(\boldsymbol{b}\) de B-matrix die in de huidige configuratie is opgebouwd (Discretisatie van interne virtuele arbeid). \(\boldsymbol{\sigma}_9, \boldsymbol{f}_9\) worden verkregen uit \(\boldsymbol{S}_9, \boldsymbol{F}_9\), die voor de TL-formulering zijn gedefinieerd, door de tweede PK-spanning (2) \(\boldsymbol{S}\) te vervangen door de Cauchy-spanning \(\boldsymbol{\sigma}\) en de gradiënt in de referentieconfiguratie \(\partial N_\alpha^e/\partial X_i\) door de gradiënt in de huidige configuratie \(\partial N_\alpha^e/\partial x_i\).

\(\boldsymbol{G}\) is een van de Cauchy-spanning afhankelijke correctiematrix die nodig is om de hypo-elastische constitutieve wet \(\hat{\boldsymbol{\sigma}}^{\nabla J} = \boldsymbol{\mathsf{C}}:\boldsymbol{D}\) consistent te maken met het raamwerk van tangentiële stijfheid als een constitutieve wet op basis van de Truesdell-snelheid. Zij wordt verkregen door de componenten van de tensor van orde 4 \(G_{ijkl} = \delta_{il}\sigma_{kj} + \delta_{jl}\sigma_{ik}\) als volgt in \(6\times 6\)-Voigt-vorm te rangschikken

\[ \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} \]

Dit is de resulterende matrix.

Assemblage van de globale stijfheidsmatrix

De globale tangentiële stijfheid \(\boldsymbol{K}\) wordt verkregen door elke elementstijfheid \(\boldsymbol{K}^e\) voor elk paar knooppunten op te delen in \(d\times d\)-blokken \(\boldsymbol{K}^e_{\alpha\beta}\) en de assemblageverzameling voor tensoren van orde 2 \(\mathcal{E}^2(i_g, i_h)\) te gebruiken die is geïntroduceerd in Assemblage van fysische grootheden op elementknooppunten:

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

De resulterende waarden worden gerangschikt als een matrix met rij \(i_g\) en kolom \(i_h\). In de implementatie wordt de verzameling \(\mathcal{E}^2\) niet expliciet opgebouwd; in plaats daarvan worden de overeenkomstige blokken rechtstreeks binnen de elementlus opgeteld. De matrix is vierkant met een dimensie die gelijk is aan het aantal vrijheidsgraden per knooppunt \(\times\) het totale aantal knooppunten \(n_g\), maar omdat componenten anders dan die tussen via elementen verbonden knooppunten \(0\) zijn, wordt zij opgeslagen in sparse-matrixvorm.

De elementstijfheidsmatrices voor de TL- en UL-formuleringen hebben dezelfde vorm, afgezien van het wisselen van de referentieconfiguratie (knooppuntcoördinaten en de bron waarmee de B-matrix wordt opgebouwd) en de aanwezigheid of afwezigheid van de matrix \(\boldsymbol{G}\). FrontISTR implementeert beide formuleringen daarom in een gemeenschappelijke subroutine.

Iteratie-algoritme

Samengevat: stel aan het begin van de iteratie \(\Delta\boldsymbol{u} = \boldsymbol{0}\) in en bereken het initiële residu \(\boldsymbol{R}_0 = \boldsymbol{F}_{n+1} - \boldsymbol{Q}(\boldsymbol{u}_n)\). Voer vervolgens in de \(i\)-de iteratie de volgende procedure uit.

  1. Bereken bij de huidige verplaatsing \(\boldsymbol{u}_n + \Delta\boldsymbol{u}\) de tangentiële stijfheid \(\boldsymbol{K}_i = \partial\boldsymbol{Q}/\partial\boldsymbol{u}|_{\boldsymbol{u}_n + \Delta\boldsymbol{u}}\) volgens de procedure in Opbouw van de tangentiële stijfheidsmatrix.
  2. Om geometrische randvoorwaarden op te leggen, wijzigt u de tangentiële stijfheidsmatrix en de residuvector voor de vrijheidsgraden waarop verplaatsingsbeperkingen van toepassing zijn, zodat \(\tilde{\boldsymbol{K}}_i, \tilde{\boldsymbol{R}}_{i-1}\) worden verkregen (Behandeling van geometrische randvoorwaarden).
  3. Los de lineaire vergelijking \(\tilde{\boldsymbol{K}}_i\, d\boldsymbol{u}_i = \tilde{\boldsymbol{R}}_{i-1}\) op om de correctie \(d\boldsymbol{u}_i\) te verkrijgen. Deze procedure neemt vaak het grootste deel van de rekenkosten van de iteratieve berekening voor haar rekening.
  4. Werk het verplaatsingsincrement bij als \(\Delta\boldsymbol{u} \leftarrow \Delta\boldsymbol{u} + d\boldsymbol{u}_i\) en bereken overeenkomstig de interne-krachtvector \(\boldsymbol{Q}(\boldsymbol{u}_n + \Delta\boldsymbol{u})\) en het residu \(\boldsymbol{R}_i = \boldsymbol{F}_{n+1} - \boldsymbol{Q}(\boldsymbol{u}_n + \Delta\boldsymbol{u})\).
  5. Controleer de convergentie en beëindig de iteratie als convergentie is bereikt. Componenten die overeenkomen met reactiekrachten op begrensde vrijheidsgraden verschijnen in het residu \(\boldsymbol{R}_i\) voor vrijheidsgraden waarop geometrische randvoorwaarden van toepassing zijn; daarom wordt de convergentie-indicator opgebouwd uit \(\tilde{\boldsymbol{R}}_i\) nadat deze componenten zijn uitgesloten. Specifieke convergentie-indicatoren en drempelwaarden worden beschreven in Convergentiecriteria. Als geen convergentie wordt bereikt en de iteratielimiet is bereikt, wordt de iteratie als mislukt beschouwd.

Wanneer de iteratie convergeert, voegt u de geconvergeerde \(\Delta\boldsymbol{u}\) toe aan \(\boldsymbol{u}_n\) om de geaccumuleerde verplaatsing op tijdstip \(t_{n+1}\), \(\boldsymbol{u}_{n+1} = \boldsymbol{u}_n + \Delta\boldsymbol{u}\), te verkrijgen en gaat u verder naar de volgende tijdstap.

Gerelateerde onderwerpen