Gå til indholdet

Newton-Raphson-metoden

Linearisering og iterativ rekursion

Virtuelt arbejde fra ydre kræfter og samling af den globale ligning giver en ikke-lineær ligning for knudeforskydningen ved tiden \(t_{n+1}\), \(\boldsymbol{u}_{n+1}\), som løses med Newton-Raphson-metoden. Knudeforskydningen frem til tiden \(t_n\), \(\boldsymbol{u}_n\), antages kendt, og forskydningsinkrementet \(\Delta\boldsymbol{u}\) betragtes som den ukendte variabel, der skal bestemmes

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

I det følgende negligeres den ydre kraftvektors afhængighed af knudeforskydningen, og med \(\boldsymbol{F}(\boldsymbol{u}_{n+1}) = \boldsymbol{F}_{n+1}\)

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

løses ligningen.

Ved den aktuelle løsning \(\Delta\boldsymbol{u}\) defineres den tangentielle stivhed

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

Ved hjælp af denne giver linearisering af den ikke-lineære ligning

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

Lad korrektionen i den \(i\)-te iteration være \(d\boldsymbol{u}_i\), og lad residualvektoren ved iterationens start være

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

Da er den 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 \]

Residualet \(\boldsymbol{R}_i\) er således en størrelse, der svarer til kraftubalancen fra ligevægt.

Konstruktion af den tangentielle stivhedsmatrix

Den tangentielle stivhed \(\boldsymbol{K} = \partial\boldsymbol{Q}/\partial\boldsymbol{u}\) konstrueres ved at partieltdifferentiere elementets indre kraftvektor, som er udledt i Diskretisering af det virtuelle arbejde fra indre kræfter, med hensyn til knudeforskydningen, integrere de resulterende integrander på elementniveau over hvert elementdomæne og samle dem. Betegnes integranden på elementniveau med \(\boldsymbol{K}^e_X\) (notation i referencekonfigurationen, TL-formulering) eller \(\boldsymbol{K}^e_x\) (notation i den aktuelle konfiguration, UL-formulering), er elementets tangentielle stivhed

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

Nedenfor gives de endelige former for TL/UL-integranderne. I begge tilfælde opdeles de i summen af et materialestivhedsled (initialforskydningsled) og et geometrisk stivhedsled (initialspændingsled).

Total Lagrange-formulering

I Total Lagrange-formuleringen antages en lineær relation mellem hastigheden af den anden Piola-Kirchhoff-spænding \(\dot{\boldsymbol{S}}\) og Green-Lagrange-tøjningshastigheden \(\dot{\boldsymbol{E}}\), nemlig \(\dot{\boldsymbol{S}} = \boldsymbol{\mathsf{C}}:\dot{\boldsymbol{E}}\). Dette svarer til konstitutive love for lineærelastiske materialer (St. Venant-Kirchhoff-materialer) og hyperelastiske materialer, og FrontISTR anvender Total Lagrange-formuleringen for disse materialer. Integranden for elementets tangentielle stivhed skrives da i tensorform som

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

Det første led på højre side er materialestivhedsleddet (initialforskydningsleddet), og det andet led er det geometriske stivhedsled (initialspændingsleddet).

I FrontISTR-implementeringen evalueres denne integrand i matrixform ved hjælp af Voigt-notation:

\[ \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 enkelte matricer er som følger. \(\boldsymbol{B}_L, \boldsymbol{B}_{NL}\) er B-matricerne introduceret i Diskretisering af det virtuelle arbejde fra indre kræfter, og \(\tilde{\boldsymbol{C}}\) er Voigt-repræsentationen af den konstitutive tensor \(\boldsymbol{\mathsf{C}}\), dvs. en \(6\times 6\) materialestivhedsmatrix (Tensor-notation og matematiske grundlag). \(\boldsymbol{S}_9, \boldsymbol{F}_9\) er de følgende omordningsmatricer, der bruges til at udtrykke det geometriske stivhedsled som et matrixprodukt. For en \(3\times 3\) tensor af anden orden \(\boldsymbol{A}\) defineres først notationen \([\,\cdot\,]\), som omordner den til en vektor med 9 komponenter, som

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

Med denne definition udtrykker \(\boldsymbol{F}_9\) variationen af deformationsgradienten på formen \([\delta\boldsymbol{F}] = \boldsymbol{F}_9\, \delta\boldsymbol{u}^e\) og er en \(9\times d n_e\)-matrix. For elementknuden \(\alpha = 1, \ldots, n_e\) er den tilsvarende \(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{ er } 3\times 3 \text{ identitetsmatricen}) \]

og gives ved \(\boldsymbol{F}_9 = [[\boldsymbol{F}_9]_1, \ldots, [\boldsymbol{F}_9]_{n_e}]\), hvor blokkene er arrangeret vandret i elementknudernes rækkefølge. \(\boldsymbol{S}_9\) vælges således, at det geometriske stivhedsled, når den kombineres med denne matrix, udtrykkes som \(\delta\boldsymbol{u}^{eT}\, \boldsymbol{F}_9^T \boldsymbol{S}_9 \boldsymbol{F}_9\, \dot{\boldsymbol{u}}^e\); den er følgende \(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} \]

Dette er den resulterende matrix.

Updated Lagrange-formulering

I Updated Lagrange-formuleringen antages en lineær relation mellem Jaumann-hastigheden af den relative Kirchhoff-spændingstensor \(\hat{\boldsymbol{\sigma}}^{\nabla J}\) og deformationshastighedstensoren \(\boldsymbol{D}\), nemlig \(\hat{\boldsymbol{\sigma}}^{\nabla J} = \boldsymbol{\mathsf{C}}:\boldsymbol{D}\). Dette er formen af en hypoelastisk konstitutiv lov, der er fælles for lineærelastiske, elastoplastiske og krybende materialer, og FrontISTR anvender Updated Lagrange-formuleringen for disse materialer. Integranden for elementets tangentielle stivhed udtrykt i den aktuelle konfiguration skrives da i tensorform som

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

hvor \(\boldsymbol{\sigma}^{\nabla T}\) er Truesdell-hastigheden, \(\boldsymbol{A}_{(L)}\) er den lineære del af Almansi-tøjningen, \(\boldsymbol{F}_t = \partial\boldsymbol{u}/\partial\boldsymbol{x}\) er forskydningsgradienten med hensyn til den aktuelle konfiguration, og \(\boldsymbol{L}\) er hastighedsgradienttensoren. Det første led på højre side er materialestivhedsleddet, og det andet led er det geometriske stivhedsled.

I FrontISTR-implementeringen evalueres denne integrand i matrixform ved hjælp af Voigt-notation:

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

Her er \(\boldsymbol{b}\) B-matricen konstrueret i den aktuelle konfiguration (Diskretisering af det virtuelle arbejde fra indre kræfter). \(\boldsymbol{\sigma}_9, \boldsymbol{f}_9\) fås fra \(\boldsymbol{S}_9, \boldsymbol{F}_9\), som er defineret for TL-formuleringen, ved at erstatte den anden PK-spænding \(\boldsymbol{S}\) med Cauchy-spændingen \(\boldsymbol{\sigma}\) og gradienten i referencekonfigurationen \(\partial N_\alpha^e/\partial X_i\) med gradienten i den aktuelle konfiguration \(\partial N_\alpha^e/\partial x_i\).

\(\boldsymbol{G}\) er en korrektionsmatrix, der afhænger af Cauchy-spændingen, og som er nødvendig for at gøre den hypoelastiske konstitutive lov \(\hat{\boldsymbol{\sigma}}^{\nabla J} = \boldsymbol{\mathsf{C}}:\boldsymbol{D}\) konsistent med rammen for tangentiel stivhed som en konstitutiv lov baseret på Truesdell-hastigheden. Den fås ved at arrangere fjerdeordens-tensorkomponenterne \(G_{ijkl} = \delta_{il}\sigma_{kj} + \delta_{jl}\sigma_{ik}\) i \(6\times 6\) Voigt-form som

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

Dette er den resulterende matrix.

Samling af den globale stivhedsmatrix

Den globale tangentielle stivhed \(\boldsymbol{K}\) fås ved at opdele hver elementstivhed \(\boldsymbol{K}^e\) i \(d\times d\)-blokke \(\boldsymbol{K}^e_{\alpha\beta}\) for hvert knudepar og anvende samlingsmængden for andenordens-tensorer \(\mathcal{E}^2(i_g, i_h)\), som introduceres i Samling af fysiske størrelser ved elementknuder:

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

De resulterende værdier arrangeres som en matrix med række \(i_g\) og kolonne \(i_h\). I implementeringen konstrueres mængden \(\mathcal{E}^2\) ikke eksplicit; i stedet lægges de tilsvarende blokke direkte til inde i elementløkken. Matricen er kvadratisk med en dimension svarende til frihedsgrader pr. knude \(\times\) det samlede antal knuder \(n_g\), men da komponenter bortset fra dem mellem knuder, der er forbundet gennem elementer, er \(0\), lagres den i sparsom matrixform.

Elementstivhedsmatricerne for TL- og UL-formuleringerne har samme form bortset fra skiftet af referencekonfiguration (knudekoordinater og det grundlag, der bruges til at konstruere B-matricen) og tilstedeværelsen eller fraværet af \(\boldsymbol{G}\)-matricen. FrontISTR implementerer derfor begge formuleringer i en fælles subrutine.

Iterationsalgoritme

Sammenfattende sættes \(\Delta\boldsymbol{u} = \boldsymbol{0}\) ved iterationens start, og det indledende residual \(\boldsymbol{R}_0 = \boldsymbol{F}_{n+1} - \boldsymbol{Q}(\boldsymbol{u}_n)\) beregnes. Derefter udføres følgende procedure i den \(i\)-te iteration.

  1. Ved den aktuelle forskydning \(\boldsymbol{u}_n + \Delta\boldsymbol{u}\) beregnes den tangentielle stivhed \(\boldsymbol{K}_i = \partial\boldsymbol{Q}/\partial\boldsymbol{u}|_{\boldsymbol{u}_n + \Delta\boldsymbol{u}}\) ved hjælp af proceduren i Konstruktion af den tangentielle stivhedsmatrix.
  2. For at pålægge geometriske randbetingelser modificeres den tangentielle stivhedsmatrix og residualvektoren for de frihedsgrader, der er underlagt forskydningsbetingelser, så \(\tilde{\boldsymbol{K}}_i, \tilde{\boldsymbol{R}}_{i-1}\) opnås (Behandling af geometriske randbetingelser).
  3. Løs den lineære ligning \(\tilde{\boldsymbol{K}}_i\, d\boldsymbol{u}_i = \tilde{\boldsymbol{R}}_{i-1}\) for at få korrektionen \(d\boldsymbol{u}_i\). Denne procedure udgør ofte størstedelen af beregningsomkostningen ved den iterative beregning.
  4. Opdater forskydningsinkrementet som \(\Delta\boldsymbol{u} \leftarrow \Delta\boldsymbol{u} + d\boldsymbol{u}_i\), og beregn tilsvarende den indre kraftvektor \(\boldsymbol{Q}(\boldsymbol{u}_n + \Delta\boldsymbol{u})\) og residualet \(\boldsymbol{R}_i = \boldsymbol{F}_{n+1} - \boldsymbol{Q}(\boldsymbol{u}_n + \Delta\boldsymbol{u})\).
  5. Kontroller konvergensen, og afslut iterationen, hvis der er opnået konvergens. Komponenter svarende til fastholdelsesreaktioner optræder i residualet \(\boldsymbol{R}_i\) ved frihedsgrader, der er underlagt geometriske randbetingelser, så konvergensindikatoren konstrueres ud fra \(\tilde{\boldsymbol{R}}_i\) efter udelukkelse af disse komponenter. Konkrete konvergensindikatorer og tærskler beskrives i Konvergenskriterier. Hvis der ikke opnås konvergens, og iterationsgrænsen nås, behandles iterationen som mislykket.

Når iterationen konvergerer, lægges det konvergerede \(\Delta\boldsymbol{u}\) til \(\boldsymbol{u}_n\) for at få den akkumulerede forskydning ved tiden \(t_{n+1}\), \(\boldsymbol{u}_{n+1} = \boldsymbol{u}_n + \Delta\boldsymbol{u}\), hvorefter der fortsættes til næste tidstrin.

Relaterede emner

AI-assisted translation May contain errors Official docs Status