Preskoči na sadržaj

Newton-Raphsonova metoda

Linearizacija i iteracijska rekurzija

Virtualni rad vanjskih sila i sastavljanje globalne jednadžbe daje nelinearnu jednadžbu za pomak čvorova u trenutku \(t_{n+1}\), \(\boldsymbol{u}_{n+1}\), koja se rješava Newton-Raphsonovom metodom. Pretpostavlja se da je pomak čvorova do trenutka \(t_n\), \(\boldsymbol{u}_n\), poznat, a prirast pomaka \(\Delta\boldsymbol{u}\) uzima se kao nepoznata varijabla koju treba odrediti

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

U nastavku se zanemaruje ovisnost vektora vanjskih sila o pomaku čvorova te, uz \(\boldsymbol{F}(\boldsymbol{u}_{n+1}) = \boldsymbol{F}_{n+1}\),

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

se rješava.

Za trenutačno rješenje \(\Delta\boldsymbol{u}\) definirajmo tangentnu krutost

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

Pomoću nje linearizacija nelinearne jednadžbe daje

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

Neka je korekcija u \(i\)-toj iteraciji \(d\boldsymbol{u}_i\), a vektor reziduala na početku iteracije

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

Tada je iteracijska rekurzija

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

Prema tome, rezidual \(\boldsymbol{R}_i\) predstavlja veličinu koja odgovara neravnoteži sila u odnosu na stanje ravnoteže.

Konstrukcija tangentne matrice krutosti

Tangentna krutost \(\boldsymbol{K} = \partial\boldsymbol{Q}/\partial\boldsymbol{u}\) konstruira se parcijalnim deriviranjem vektora unutarnjih sila elementa dobivenog u Diskretizaciji virtualnog rada unutarnjih sila po pomaku čvorova, integriranjem tako dobivenih integranada na razini elementa po domeni svakog elementa te njihovim sastavljanjem. Označimo li integrand na razini elementa s \(\boldsymbol{K}^e_X\) (notacija referentne konfiguracije, TL formulacija) ili \(\boldsymbol{K}^e_x\) (notacija trenutačne konfiguracije, UL formulacija), tangentna krutost elementa je

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

U nastavku su dane konačne forme TL/UL integranada. U oba slučaja rastavljaju se na zbroj člana materijalne krutosti (člana početnog pomaka) i člana geometrijske krutosti (člana početnog naprezanja).

Total Lagrange formulacija

U Total Lagrange formulaciji pretpostavlja se linearna veza između stope drugog Piola-Kirchhoffova naprezanja \(\dot{\boldsymbol{S}}\) i stope Green-Lagrangeove deformacije \(\dot{\boldsymbol{E}}\), odnosno \(\dot{\boldsymbol{S}} = \boldsymbol{\mathsf{C}}:\dot{\boldsymbol{E}}\). To odgovara konstitutivnim zakonima za linearno elastične materijale (St. Venant-Kirchhoffove materijale) i hiperelastične materijale, a FrontISTR za te materijale koristi Total Lagrange formulaciju. Integrand tangentne krutosti elementa tada se u tenzorskom obliku zapisuje kao

\[ \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 član na desnoj strani jest član materijalne krutosti (član početnog pomaka), a drugi je član geometrijske krutosti (član početnog naprezanja).

U implementaciji FrontISTR-a taj se integrand izračunava u matričnom obliku koristeći Voigtovu notaciju:

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

Pojedine matrice su sljedeće. \(\boldsymbol{B}_L, \boldsymbol{B}_{NL}\) su B-matrice uvedene u Diskretizaciji virtualnog rada unutarnjih sila, a \(\tilde{\boldsymbol{C}}\) je Voigtov prikaz konstitutivnog tenzora \(\boldsymbol{\mathsf{C}}\), odnosno \(6\times 6\) matrica materijalne krutosti (Tenzorska notacija i matematičke osnove). \(\boldsymbol{S}_9, \boldsymbol{F}_9\) su sljedeće matrice preslagivanja koje se koriste kako bi se član geometrijske krutosti izrazio kao matrični produkt. Najprije se za \(3\times 3\) tenzor drugog reda \(\boldsymbol{A}\) definira notacija \([\,\cdot\,]\) kojom se on preslaguje u vektor s 9 komponenti kao

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

S tom definicijom \(\boldsymbol{F}_9\) izražava varijaciju gradijenta deformacije u obliku \([\delta\boldsymbol{F}] = \boldsymbol{F}_9\, \delta\boldsymbol{u}^e\) i predstavlja matricu dimenzija \(9\times d n_e\). Za čvor elementa \(\alpha = 1, \ldots, n_e\) odgovarajući blok dimenzija \(9\times d\) je

\[ [\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 } 3\times 3 \text{ jedinična matrica}) \]

i vrijedi \(\boldsymbol{F}_9 = [[\boldsymbol{F}_9]_1, \ldots, [\boldsymbol{F}_9]_{n_e}]\), pri čemu su blokovi vodoravno poredani prema redoslijedu čvorova elementa. \(\boldsymbol{S}_9\) odabire se tako da se u kombinaciji s tom matricom član geometrijske krutosti izrazi kao \(\delta\boldsymbol{u}^{eT}\, \boldsymbol{F}_9^T \boldsymbol{S}_9 \boldsymbol{F}_9\, \dot{\boldsymbol{u}}^e\); to je sljedeća matrica dimenzija \(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} \]

Dobiva se ta matrica.

Updated Lagrange formulacija

U Updated Lagrange formulaciji pretpostavlja se linearna veza između Jaumannove stope relativnog Kirchhoffova tenzora naprezanja \(\hat{\boldsymbol{\sigma}}^{\nabla J}\) i tenzora brzine deformacije \(\boldsymbol{D}\), odnosno \(\hat{\boldsymbol{\sigma}}^{\nabla J} = \boldsymbol{\mathsf{C}}:\boldsymbol{D}\). To je oblik hipoelastičnog konstitutivnog zakona uobičajen za linearno elastične, elastoplastične i puzajuće materijale, a FrontISTR za te materijale koristi Updated Lagrange formulaciju. Integrand tangentne krutosti elementa izražen u trenutačnoj konfiguraciji tada se u tenzorskom obliku zapisuje kao

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

gdje je \(\boldsymbol{\sigma}^{\nabla T}\) Truesdellova stopa, \(\boldsymbol{A}_{(L)}\) linearni dio Almansijeve deformacije, \(\boldsymbol{F}_t = \partial\boldsymbol{u}/\partial\boldsymbol{x}\) gradijent pomaka u odnosu na trenutačnu konfiguraciju, a \(\boldsymbol{L}\) tenzor gradijenta brzine. Prvi član na desnoj strani jest član materijalne krutosti, a drugi član geometrijske krutosti.

U implementaciji FrontISTR-a taj se integrand izračunava u matričnom obliku koristeći Voigtovu notaciju:

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

Ovdje je \(\boldsymbol{b}\) B-matrica konstruirana u trenutačnoj konfiguraciji (Diskretizacija virtualnog rada unutarnjih sila). \(\boldsymbol{\sigma}_9, \boldsymbol{f}_9\) dobivaju se iz \(\boldsymbol{S}_9, \boldsymbol{F}_9\) definiranih za TL formulaciju tako da se drugo PK naprezanje \(\boldsymbol{S}\) zamijeni Cauchyjevim naprezanjem \(\boldsymbol{\sigma}\), a gradijent referentne konfiguracije \(\partial N_\alpha^e/\partial X_i\) gradijentom trenutačne konfiguracije \(\partial N_\alpha^e/\partial x_i\).

\(\boldsymbol{G}\) je korekcijska matrica ovisna o Cauchyjevu naprezanju, potrebna kako bi hipoelastični konstitutivni zakon \(\hat{\boldsymbol{\sigma}}^{\nabla J} = \boldsymbol{\mathsf{C}}:\boldsymbol{D}\) bio usklađen s okvirom tangentne krutosti kao konstitutivni zakon temeljen na Truesdellovoj stopi. Dobiva se slaganjem komponenti tenzora četvrtog reda \(G_{ijkl} = \delta_{il}\sigma_{kj} + \delta_{jl}\sigma_{ik}\) u \(6\times 6\) Voigtov oblik kao

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

Dobiva se ta matrica.

Sastavljanje globalne matrice krutosti

Globalna tangentna krutost \(\boldsymbol{K}\) dobiva se dijeljenjem svake krutosti elementa \(\boldsymbol{K}^e\) na blokove \(d\times d\) \(\boldsymbol{K}^e_{\alpha\beta}\) za svaki par čvorova i uporabom skupa za sastavljanje tenzora drugog reda \(\mathcal{E}^2(i_g, i_h)\) uvedenog u Sastavljanju fizikalnih veličina na čvorovima elemenata:

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

Dobivene vrijednosti raspoređuju se u matricu s retkom \(i_g\) i stupcem \(i_h\). U implementaciji se skup \(\mathcal{E}^2\) ne konstruira izričito; umjesto toga, odgovarajući se blokovi izravno pribrajaju unutar petlje po elementima. Matrica je kvadratna, dimenzije jednake broju stupnjeva slobode po čvoru \(\times\) ukupnom broju čvorova \(n_g\), ali budući da su komponente osim onih između čvorova povezanih elementima jednake \(0\), pohranjuje se u obliku rijetke matrice.

Matrice krutosti elemenata za TL i UL formulacije imaju isti oblik osim zamjene referentne konfiguracije (koordinate čvorova i izvor za konstrukciju B-matrice) te prisutnosti ili odsutnosti matrice \(\boldsymbol{G}\). FrontISTR zato obje formulacije implementira u zajedničkoj podrutini.

Iteracijski algoritam

Sažimajući prethodno, na početku iteracije postavi se \(\Delta\boldsymbol{u} = \boldsymbol{0}\) i izračuna početni rezidual \(\boldsymbol{R}_0 = \boldsymbol{F}_{n+1} - \boldsymbol{Q}(\boldsymbol{u}_n)\). Zatim se u \(i\)-toj iteraciji provodi sljedeći postupak.

  1. Pri trenutačnom pomaku \(\boldsymbol{u}_n + \Delta\boldsymbol{u}\) izračunaj tangentnu krutost \(\boldsymbol{K}_i = \partial\boldsymbol{Q}/\partial\boldsymbol{u}|_{\boldsymbol{u}_n + \Delta\boldsymbol{u}}\) postupkom opisanim u Konstrukciji tangentne matrice krutosti.
  2. Za nametanje geometrijskih graničnih uvjeta izmijeni matricu tangentne krutosti i vektor reziduala za stupnjeve slobode na koje se primjenjuju ograničenja pomaka, čime se dobivaju \(\tilde{\boldsymbol{K}}_i, \tilde{\boldsymbol{R}}_{i-1}\) (Obrada geometrijskih graničnih uvjeta).
  3. Riješi linearnu jednadžbu \(\tilde{\boldsymbol{K}}_i\, d\boldsymbol{u}_i = \tilde{\boldsymbol{R}}_{i-1}\) i dobije se korekcija \(d\boldsymbol{u}_i\). Taj postupak često čini najveći dio računalnog troška iteracijskog izračuna.
  4. Ažuriraj prirast pomaka kao \(\Delta\boldsymbol{u} \leftarrow \Delta\boldsymbol{u} + d\boldsymbol{u}_i\), a zatim izračunaj vektor unutarnjih sila \(\boldsymbol{Q}(\boldsymbol{u}_n + \Delta\boldsymbol{u})\) i rezidual \(\boldsymbol{R}_i = \boldsymbol{F}_{n+1} - \boldsymbol{Q}(\boldsymbol{u}_n + \Delta\boldsymbol{u})\).
  5. Provjeri konvergenciju i završi iteraciju ako je konvergencija postignuta. Komponente koje odgovaraju reakcijama na ograničenjima pojavljuju se u rezidualu \(\boldsymbol{R}_i\) na stupnjevima slobode na koje se primjenjuju geometrijski granični uvjeti, pa se pokazatelj konvergencije konstruira iz \(\tilde{\boldsymbol{R}}_i\) nakon izuzimanja tih komponenti. Konkretni pokazatelji konvergencije i pragovi opisani su u Kriterijima konvergencije. Ako konvergencija nije postignuta i dosegne se ograničenje broja iteracija, iteracija se smatra neuspjelom.

Kada iteracija konvergira, konvergirani \(\Delta\boldsymbol{u}\) pribraja se \(\boldsymbol{u}_n\) kako bi se dobio akumulirani pomak u trenutku \(t_{n+1}\), \(\boldsymbol{u}_{n+1} = \boldsymbol{u}_n + \Delta\boldsymbol{u}\), te se prelazi na sljedeći vremenski korak.

Povezane teme

AI-assisted translation May contain errors Official docs Status