Gå til innhold

Newton-Raphson-metoden

Linearisering og iterativ rekurrens

Den ikke-lineære ligningen for knutepunktsforskyvningen ved tiden \(t_{n+1}\), \(\boldsymbol{u}_{n+1}\), som ble utledet i Virtuelt arbeid fra ytre krefter og sammenstilling av den globale ligningen, løses med Newton-Raphson-metoden. Knutepunktsforskyvningen frem til tiden \(t_n\), \(\boldsymbol{u}_n\), antas kjent, og forskyvningsinkrementet \(\Delta\boldsymbol{u}\) tas som den ukjente variabelen slik at

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

I det følgende neglisjeres den ytre kraftvektorens avhengighet av knutepunktsforskyvningen, og med \(\boldsymbol{F}(\boldsymbol{u}_{n+1}) = \boldsymbol{F}_{n+1}\) løses

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

ligningen.

Ved den gjeldende løsningen \(\Delta\boldsymbol{u}\) defineres tangentstivheten

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

Når den ikke-lineære ligningen lineariseres ved hjelp av denne, fås

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

La korreksjonen i den \(i\)-te iterasjonen være \(d\boldsymbol{u}_i\), og la residualvektoren ved starten av iterasjonen være

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

Da blir den iterative rekurrensen

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

Residualen \(\boldsymbol{R}_i\) er dermed en størrelse som representerer kraftubalansen fra likevekt.

Oppbygging av tangentstivhetsmatrisen

Tangentstivheten \(\boldsymbol{K} = \partial\boldsymbol{Q}/\partial\boldsymbol{u}\) bygges opp ved å partiell-derivere elementets indre kraftvektor, utledet i Diskretisering av virtuelt arbeid fra indre krefter, med hensyn på knutepunktsforskyvningen, integrere de resulterende integrandene for hvert element over elementområdet og sammenstille dem. Betegnes integranden på elementnivå med \(\boldsymbol{K}^e_X\) (referansekonfigurasjon, TL-metoden) eller \(\boldsymbol{K}^e_x\) (gjeldende konfigurasjon, UL-metoden), er elementets tangentstivhet

\[ \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 vises sluttformene for integrandene i TL-/UL-metoden. I begge tilfeller dekomponeres de i summen av et materialstivhetsledd (initialforskyvningsledd) og et geometrisk stivhetsledd (initialspenningsledd).

Total Lagrange-metoden

I Total Lagrange-metoden antas en lineær sammenheng mellom hastigheten til den andre Piola-Kirchhoff-spenningen \(\dot{\boldsymbol{S}}\) og Green-Lagrange-tøyningshastigheten \(\dot{\boldsymbol{E}}\), nemlig \(\dot{\boldsymbol{S}} = \boldsymbol{\mathsf{C}}:\dot{\boldsymbol{E}}\). Dette tilsvarer konstitutivlover for lineærelastiske materialer (St. Venant-Kirchhoff-materialer) og hyperelastiske materialer, og FrontISTR bruker Total Lagrange-metoden for disse materialene. Integranden for elementets tangentstivhet 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 leddet på høyre side er materialstivhetsleddet (initialforskyvningsleddet), og det andre er det geometriske stivhetsleddet (initialspenningsleddet).

I FrontISTR-implementasjonen beregnes denne integranden i matriseform med Voigt-notasjon:

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

Matrisene er som følger. \(\boldsymbol{B}_L, \boldsymbol{B}_{NL}\) er B-matrisene som ble introdusert i Diskretisering av virtuelt arbeid fra indre krefter, og \(\tilde{\boldsymbol{C}}\) er Voigt-representasjonen av konstitutivtensoren \(\boldsymbol{\mathsf{C}}\), det vil si en \(6\times 6\) materialstivhetsmatrise (Tensornotasjon og matematisk grunnlag). \(\boldsymbol{S}_9, \boldsymbol{F}_9\) er omordningsmatrisene nedenfor, som brukes til å uttrykke det geometriske stivhetsleddet som et matriseprodukt. Først defineres notasjonen \([\,\cdot\,]\), som omordner en \(3\times 3\) andreordenstensor \(\boldsymbol{A}\) til en 9-komponentvektor, som

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

Med denne definisjonen er \(\boldsymbol{F}_9\) en \(9\times d n_e\)-matrise som uttrykker variasjonen av deformasjonsgradienten som \([\delta\boldsymbol{F}] = \boldsymbol{F}_9\, \delta\boldsymbol{u}^e\). For elementknutepunkt \(\alpha = 1, \ldots, n_e\) er den tilsvarende \(9\times d\)-blokken

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

og \(\boldsymbol{F}_9 = [[\boldsymbol{F}_9]_1, \ldots, [\boldsymbol{F}_9]_{n_e}]\) fås ved å plassere blokkene horisontalt i elementknutepunktsrekkefølge. \(\boldsymbol{S}_9\) er en \(9\times 9\)-matrise valgt slik at det geometriske stivhetsleddet, i kombinasjon med denne matrisen, uttrykkes som \(\delta\boldsymbol{u}^{eT}\, \boldsymbol{F}_9^T \boldsymbol{S}_9 \boldsymbol{F}_9\, \dot{\boldsymbol{u}}^e\):

\[ \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 matrisen.

Updated Lagrange-metoden

I Updated Lagrange-metoden antas en lineær sammenheng mellom Jaumann-hastigheten til den relative Kirchhoff-spenningstensoren \(\hat{\boldsymbol{\sigma}}^{\nabla J}\) og deformasjonshastighetstensoren \(\boldsymbol{D}\), nemlig \(\hat{\boldsymbol{\sigma}}^{\nabla J} = \boldsymbol{\mathsf{C}}:\boldsymbol{D}\). Dette er formen til en hypoelastisk konstitutivlov som er felles for lineærelastiske, elastoplastiske og krypende materialer, og FrontISTR bruker Updated Lagrange-metoden for disse materialene. Integranden for elementets tangentstivhet uttrykt i den gjeldende konfigurasjonen 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}) \]

der \(\boldsymbol{\sigma}^{\nabla T}\) er Truesdell-hastigheten, \(\boldsymbol{A}_{(L)}\) er den lineære delen av Almansi-tøyningen, \(\boldsymbol{F}_t = \partial\boldsymbol{u}/\partial\boldsymbol{x}\) er forskyvningsgradienten med hensyn på gjeldende konfigurasjon, og \(\boldsymbol{L}\) er hastighetsgradienttensoren. Det første leddet på høyre side er materialstivhetsleddet, og det andre er det geometriske stivhetsleddet.

I FrontISTR-implementasjonen beregnes denne integranden i matriseform med Voigt-notasjon:

\[ \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-matrisen konstruert i den gjeldende konfigurasjonen (Diskretisering av virtuelt arbeid fra indre krefter). \(\boldsymbol{\sigma}_9, \boldsymbol{f}_9\) fås fra \(\boldsymbol{S}_9, \boldsymbol{F}_9\), som ble definert for TL-metoden, ved å erstatte den andre PK-spenningen \(\boldsymbol{S}\) med Cauchy-spenningen \(\boldsymbol{\sigma}\) og gradienten i referansekonfigurasjonen \(\partial N_\alpha^e/\partial X_i\) med gradienten i gjeldende konfigurasjon \(\partial N_\alpha^e/\partial x_i\).

\(\boldsymbol{G}\) er en Cauchy-spenningsavhengig korreksjonsmatrise som er nødvendig for å gjøre den hypoelastiske konstitutivloven \(\hat{\boldsymbol{\sigma}}^{\nabla J} = \boldsymbol{\mathsf{C}}:\boldsymbol{D}\) konsistent med rammeverket for tangentstivhet som en konstitutivlov basert på Truesdell-hastigheten. Den fås ved å ordne komponentene i fjerdeordenstensoren \(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 matrisen.

Sammenstilling av den globale stivhetsmatrisen

Den globale tangentstivheten \(\boldsymbol{K}\) fås ved å dele hvert elements stivhet \(\boldsymbol{K}^e\) i \(d\times d\)-blokker \(\boldsymbol{K}^e_{\alpha\beta}\) for hvert par av knutepunkter og bruke sammenstillingsmengden for andreordenstensorer \(\mathcal{E}^2(i_g, i_h)\) som ble introdusert i Sammenstilling av fysiske størrelser i elementknutepunkter:

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

De resulterende verdiene ordnes i en matrise med rad \(i_g\) og kolonne \(i_h\). I implementasjonen konstrueres ikke mengden \(\mathcal{E}^2\) eksplisitt; de tilsvarende blokkene legges i stedet direkte til inne i elementløkken. Matrisen er kvadratisk med dimensjon lik antall frihetsgrader per knutepunkt \(\times\) totalt antall knutepunkter \(n_g\), men siden komponenter mellom andre knutepunkter enn dem som er koblet gjennom elementer, er \(0\), lagres den som en spredt matrise.

Elementstivhetene for TL- og UL-metoden har samme form bortsett fra bytte av referansekonfigurasjon (knutepunktskoordinater og grunnlaget for konstruksjon av B-matrisen) og om \(\boldsymbol{G}\)-matrisen er til stede. FrontISTR implementerer derfor begge metodene i en felles subrutine.

Iterasjonsalgoritme

Oppsummert settes \(\Delta\boldsymbol{u} = \boldsymbol{0}\) ved starten av iterasjonen, og den innledende residualen \(\boldsymbol{R}_0 = \boldsymbol{F}_{n+1} - \boldsymbol{Q}(\boldsymbol{u}_n)\) beregnes. I den \(i\)-te iterasjonen utføres deretter følgende trinn.

  1. Ved den gjeldende forskyvningen \(\boldsymbol{u}_n + \Delta\boldsymbol{u}\) beregnes tangentstivheten \(\boldsymbol{K}_i = \partial\boldsymbol{Q}/\partial\boldsymbol{u}|_{\boldsymbol{u}_n + \Delta\boldsymbol{u}}\) etter fremgangsmåten i Oppbygging av tangentstivhetsmatrisen.
  2. For å innføre geometriske randbetingelser modifiseres tangentstivhetsmatrisen og residualvektoren for frihetsgradene som har forskyvningsbetingelser, slik at \(\tilde{\boldsymbol{K}}_i, \tilde{\boldsymbol{R}}_{i-1}\) fås (Behandling av geometriske randbetingelser).
  3. Løs den lineære ligningen \(\tilde{\boldsymbol{K}}_i\, d\boldsymbol{u}_i = \tilde{\boldsymbol{R}}_{i-1}\) for å finne korreksjonen \(d\boldsymbol{u}_i\). Dette trinnet står ofte for størstedelen av beregningskostnaden i iterasjonen.
  4. Oppdater forskyvningsinkrementet som \(\Delta\boldsymbol{u} \leftarrow \Delta\boldsymbol{u} + d\boldsymbol{u}_i\), og beregn deretter den indre kraftvektoren \(\boldsymbol{Q}(\boldsymbol{u}_n + \Delta\boldsymbol{u})\) og residualen \(\boldsymbol{R}_i = \boldsymbol{F}_{n+1} - \boldsymbol{Q}(\boldsymbol{u}_n + \Delta\boldsymbol{u})\).
  5. Kontroller konvergens og avslutt iterasjonen dersom konvergens er oppnådd. I residualen \(\boldsymbol{R}_i\) opptrer komponenter som tilsvarer reaksjonskrefter ved frihetsgrader med geometriske randbetingelser, så konvergensindikatoren bygges fra \(\tilde{\boldsymbol{R}}_i\) etter at disse komponentene er utelatt. De konkrete konvergensindikatorene og tersklene beskrives i Konvergenskriterier. Hvis konvergens ikke oppnås før iterasjonsgrensen nås, regnes iterasjonen som mislykket.

Når iterasjonen har konvergert, legges den konvergerte \(\Delta\boldsymbol{u}\) til \(\boldsymbol{u}_n\) for å få den akkumulerte forskyvningen ved tiden \(t_{n+1}\), \(\boldsymbol{u}_{n+1} = \boldsymbol{u}_n + \Delta\boldsymbol{u}\), og beregningen går videre til neste tidssteg.

Relaterte emner

AI-assisted translation May contain errors Official docs Status