Preskočiť na obsah

Newtonova-Raphsonova metóda

Linearizácia a iteračná rekurencia

Virtuálna práca vonkajších síl a zostavenie globálnej rovnice dáva nelineárnu rovnicu pre uzlové posunutie v čase \(t_{n+1}\), \(\boldsymbol{u}_{n+1}\), ktorá sa rieši Newtonovou-Raphsonovou metódou. Uzlové posunutie do času \(t_n\), \(\boldsymbol{u}_n\), sa považuje za známe a prírastok posunutia \(\Delta\boldsymbol{u}\) sa berie ako neznáma premenná, ktorú treba určiť

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

Ďalej sa zanedbáva závislosť vektora vonkajších síl od uzlového posunutia a pri \(\boldsymbol{F}(\boldsymbol{u}_{n+1}) = \boldsymbol{F}_{n+1}\)

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

sa rieši.

Pri aktuálnom riešení \(\Delta\boldsymbol{u}\) definujme dotykovú tuhosť

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

Použitím tejto veličiny dá linearizácia nelineárnej rovnice

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

Nech korekcia v \(i\)-tej iterácii je \(d\boldsymbol{u}_i\) a vektor rezíduí na začiatku iterácie je

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

Potom je iteračná rekurencia

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

Rezíduum \(\boldsymbol{R}_i\) je teda veličina zodpovedajúca nevyváženosti síl vzhľadom na rovnováhu.

Zostavenie dotykovej matice tuhosti

Dotyková tuhosť \(\boldsymbol{K} = \partial\boldsymbol{Q}/\partial\boldsymbol{u}\) sa zostaví parciálnym derivovaním vektora vnútorných síl prvku získaného v Diskretizácii virtuálnej práce vnútorných síl podľa uzlového posunutia, integráciou výsledných integrandov na úrovni prvku cez každú oblasť prvku a ich zostavením. Ak označíme integrand na úrovni prvku ako \(\boldsymbol{K}^e_X\) (zápis v referenčnej konfigurácii, formulácia TL) alebo \(\boldsymbol{K}^e_x\) (zápis v aktuálnej konfigurácii, formulácia UL), dotyková tuhosť prvku 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}) \]

Nasledujú konečné tvary integrandov TL/UL. V oboch prípadoch sa rozkladajú na súčet člena materiálovej tuhosti (člen počiatočného posunutia) a člena geometrickej tuhosti (člen počiatočného napätia).

Úplná Lagrangeova formulácia

V úplnej Lagrangeovej formulácii sa predpokladá lineárny vzťah medzi rýchlosťou druhého Piolovho-Kirchhoffovho napätia \(\dot{\boldsymbol{S}}\) a rýchlosťou Greenovej-Lagrangeovej deformácie \(\dot{\boldsymbol{E}}\), konkrétne \(\dot{\boldsymbol{S}} = \boldsymbol{\mathsf{C}}:\dot{\boldsymbol{E}}\). To zodpovedá konštitutívnym zákonom lineárne elastických materiálov (materiály St. Venant-Kirchhoff) a hyperelastických materiálov a FrontISTR pre tieto materiály používa úplnú Lagrangeovu formuláciu. Integrand dotykovej tuhosti prvku sa potom v tenzorovom tvare zapíše ako

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

Prvý člen na pravej strane je člen materiálovej tuhosti (člen počiatočného posunutia) a druhý člen je člen geometrickej tuhosti (člen počiatočného napätia).

V implementácii FrontISTR sa tento integrand vyhodnocuje v maticovom tvare pomocou Voigtovho zápisu:

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

Jednotlivé matice sú nasledujúce. \(\boldsymbol{B}_L, \boldsymbol{B}_{NL}\) sú B matice zavedené v Diskretizácii virtuálnej práce vnútorných síl a \(\tilde{\boldsymbol{C}}\) je Voigtova reprezentácia konštitutívneho tenzora \(\boldsymbol{\mathsf{C}}\), teda materiálová matica tuhosti \(6\times 6\) (Tenzorový zápis a matematické základy). \(\boldsymbol{S}_9, \boldsymbol{F}_9\) sú nasledujúce preusporiadavacie matice používané na vyjadrenie člena geometrickej tuhosti ako maticového súčinu. Najprv pre tenzor druhého rádu \(3\times 3\) \(\boldsymbol{A}\) definujme zápis \([\,\cdot\,]\), ktorý ho preusporiada na 9-zložkový vektor, ako

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

Pri tejto definícii \(\boldsymbol{F}_9\) vyjadruje variáciu gradientu deformácie v tvare \([\delta\boldsymbol{F}] = \boldsymbol{F}_9\, \delta\boldsymbol{u}^e\) a je maticou \(9\times d n_e\). Pre uzol prvku \(\alpha = 1, \ldots, n_e\) je zodpovedajúci blok \(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 jednotková matica } 3\times 3) \]

a je daný vzťahom \(\boldsymbol{F}_9 = [[\boldsymbol{F}_9]_1, \ldots, [\boldsymbol{F}_9]_{n_e}]\), pričom bloky sú usporiadané vodorovne v poradí uzlov prvku. \(\boldsymbol{S}_9\) sa volí tak, aby sa v kombinácii s touto maticou člen geometrickej tuhosti vyjadril ako \(\delta\boldsymbol{u}^{eT}\, \boldsymbol{F}_9^T \boldsymbol{S}_9 \boldsymbol{F}_9\, \dot{\boldsymbol{u}}^e\); ide o nasledujúcu maticu \(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} \]

Toto je výsledná matica.

Aktualizovaná Lagrangeova formulácia

V aktualizovanej Lagrangeovej formulácii sa predpokladá lineárny vzťah medzi Jaumannovou rýchlosťou relatívneho Kirchhoffovho tenzora napätia \(\hat{\boldsymbol{\sigma}}^{\nabla J}\) a tenzorom rýchlosti deformácie \(\boldsymbol{D}\), konkrétne \(\hat{\boldsymbol{\sigma}}^{\nabla J} = \boldsymbol{\mathsf{C}}:\boldsymbol{D}\). Ide o tvar hypoelastického konštitutívneho zákona spoločného pre lineárne elastické, elastoplastické a creepové materiály a FrontISTR pre tieto materiály používa aktualizovanú Lagrangeovu formuláciu. Integrand dotykovej tuhosti prvku vyjadrený v aktuálnej konfigurácii sa potom v tenzorovom tvare zapíše ako

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

kde \(\boldsymbol{\sigma}^{\nabla T}\) je Truesdellova rýchlosť, \(\boldsymbol{A}_{(L)}\) je lineárna časť Almansiho deformácie, \(\boldsymbol{F}_t = \partial\boldsymbol{u}/\partial\boldsymbol{x}\) je gradient posunutia vzhľadom na aktuálnu konfiguráciu a \(\boldsymbol{L}\) je tenzor gradientu rýchlosti. Prvý člen na pravej strane je člen materiálovej tuhosti a druhý člen je člen geometrickej tuhosti.

V implementácii FrontISTR sa tento integrand vyhodnocuje v maticovom tvare pomocou Voigtovho zápisu:

\[ \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}\) B matica zostavená v aktuálnej konfigurácii (Diskretizácia virtuálnej práce vnútorných síl). \(\boldsymbol{\sigma}_9, \boldsymbol{f}_9\) sa získajú z \(\boldsymbol{S}_9, \boldsymbol{F}_9\) definovaných pre formuláciu TL nahradením druhého PK napätia \(\boldsymbol{S}\) Cauchyho napätím \(\boldsymbol{\sigma}\) a gradientu v referenčnej konfigurácii \(\partial N_\alpha^e/\partial X_i\) gradientom v aktuálnej konfigurácii \(\partial N_\alpha^e/\partial x_i\).

\(\boldsymbol{G}\) je korekčná matica závislá od Cauchyho napätia, potrebná na zosúladenie hypoelastického konštitutívneho zákona \(\hat{\boldsymbol{\sigma}}^{\nabla J} = \boldsymbol{\mathsf{C}}:\boldsymbol{D}\) s rámcom dotykovej tuhosti ako konštitutívneho zákona založeného na Truesdellovej rýchlosti. Získa sa usporiadaním zložiek tenzora štvrtého rádu \(G_{ijkl} = \delta_{il}\sigma_{kj} + \delta_{jl}\sigma_{ik}\) do Voigtovho tvaru \(6\times 6\) ako

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

Toto je výsledná matica.

Zostavenie globálnej matice tuhosti

Globálna dotyková tuhosť \(\boldsymbol{K}\) sa získa rozdelením každej tuhosti prvku \(\boldsymbol{K}^e\) na bloky \(d\times d\) \(\boldsymbol{K}^e_{\alpha\beta}\) pre každý pár uzlov a použitím množiny zostavenia tenzora druhého rádu \(\mathcal{E}^2(i_g, i_h)\) zavedenej v Zostavení fyzikálnych veličín uzlov prvkov:

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

Výsledné hodnoty sa usporiadajú do matice s riadkom \(i_g\) a stĺpcom \(i_h\). V implementácii sa množina \(\mathcal{E}^2\) explicitne nevytvára; namiesto toho sa príslušné bloky pridávajú priamo v cykle prvkov. Matica je štvorcová s rozmerom počet stupňov voľnosti na uzol \(\times\) celkový počet uzlov \(n_g\), ale keďže zložky mimo uzlov spojených prvkami sú \(0\), ukladá sa v riedkom tvare.

Matice tuhosti prvkov pre formulácie TL a UL majú rovnaký tvar s výnimkou zmeny referenčnej konfigurácie (uzlových súradníc a zdroja použitého na zostavenie B matice) a prítomnosti alebo neprítomnosti matice \(\boldsymbol{G}\). FrontISTR preto implementuje obe formulácie v spoločnom podprograme.

Iteračný algoritmus

Zhrnutím vyššie uvedeného: na začiatku iterácie nastavte \(\Delta\boldsymbol{u} = \boldsymbol{0}\) a vypočítajte počiatočné rezíduum \(\boldsymbol{R}_0 = \boldsymbol{F}_{n+1} - \boldsymbol{Q}(\boldsymbol{u}_n)\). Potom v \(i\)-tej iterácii vykonajte nasledujúci postup.

  1. Pri aktuálnom posunutí \(\boldsymbol{u}_n + \Delta\boldsymbol{u}\) vypočítajte dotykovú tuhosť \(\boldsymbol{K}_i = \partial\boldsymbol{Q}/\partial\boldsymbol{u}|_{\boldsymbol{u}_n + \Delta\boldsymbol{u}}\) postupom uvedeným v Zostavení dotykovej matice tuhosti.
  2. Na vynútenie geometrických okrajových podmienok upravte dotykovú maticu tuhosti a vektor rezíduí pre stupne voľnosti podliehajúce obmedzeniam posunutia, čím získate \(\tilde{\boldsymbol{K}}_i, \tilde{\boldsymbol{R}}_{i-1}\) (Spracovanie geometrických okrajových podmienok).
  3. Vyriešte lineárnu rovnicu \(\tilde{\boldsymbol{K}}_i\, d\boldsymbol{u}_i = \tilde{\boldsymbol{R}}_{i-1}\) a získajte korekciu \(d\boldsymbol{u}_i\). Tento postup často tvorí väčšinu výpočtových nákladov iteračného výpočtu.
  4. Aktualizujte prírastok posunutia ako \(\Delta\boldsymbol{u} \leftarrow \Delta\boldsymbol{u} + d\boldsymbol{u}_i\) a zodpovedajúcim spôsobom vypočítajte vektor vnútorných síl \(\boldsymbol{Q}(\boldsymbol{u}_n + \Delta\boldsymbol{u})\) a rezíduum \(\boldsymbol{R}_i = \boldsymbol{F}_{n+1} - \boldsymbol{Q}(\boldsymbol{u}_n + \Delta\boldsymbol{u})\).
  5. Skontrolujte konvergenciu a po jej dosiahnutí iteráciu ukončite. V rezíduu \(\boldsymbol{R}_i\) sa na stupňoch voľnosti podliehajúcich geometrickým okrajovým podmienkam objavujú zložky zodpovedajúce reakciám väzieb, preto sa indikátor konvergencie zostavuje z \(\tilde{\boldsymbol{R}}_i\) po vylúčení týchto zložiek. Konkrétne indikátory konvergencie a prahy sú opísané v Kritériách konvergencie. Ak sa konvergencia nedosiahne a dosiahne sa limit iterácií, iterácia sa považuje za neúspešnú.

Keď iterácia skonverguje, pripočítajte skonvergované \(\Delta\boldsymbol{u}\) k \(\boldsymbol{u}_n\), čím získate akumulované posunutie v čase \(t_{n+1}\), \(\boldsymbol{u}_{n+1} = \boldsymbol{u}_n + \Delta\boldsymbol{u}\), a pokračujte ďalším časovým krokom.

Súvisiace témy

AI-assisted translation May contain errors Official docs Status