Skoči na vsebino

Newton-Raphsonova metoda

Linearizacija in iterativna rekurzija

Virtualno delo zunanjih sil in sestavljanje globalne enačbe podaja nelinearno enačbo za vozliščni pomik v času \(t_{n+1}\), \(\boldsymbol{u}_{n+1}\), ki se rešuje z Newton-Raphsonovo metodo. Predpostavimo, da je vozliščni pomik do časa \(t_n\), \(\boldsymbol{u}_n\), znan, prirastek pomika \(\Delta\boldsymbol{u}\) pa vzamemo kot neznano spremenljivko, ki jo je treba določiti

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

V nadaljevanju zanemarimo odvisnost vektorja zunanjih sil od vozliščnega pomika in ob \(\boldsymbol{F}(\boldsymbol{u}_{n+1}) = \boldsymbol{F}_{n+1}\) rešujemo

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

Pri trenutni rešitvi \(\Delta\boldsymbol{u}\) definiramo tangentno togost

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

Z njeno uporabo linearizacija nelinearne enačbe da

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

Naj bo popravek v \(i\)-ti iteraciji \(d\boldsymbol{u}_i\), rezidualni vektor na začetku iteracije pa

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

Tedaj je iterativna rekurzija

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

Rezidual \(\boldsymbol{R}_i\) je torej količina, ki ustreza neravnotežju sil glede na ravnotežno stanje.

Konstrukcija tangentne matrike togosti

Tangentna togost \(\boldsymbol{K} = \partial\boldsymbol{Q}/\partial\boldsymbol{u}\) se sestavi tako, da se elementni vektor notranjih sil, dobljen v Diskretizaciji virtualnega dela notranjih sil, parcialno odvaja po vozliščnem pomiku, nastali integrandi na ravni elementa integrirajo po domeni vsakega elementa in sestavijo. Če integrand na ravni elementa označimo z \(\boldsymbol{K}^e_X\) (zapis v referenčni konfiguraciji, formulacija TL) ali \(\boldsymbol{K}^e_x\) (zapis v trenutni konfiguraciji, formulacija UL), je elementna tangentna togost

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

V nadaljevanju so podane končne oblike integrandov TL/UL. V obeh primerih sta razstavljena na vsoto člena materialne togosti (člen začetnega pomika) in člena geometrijske togosti (člen začetne napetosti).

Totalna Lagrangeeva formulacija

V totalni Lagrangeevi formulaciji se predpostavi linearna zveza med hitrostjo druge Piola-Kirchhoffove napetosti \(\dot{\boldsymbol{S}}\) in hitrostjo Green-Lagrangeeve deformacije \(\dot{\boldsymbol{E}}\), in sicer \(\dot{\boldsymbol{S}} = \boldsymbol{\mathsf{C}}:\dot{\boldsymbol{E}}\). To ustreza konstitutivnim zakonom za linearno elastične materiale (St. Venant-Kirchhoffove materiale) in hiperelastične materiale; FrontISTR za te materiale uporablja totalno Lagrangeevo formulacijo. Integrand elementne tangentne togosti je nato v tenzorski obliki zapisan kot

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

Prvi člen na desni je člen materialne togosti (člen začetnega pomika), drugi pa člen geometrijske togosti (člen začetne napetosti).

V implementaciji FrontISTR se ta integrand ovrednoti v matrični obliki z Voigtovim zapisom:

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

Matrike so naslednje. \(\boldsymbol{B}_L, \boldsymbol{B}_{NL}\) sta matriki B, uvedeni v Diskretizaciji virtualnega dela notranjih sil, \(\tilde{\boldsymbol{C}}\) pa je Voigtov zapis konstitutivnega tenzorja \(\boldsymbol{\mathsf{C}}\), torej materialna matrika togosti velikosti \(6\times 6\) (Tenzorski zapis in matematične osnove). \(\boldsymbol{S}_9, \boldsymbol{F}_9\) sta preureditveni matriki, uporabljeni za zapis člena geometrijske togosti kot matričnega produkta. Najprej za tenzor drugega reda \(\boldsymbol{A}\) velikosti \(3\times 3\) definiramo zapis \([\,\cdot\,]\), ki ga preuredi v 9-komponentni vektor, kot

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

S to definicijo \(\boldsymbol{F}_9\) izraža variacijo deformacijskega gradienta v obliki \([\delta\boldsymbol{F}] = \boldsymbol{F}_9\, \delta\boldsymbol{u}^e\) in je matrika velikosti \(9\times d n_e\). Za elementno vozlišče \(\alpha = 1, \ldots, n_e\) je pripadajoči blok velikosti \(9\times d\)

\[ [\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{ je identična matrika } 3\times 3) \]

in velja \(\boldsymbol{F}_9 = [[\boldsymbol{F}_9]_1, \ldots, [\boldsymbol{F}_9]_{n_e}]\), pri čemer so bloki vodoravno razporejeni po vrstnem redu elementnih vozlišč. \(\boldsymbol{S}_9\) je izbrana tako, da je člen geometrijske togosti v kombinaciji s to matriko izražen kot \(\delta\boldsymbol{u}^{eT}\, \boldsymbol{F}_9^T \boldsymbol{S}_9 \boldsymbol{F}_9\, \dot{\boldsymbol{u}}^e\); gre za naslednjo matriko velikosti \(9\times 9\)

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

To je dobljena matrika.

Posodobljena Lagrangeeva formulacija

V posodobljeni Lagrangeevi formulaciji se predpostavi linearna zveza med Jaumannovo hitrostjo relativnega Kirchhoffovega napetostnega tenzorja \(\hat{\boldsymbol{\sigma}}^{\nabla J}\) in tenzorjem hitrosti deformacije \(\boldsymbol{D}\), in sicer \(\hat{\boldsymbol{\sigma}}^{\nabla J} = \boldsymbol{\mathsf{C}}:\boldsymbol{D}\). To je oblika hipoelastičnega konstitutivnega zakona, skupna linearno elastičnim, elastoplastičnim in materialom z lezenjem, FrontISTR pa za te materiale uporablja posodobljeno Lagrangeevo formulacijo. Integrand elementne tangentne togosti, izražen v trenutni konfiguraciji, je nato v tenzorski obliki zapisan kot

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

kjer je \(\boldsymbol{\sigma}^{\nabla T}\) Truesdellova hitrost, \(\boldsymbol{A}_{(L)}\) linearni del Almansijeve deformacije, \(\boldsymbol{F}_t = \partial\boldsymbol{u}/\partial\boldsymbol{x}\) gradient pomika glede na trenutno konfiguracijo, \(\boldsymbol{L}\) pa tenzor gradienta hitrosti. Prvi člen na desni je člen materialne togosti, drugi pa člen geometrijske togosti.

V implementaciji FrontISTR se ta integrand ovrednoti v matrični obliki z Voigtovim zapisom:

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

Tu je \(\boldsymbol{b}\) matrika B, sestavljena v trenutni konfiguraciji (Diskretizacija virtualnega dela notranjih sil). \(\boldsymbol{\sigma}_9, \boldsymbol{f}_9\) dobimo iz \(\boldsymbol{S}_9, \boldsymbol{F}_9\), definiranih za formulacijo TL, tako da drugo PK-napetost \(\boldsymbol{S}\) nadomestimo s Cauchyjevo napetostjo \(\boldsymbol{\sigma}\), gradient v referenčni konfiguraciji \(\partial N_\alpha^e/\partial X_i\) pa z gradientom v trenutni konfiguraciji \(\partial N_\alpha^e/\partial x_i\).

\(\boldsymbol{G}\) je korekcijska matrika, odvisna od Cauchyjeve napetosti, potrebna zato, da je hipoelastični konstitutivni zakon \(\hat{\boldsymbol{\sigma}}^{\nabla J} = \boldsymbol{\mathsf{C}}:\boldsymbol{D}\) skladen z ogrodjem tangentne togosti kot konstitutivni zakon na osnovi Truesdellove hitrosti. Dobimo jo z razporeditvijo komponent tenzorja četrtega reda \(G_{ijkl} = \delta_{il}\sigma_{kj} + \delta_{jl}\sigma_{ik}\) v Voigtovo obliko \(6\times 6\) kot

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

To je dobljena matrika.

Sestavljanje globalne matrike togosti

Globalna tangentna togost \(\boldsymbol{K}\) se dobi tako, da se vsaka elementna togost \(\boldsymbol{K}^e\) razdeli na bloke \(\boldsymbol{K}^e_{\alpha\beta}\) velikosti \(d\times d\) za vsak par vozlišč in uporabi množica za sestavljanje tenzorjev drugega reda \(\mathcal{E}^2(i_g, i_h)\), uvedena v Sestavljanju fizikalnih količin elementnih vozlišč:

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

Dobljene vrednosti so razporejene kot matrika z vrstico \(i_g\) in stolpcem \(i_h\). V implementaciji se množica \(\mathcal{E}^2\) ne sestavi eksplicitno; namesto tega se ustrezni bloki dodajo neposredno znotraj elementne zanke. Matrika je kvadratna z dimenzijo, enako številu prostostnih stopenj na vozlišče \(\times\) skupno število vozlišč \(n_g\), vendar so komponente razen tistih med vozlišči, povezanimi prek elementov, enake \(0\), zato je shranjena v obliki redke matrike.

Elementne matrike togosti za formulaciji TL in UL imajo enako obliko, razen preklopa referenčne konfiguracije (vozliščne koordinate in vir, uporabljen za sestavo matrike B) ter prisotnosti oziroma odsotnosti matrike \(\boldsymbol{G}\). FrontISTR zato obe formulaciji implementira v skupnem podprogramu.

Iteracijski algoritem

Če povzamemo zgoraj navedeno: na začetku iteracije nastavite \(\Delta\boldsymbol{u} = \boldsymbol{0}\) in izračunajte začetni rezidual \(\boldsymbol{R}_0 = \boldsymbol{F}_{n+1} - \boldsymbol{Q}(\boldsymbol{u}_n)\). Nato v \(i\)-ti iteraciji izvedite naslednji postopek.

  1. Pri trenutnem pomiku \(\boldsymbol{u}_n + \Delta\boldsymbol{u}\) izračunajte tangentno togost \(\boldsymbol{K}_i = \partial\boldsymbol{Q}/\partial\boldsymbol{u}|_{\boldsymbol{u}_n + \Delta\boldsymbol{u}}\) po postopku iz Konstrukcije tangentne matrike togosti.
  2. Za uveljavitev geometrijskih robnih pogojev spremenite tangentno matriko togosti in rezidualni vektor za prostostne stopnje, na katerih veljajo omejitve pomika, ter dobite \(\tilde{\boldsymbol{K}}_i, \tilde{\boldsymbol{R}}_{i-1}\) (Obravnava geometrijskih robnih pogojev).
  3. Rešite linearno enačbo \(\tilde{\boldsymbol{K}}_i\, d\boldsymbol{u}_i = \tilde{\boldsymbol{R}}_{i-1}\) in določite popravek \(d\boldsymbol{u}_i\). Ta postopek pogosto predstavlja večino računskega stroška iterativnega izračuna.
  4. Posodobite prirastek pomika kot \(\Delta\boldsymbol{u} \leftarrow \Delta\boldsymbol{u} + d\boldsymbol{u}_i\) in skladno s tem izračunajte vektor notranjih sil \(\boldsymbol{Q}(\boldsymbol{u}_n + \Delta\boldsymbol{u})\) ter rezidual \(\boldsymbol{R}_i = \boldsymbol{F}_{n+1} - \boldsymbol{Q}(\boldsymbol{u}_n + \Delta\boldsymbol{u})\).
  5. Preverite konvergenco in končajte iteracijo, če je dosežena. Komponente, ki ustrezajo reakcijam vpetja, se pojavijo v rezidualu \(\boldsymbol{R}_i\) na prostostnih stopnjah, na katerih veljajo geometrijski robni pogoji, zato se konvergenčni kazalnik sestavi iz \(\tilde{\boldsymbol{R}}_i\) po izločitvi teh komponent. Posamezni konvergenčni kazalniki in pragovi so opisani v Konvergenčnih kriterijih. Če konvergenca ni dosežena in je dosežena meja iteracij, se iteracija obravnava kot neuspešna.

Ko iteracija konvergira, prištejte konvergirani \(\Delta\boldsymbol{u}\) k \(\boldsymbol{u}_n\), da dobite akumulirani pomik v času \(t_{n+1}\), \(\boldsymbol{u}_{n+1} = \boldsymbol{u}_n + \Delta\boldsymbol{u}\), in nadaljujte z naslednjim časovnim korakom.

Sorodne teme

AI-assisted translation May contain errors Official docs Status