Slaan oor na inhoud

Newton-Raphson-metode

Linearisering en iteratiewe rekurrensie

Virtuele arbeid van eksterne kragte en samestelling van die globale vergelyking gee ’n nie-lineêre vergelyking vir die nodusverplasing op tyd \(t_{n+1}\), \(\boldsymbol{u}_{n+1}\), wat met die Newton-Raphson-metode opgelos word. Die nodusverplasing tot tyd \(t_n\), \(\boldsymbol{u}_n\), word as bekend aanvaar, en die verplasingstoename \(\Delta\boldsymbol{u}\) word as die onbekende veranderlike geneem om te bepaal

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

Hierna word die afhanklikheid van die eksternekragvektor van nodusverplasing verwaarloos, en 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} \]

word opgelos.

Definieer by die huidige oplossing \(\Delta\boldsymbol{u}\) die tangensiële styfheid

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

Met behulp hiervan gee linearisering van die nie-lineêre vergelyking

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

Laat die korreksie in die \(i\)-de iterasie \(d\boldsymbol{u}_i\) wees, en laat die residuvektor aan die begin van die iterasie wees

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

Dan is die iteratiewe rekurrensie

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

Dus is die residu \(\boldsymbol{R}_i\) ’n grootheid wat ooreenstem met die kragwanbalans vanaf ewewig.

Konstruksie van die tangensiële styfheidsmatriks

Die tangensiële styfheid \(\boldsymbol{K} = \partial\boldsymbol{Q}/\partial\boldsymbol{u}\) word gekonstrueer deur die element-internekragvektor wat in Diskretisering van die virtuele arbeid van interne kragte verkry is gedeeltelik ten opsigte van nodusverplasing te differensieer, die resulterende elementvlak-integrande oor elke elementdomein te integreer en hulle saam te stel. As die elementvlak-integrand aangedui word deur \(\boldsymbol{K}^e_X\) (verwysingskonfigurasienotasie, TL-formulering) of \(\boldsymbol{K}^e_x\) (huidige-konfigurasienotasie, UL-formulering), is die element se tangensiële styfheid

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

Die volgende gee die finale vorme van die TL/UL-integrande. In albei gevalle word hulle ontbind in die som van ’n materiaalstyfheidsterm (aanvanklike-verplasingsterm) en ’n geometriese styfheidsterm (aanvanklike-spanningsterm).

Totale Lagrange-formulering

In die Totale Lagrange-formulering word ’n lineêre verband aanvaar tussen die tempo van die tweede Piola-Kirchhoff-spanning \(\dot{\boldsymbol{S}}\) en die Green-Lagrange-rektempo \(\dot{\boldsymbol{E}}\), naamlik \(\dot{\boldsymbol{S}} = \boldsymbol{\mathsf{C}}:\dot{\boldsymbol{E}}\). Dit stem ooreen met konstitutiewe wette vir lineêr-elastiese materiale (St. Venant-Kirchhoff-materiale) en hiperelastiese materiale, en FrontISTR gebruik die Totale Lagrange-formulering vir hierdie materiale. Die element se tangensiële-styfheidsintegrand word dan in tensorvorm geskryf as

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

Die eerste term aan die regterkant is die materiaalstyfheidsterm (aanvanklike-verplasingsterm), en die tweede term is die geometriese styfheidsterm (aanvanklike-spanningsterm).

In die FrontISTR-implementering word hierdie integrand in matriksvorm met Voigt-notasie geëvalueer:

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

Elke matriks is soos volg. \(\boldsymbol{B}_L, \boldsymbol{B}_{NL}\) is die B-matrikse wat in Diskretisering van die virtuele arbeid van interne kragte bekendgestel is, en \(\tilde{\boldsymbol{C}}\) is die Voigt-voorstelling van die konstitutiewe tensor \(\boldsymbol{\mathsf{C}}\), naamlik ’n \(6\times 6\)-materiaalstyfheidsmatriks (Tensornotasie en wiskundige grondbeginsels). \(\boldsymbol{S}_9, \boldsymbol{F}_9\) is die volgende herrangskikkingsmatrikse wat gebruik word om die geometriese styfheidsterm as ’n matriksproduk uit te druk. Definieer eers, vir ’n \(3\times 3\) tweede-orde-tensor \(\boldsymbol{A}\), die notasie \([\,\cdot\,]\) wat dit in ’n 9-komponentvektor herrangskik as

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

Met hierdie definisie druk \(\boldsymbol{F}_9\) die variasie van die vervormingsgradiënt uit in die vorm \([\delta\boldsymbol{F}] = \boldsymbol{F}_9\, \delta\boldsymbol{u}^e\) en is dit ’n \(9\times d n_e\)-matriks. Vir elementnodus \(\alpha = 1, \ldots, n_e\) is die ooreenstemmende \(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 die } 3\times 3 \text{ identiteitsmatriks}) \]

en word gegee deur \(\boldsymbol{F}_9 = [[\boldsymbol{F}_9]_1, \ldots, [\boldsymbol{F}_9]_{n_e}]\), met die blokke horisontaal in elementnodusvolgorde gerangskik. \(\boldsymbol{S}_9\) word so gekies dat, wanneer dit met hierdie matriks gekombineer word, die geometriese styfheidsterm uitgedruk word as \(\delta\boldsymbol{u}^{eT}\, \boldsymbol{F}_9^T \boldsymbol{S}_9 \boldsymbol{F}_9\, \dot{\boldsymbol{u}}^e\); dit is die volgende \(9\times 9\)-matriks

\[ \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 die resulterende matriks.

Bygewerkte Lagrange-formulering

In die Bygewerkte Lagrange-formulering word ’n lineêre verband aanvaar tussen die Jaumann-tempo van die relatiewe Kirchhoff-spanningstensor \(\hat{\boldsymbol{\sigma}}^{\nabla J}\) en die vervormingstempotensor \(\boldsymbol{D}\), naamlik \(\hat{\boldsymbol{\sigma}}^{\nabla J} = \boldsymbol{\mathsf{C}}:\boldsymbol{D}\). Dit is die vorm van ’n hipoelastiese konstitutiewe wet wat gemeenskaplik is aan lineêr-elastiese, elastoplastiese en kruipmateriale, en FrontISTR gebruik die Bygewerkte Lagrange-formulering vir hierdie materiale. Die element se tangensiële-styfheidsintegrand wat in die huidige konfigurasie uitgedruk word, word dan in tensorvorm geskryf as

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

waar \(\boldsymbol{\sigma}^{\nabla T}\) die Truesdell-tempo is, \(\boldsymbol{A}_{(L)}\) die lineêre deel van die Almansi-rek is, \(\boldsymbol{F}_t = \partial\boldsymbol{u}/\partial\boldsymbol{x}\) die verplasingsgradiënt ten opsigte van die huidige konfigurasie is, en \(\boldsymbol{L}\) die snelheidsgradiënttensor is. Die eerste term aan die regterkant is die materiaalstyfheidsterm, en die tweede term is die geometriese styfheidsterm.

In die FrontISTR-implementering word hierdie integrand in matriksvorm met Voigt-notasie geëvalueer:

\[ \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 is \(\boldsymbol{b}\) die B-matriks wat in die huidige konfigurasie gekonstrueer is (Diskretisering van die virtuele arbeid van interne kragte). \(\boldsymbol{\sigma}_9, \boldsymbol{f}_9\) word verkry uit \(\boldsymbol{S}_9, \boldsymbol{F}_9\) wat vir die TL-formulering gedefinieer is deur die tweede PK-spanning \(\boldsymbol{S}\) met die Cauchy-spanning \(\boldsymbol{\sigma}\) te vervang en die verwysingskonfigurasiegradiënt \(\partial N_\alpha^e/\partial X_i\) met die huidige-konfigurasiegradiënt \(\partial N_\alpha^e/\partial x_i\).

\(\boldsymbol{G}\) is ’n Cauchy-spanning-afhanklike korreksiematriks wat nodig is om die hipoelastiese konstitutiewe wet \(\hat{\boldsymbol{\sigma}}^{\nabla J} = \boldsymbol{\mathsf{C}}:\boldsymbol{D}\) met die tangensiële-styfheidsraamwerk as ’n Truesdell-tempo-gebaseerde konstitutiewe wet konsekwent te maak. Dit word verkry deur die vierde-orde-tensorkomponente \(G_{ijkl} = \delta_{il}\sigma_{kj} + \delta_{jl}\sigma_{ik}\) in \(6\times 6\) Voigt-vorm te rangskik as

\[ \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 die resulterende matriks.

Samestelling van die globale styfheidsmatriks

Die globale tangensiële styfheid \(\boldsymbol{K}\) word verkry deur elke elementstyfheid \(\boldsymbol{K}^e\) in \(d\times d\)-blokke \(\boldsymbol{K}^e_{\alpha\beta}\) vir elke paar nodusse te verdeel en die tweede-orde-tensor-samestellingstel \(\mathcal{E}^2(i_g, i_h)\) te gebruik wat in Samestelling van element-nodusfisiese groothede bekendgestel is:

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

Die resulterende waardes word as ’n matriks met ry \(i_g\) en kolom \(i_h\) gerangskik. In die implementering word die stel \(\mathcal{E}^2\) nie eksplisiet gekonstrueer nie; die ooreenstemmende blokke word eerder direk binne die elementlus bygevoeg. Die matriks is vierkantig met dimensie gelyk aan vryheidsgrade per nodus \(\times\) totale aantal nodusse \(n_g\), maar omdat komponente anders as dié tussen nodusse wat deur elemente verbind is \(0\) is, word dit in ylmatriksvorm gestoor.

Die elementstyfheidsmatrikse vir die TL- en UL-formulerings het dieselfde vorm, behalwe vir die omskakeling van die verwysingskonfigurasie (noduskoördinate en die bron wat gebruik word om die B-matriks te konstrueer) en die teenwoordigheid of afwesigheid van die \(\boldsymbol{G}\)-matriks. FrontISTR implementeer daarom albei formulerings in ’n gemeenskaplike subroetine.

Iterasie-algoritme

Opsommend, stel aan die begin van die iterasie \(\Delta\boldsymbol{u} = \boldsymbol{0}\) en bereken die aanvanklike residu \(\boldsymbol{R}_0 = \boldsymbol{F}_{n+1} - \boldsymbol{Q}(\boldsymbol{u}_n)\). Voer dan in die \(i\)-de iterasie die volgende prosedure uit.

  1. Bereken by die huidige verplasing \(\boldsymbol{u}_n + \Delta\boldsymbol{u}\) die tangensiële styfheid \(\boldsymbol{K}_i = \partial\boldsymbol{Q}/\partial\boldsymbol{u}|_{\boldsymbol{u}_n + \Delta\boldsymbol{u}}\) met die prosedure in Konstruksie van die tangensiële styfheidsmatriks.
  2. Om geometriese randvoorwaardes af te dwing, wysig die tangensiële styfheidsmatriks en residuvektor vir die vryheidsgrade wat aan verplasingsbeperkings onderworpe is, en verkry \(\tilde{\boldsymbol{K}}_i, \tilde{\boldsymbol{R}}_{i-1}\) (Behandeling van geometriese randvoorwaardes).
  3. Los die lineêre vergelyking \(\tilde{\boldsymbol{K}}_i\, d\boldsymbol{u}_i = \tilde{\boldsymbol{R}}_{i-1}\) op om die korreksie \(d\boldsymbol{u}_i\) te verkry. Hierdie prosedure maak dikwels die grootste deel van die berekeningskoste van die iteratiewe berekening uit.
  4. Werk die verplasingstoename by as \(\Delta\boldsymbol{u} \leftarrow \Delta\boldsymbol{u} + d\boldsymbol{u}_i\), en bereken dienooreenkomstig die internekragvektor \(\boldsymbol{Q}(\boldsymbol{u}_n + \Delta\boldsymbol{u})\) en residu \(\boldsymbol{R}_i = \boldsymbol{F}_{n+1} - \boldsymbol{Q}(\boldsymbol{u}_n + \Delta\boldsymbol{u})\).
  5. Kontroleer konvergensie en beëindig die iterasie indien konvergensie bereik is. Komponente wat met beperkte reaksies ooreenstem, verskyn in die residu \(\boldsymbol{R}_i\) by vryheidsgrade wat aan geometriese randvoorwaardes onderworpe is, dus word die konvergensie-aanwyser uit \(\tilde{\boldsymbol{R}}_i\) gekonstrueer nadat hierdie komponente uitgesluit is. Spesifieke konvergensie-aanwysers en drempels word in Konvergensiekriteria beskryf. Indien konvergensie nie bereik word nie en die iterasielimiet bereik word, word die iterasie as misluk behandel.

Wanneer die iterasie konvergeer, voeg die gekonvergeerde \(\Delta\boldsymbol{u}\) by \(\boldsymbol{u}_n\) om die opgehoopte verplasing op tyd \(t_{n+1}\), \(\boldsymbol{u}_{n+1} = \boldsymbol{u}_n + \Delta\boldsymbol{u}\), te verkry en gaan voort na die volgende tydstap.

Verwante onderwerpe