Přeskočit obsah

Newtonova-Raphsonova metoda

Linearizace a iterační rekurence

Virtuální práce vnějších sil a sestavení globální rovnice poskytuje nelineární rovnici pro uzlové posunutí v čase \(t_{n+1}\), \(\boldsymbol{u}_{n+1}\), která se řeší Newtonovou-Raphsonovou metodou. Uzlové posunutí do času \(t_n\), \(\boldsymbol{u}_n\), se považuje za známé a přírůstek posunutí \(\Delta\boldsymbol{u}\) se bere jako neznámá proměnná, kterou je třeba určit

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

Dále zanedbáme závislost vektoru vnějších sil na uzlovém posunutí a položíme \(\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 řeší.

V aktuálním řešení \(\Delta\boldsymbol{u}\) definujme tečnou tuhost

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

Pomocí ní linearizace nelineární rovnice dává

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

Nechť korekce v \(i\)-té iteraci je \(d\boldsymbol{u}_i\) a vektor rezidua na začátku iterace je

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

Potom má iterační rekurence tvar

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

Reziduum \(\boldsymbol{R}_i\) tedy představuje veličinu odpovídající nerovnováze sil vůči rovnovážnému stavu.

Konstrukce tečné matice tuhosti

Tečná tuhost \(\boldsymbol{K} = \partial\boldsymbol{Q}/\partial\boldsymbol{u}\) se sestaví parciální derivací vektoru vnitřních sil prvku získaného v části Diskretizace virtuální práce vnitřních sil podle uzlového posunutí, integrací vzniklých integrandů na úrovni prvku přes každou oblast prvku a jejich sestavením. Označíme-li integrand na úrovni prvku jako \(\boldsymbol{K}^e_X\) (zápis v referenční konfiguraci, formulace TL) nebo \(\boldsymbol{K}^e_x\) (zápis v aktuální konfiguraci, formulace UL), je tečná tuhost prvku

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

Níže jsou uvedeny výsledné tvary integrandů TL/UL. V obou případech se rozkládají na součet členu materiálové tuhosti (člen počátečního posunutí) a členu geometrické tuhosti (člen počátečního napětí).

Formulace Total Lagrange

Ve formulaci Total Lagrange se předpokládá lineární vztah mezi rychlostí druhého Piolova-Kirchhoffova napětí \(\dot{\boldsymbol{S}}\) a rychlostí Greenova-Lagrangeova přetvoření \(\dot{\boldsymbol{E}}\), tedy \(\dot{\boldsymbol{S}} = \boldsymbol{\mathsf{C}}:\dot{\boldsymbol{E}}\). To odpovídá konstitutivním zákonům pro lineárně elastické materiály (St. Venantovy-Kirchhoffovy materiály) a hyperelastické materiály a FrontISTR pro tyto materiály používá formulaci Total Lagrange. Integrand tečné tuhosti prvku lze potom v tenzorovém tvaru zapsat jako

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

První člen na pravé straně je člen materiálové tuhosti (člen počátečního posunutí) a druhý člen je člen geometrické tuhosti (člen počátečního napětí).

V implementaci FrontISTR se tento integrand vyhodnocuje v maticovém tvaru pomocí Voigtova 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 jsou následující. \(\boldsymbol{B}_L, \boldsymbol{B}_{NL}\) jsou B-matice zavedené v části Diskretizace virtuální práce vnitřních sil a \(\tilde{\boldsymbol{C}}\) je Voigtova reprezentace konstitutivního tenzoru \(\boldsymbol{\mathsf{C}}\), tedy materiálová matice tuhosti \(6\times 6\) (Tenzorový zápis a matematické základy). \(\boldsymbol{S}_9, \boldsymbol{F}_9\) jsou následující přeskupovací matice používané k vyjádření členu geometrické tuhosti jako maticového součinu. Nejprve pro tenzor druhého řádu \(\boldsymbol{A}\) o rozměru \(3\times 3\) definujeme zápis \([\,\cdot\,]\), který jej přeskupí do devítisložkového vektoru, jako

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

S touto definicí \(\boldsymbol{F}_9\) vyjadřuje variaci deformačního gradientu ve tvaru \([\delta\boldsymbol{F}] = \boldsymbol{F}_9\, \delta\boldsymbol{u}^e\) a je maticí \(9\times d n_e\). Pro uzel prvku \(\alpha = 1, \ldots, n_e\) má odpovídající blok \(9\times d\) tvar

\[ [\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á matice } 3\times 3) \]

a platí \(\boldsymbol{F}_9 = [[\boldsymbol{F}_9]_1, \ldots, [\boldsymbol{F}_9]_{n_e}]\), přičemž bloky jsou uspořádány vodorovně podle pořadí uzlů prvku. \(\boldsymbol{S}_9\) je zvolena tak, aby se při kombinaci s touto maticí člen geometrické tuhosti vyjádřil jako \(\delta\boldsymbol{u}^{eT}\, \boldsymbol{F}_9^T \boldsymbol{S}_9 \boldsymbol{F}_9\, \dot{\boldsymbol{u}}^e\); jde o následující matici \(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} \]

Tím získáme výslednou matici.

Formulace Updated Lagrange

Ve formulaci Updated Lagrange se předpokládá lineární vztah mezi Jaumannovou rychlostí relativního Kirchhoffova tenzoru napětí \(\hat{\boldsymbol{\sigma}}^{\nabla J}\) a tenzorem rychlosti deformace \(\boldsymbol{D}\), tedy \(\hat{\boldsymbol{\sigma}}^{\nabla J} = \boldsymbol{\mathsf{C}}:\boldsymbol{D}\). Jde o tvar hypoelastického konstitutivního zákona běžný pro lineárně elastické, elastoplastické a creepové materiály a FrontISTR pro tyto materiály používá formulaci Updated Lagrange. Integrand tečné tuhosti prvku vyjádřený v aktuální konfiguraci lze potom v tenzorovém tvaru zapsat jako

\[ \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 rychlost, \(\boldsymbol{A}_{(L)}\) je lineární část Almansiho přetvoření, \(\boldsymbol{F}_t = \partial\boldsymbol{u}/\partial\boldsymbol{x}\) je gradient posunutí vzhledem k aktuální konfiguraci a \(\boldsymbol{L}\) je tenzor gradientu rychlosti. První člen na pravé straně je člen materiálové tuhosti a druhý člen je člen geometrické tuhosti.

V implementaci FrontISTR se tento integrand vyhodnocuje v maticovém tvaru pomocí Voigtova 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 \]

Zde je \(\boldsymbol{b}\) B-matice sestavená v aktuální konfiguraci (Diskretizace virtuální práce vnitřních sil). \(\boldsymbol{\sigma}_9, \boldsymbol{f}_9\) se získají z \(\boldsymbol{S}_9, \boldsymbol{F}_9\) definovaných pro formulaci TL nahrazením druhého PK napětí \(\boldsymbol{S}\) Cauchyho napětím \(\boldsymbol{\sigma}\) a gradientu v referenční konfiguraci \(\partial N_\alpha^e/\partial X_i\) gradientem v aktuální konfiguraci \(\partial N_\alpha^e/\partial x_i\).

\(\boldsymbol{G}\) je korekční matice závislá na Cauchyho napětí, která je potřebná k uvedení hypoelastického konstitutivního zákona \(\hat{\boldsymbol{\sigma}}^{\nabla J} = \boldsymbol{\mathsf{C}}:\boldsymbol{D}\) do souladu s rámcem tečné tuhosti jako konstitutivního zákona založeného na Truesdellově rychlosti. Získá se uspořádáním složek tenzoru čtvrtého řádu \(G_{ijkl} = \delta_{il}\sigma_{kj} + \delta_{jl}\sigma_{ik}\) do Voigtova tvaru \(6\times 6\) jako

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

Tím získáme výslednou matici.

Sestavení globální matice tuhosti

Globální tečná tuhost \(\boldsymbol{K}\) se získá rozdělením každé tuhosti prvku \(\boldsymbol{K}^e\) na bloky \(d\times d\) \(\boldsymbol{K}^e_{\alpha\beta}\) pro každou dvojici uzlů a použitím množiny sestavení tenzoru druhého řádu \(\mathcal{E}^2(i_g, i_h)\) zavedené v části Sestavení fyzikálních veličin v uzlech prvků:

\[ \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 se uspořádají do matice s řádkem \(i_g\) a sloupcem \(i_h\). V implementaci se množina \(\mathcal{E}^2\) explicitně nekonstruuje; odpovídající bloky se místo toho přímo přičítají uvnitř smyčky přes prvky. Matice je čtvercová s rozměrem rovným počtu stupňů volnosti na uzel \(\times\) celkovému počtu uzlů \(n_g\), avšak protože složky jiné než mezi uzly propojenými prvky jsou \(0\), ukládá se ve formátu řídké matice.

Matice tuhosti prvků pro formulace TL a UL mají stejný tvar s výjimkou přepnutí referenční konfigurace (uzlové souřadnice a zdroj použitý ke konstrukci B-matice) a přítomnosti či nepřítomnosti matice \(\boldsymbol{G}\). FrontISTR proto implementuje obě formulace ve společném podprogramu.

Iterační algoritmus

Shrneme-li výše uvedené, na začátku iterace nastavíme \(\Delta\boldsymbol{u} = \boldsymbol{0}\) a vypočteme počáteční reziduum \(\boldsymbol{R}_0 = \boldsymbol{F}_{n+1} - \boldsymbol{Q}(\boldsymbol{u}_n)\). Poté v \(i\)-té iteraci provedeme následující postup.

  1. Při aktuálním posunutí \(\boldsymbol{u}_n + \Delta\boldsymbol{u}\) vypočtěte tečnou tuhost \(\boldsymbol{K}_i = \partial\boldsymbol{Q}/\partial\boldsymbol{u}|_{\boldsymbol{u}_n + \Delta\boldsymbol{u}}\) postupem uvedeným v části Konstrukce tečné matice tuhosti.
  2. Pro uložení geometrických okrajových podmínek upravte tečnou matici tuhosti a vektor rezidua pro stupně volnosti podléhající omezením posunutí a získejte \(\tilde{\boldsymbol{K}}_i, \tilde{\boldsymbol{R}}_{i-1}\) (Zpracování geometrických okrajových podmínek).
  3. Vyřešte lineární rovnici \(\tilde{\boldsymbol{K}}_i\, d\boldsymbol{u}_i = \tilde{\boldsymbol{R}}_{i-1}\) a získejte korekci \(d\boldsymbol{u}_i\). Tento postup často představuje většinu výpočetních nákladů iteračního výpočtu.
  4. Aktualizujte přírůstek posunutí jako \(\Delta\boldsymbol{u} \leftarrow \Delta\boldsymbol{u} + d\boldsymbol{u}_i\) a následně vypočtěte vektor vnitřních sil \(\boldsymbol{Q}(\boldsymbol{u}_n + \Delta\boldsymbol{u})\) a reziduum \(\boldsymbol{R}_i = \boldsymbol{F}_{n+1} - \boldsymbol{Q}(\boldsymbol{u}_n + \Delta\boldsymbol{u})\).
  5. Zkontrolujte konvergenci a při dosažení konvergence iteraci ukončete. Ve stupních volnosti podléhajících geometrickým okrajovým podmínkám se v reziduu \(\boldsymbol{R}_i\) objevují složky odpovídající reakcím omezení, proto se ukazatel konvergence sestavuje z \(\tilde{\boldsymbol{R}}_i\) po vyloučení těchto složek. Konkrétní ukazatele konvergence a prahové hodnoty jsou popsány v části Kritéria konvergence. Pokud konvergence není dosažena a je dosažen limit iterací, iterace se považuje za neúspěšnou.

Po konvergenci iterace přičtěte konvergovaný \(\Delta\boldsymbol{u}\) k \(\boldsymbol{u}_n\), čímž získáte akumulované posunutí v čase \(t_{n+1}\), \(\boldsymbol{u}_{n+1} = \boldsymbol{u}_n + \Delta\boldsymbol{u}\), a pokračujte k dalšímu časovému kroku.

Související témata

AI-assisted translation May contain errors Official docs Status