Kihagyás

Newton-Raphson-módszer

Linearizálás és iterációs rekurzió

A Külső erők virtuális munkája és a globális egyenlet összeállítása fejezetben kapott, a \(t_{n+1}\) időpontbeli \(\boldsymbol{u}_{n+1}\) csomóponti elmozdulásra vonatkozó nemlineáris egyenletet Newton-Raphson-módszerrel oldjuk meg. A \(t_n\) időpontig ismertnek tekintjük a \(\boldsymbol{u}_n\) csomóponti elmozdulást, és a \(\Delta\boldsymbol{u}\) elmozdulásnövekményt ismeretlen változóként határozzuk meg úgy, hogy

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

teljesüljön. A továbbiakban elhanyagoljuk a külsőerő-vektor csomóponti elmozdulástól való függését, és \(\boldsymbol{F}(\boldsymbol{u}_{n+1}) = \boldsymbol{F}_{n+1}\) mellett a

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

egyenletet oldjuk meg.

A jelenlegi \(\Delta\boldsymbol{u}\) megoldásnál az érintőmerevséget

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

definiálva a nemlineáris egyenlet linearizált alakja

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

lesz. Jelölje az \(i\)-edik iteráció korrekcióját \(d\boldsymbol{u}_i\), az iteráció kezdetén vett maradékvektort pedig

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

Ekkor az iterációs rekurzió

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

alakú. A \(\boldsymbol{R}_i\) maradék a mechanikai egyensúlytól való erőeltérést jellemző mennyiség.

Az érintőmerevségi mátrix felépítése

Az \(\boldsymbol{K} = \partial\boldsymbol{Q}/\partial\boldsymbol{u}\) érintőmerevséget úgy képezzük, hogy a Belső erők virtuális munkájának diszkretizálása fejezetben kapott elem-belsőerővektort a csomóponti elmozdulás szerint parciálisan differenciáljuk, az így kapott elemszintű integrandusokat az egyes elemtartományokon integráljuk, majd összeállítjuk. Az elemszintű integrandust \(\boldsymbol{K}^e_X\)-szel (referencia-konfigurációs jelölés, TL-módszer) vagy \(\boldsymbol{K}^e_x\)-szel (aktuális konfigurációs jelölés, UL-módszer) jelölve az elem érintőmerevsége

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

alakban adódik. Az alábbiakban a TL/UL integrandusok végső alakját adjuk meg. Mindkettő egy anyagmerevségi tag (kezdeti elmozdulási tag) és egy geometriai merevségi tag (kezdeti feszültségi tag) összegeként bontható fel.

Total Lagrange-módszer

A Total Lagrange-módszerben a második Piola-Kirchhoff-feszültség sebessége \(\dot{\boldsymbol{S}}\) és a Green-Lagrange-alakváltozás sebessége \(\dot{\boldsymbol{E}}\) között a \(\dot{\boldsymbol{S}} = \boldsymbol{\mathsf{C}}:\dot{\boldsymbol{E}}\) lineáris kapcsolatot feltételezzük. Ez a lineárisan rugalmas anyagok (St. Venant-Kirchhoff-anyagok) és a hiperelasztikus anyagok konstitutív törvényeinek felel meg; a FrontISTR ezekhez az anyagokhoz a Total Lagrange-módszert használja. Ekkor az elem érintőmerevségének integrandusa tenzoralakban

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

alakú. A jobb oldal első tagja az anyagmerevségi tag (kezdeti elmozdulási tag), a második tag pedig a geometriai merevségi tag (kezdeti feszültségi tag).

A FrontISTR megvalósításában ezt az integrandust Voigt-jelölésű mátrixalakban számítjuk:

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

Az egyes mátrixok a következők. A \(\boldsymbol{B}_L, \boldsymbol{B}_{NL}\) a Belső erők virtuális munkájának diszkretizálása fejezetben bevezetett B-mátrixok, a \(\tilde{\boldsymbol{C}}\) pedig a \(\boldsymbol{\mathsf{C}}\) konstitutív tenzor Voigt-reprezentációja, azaz egy \(6\times 6\) anyagmerevségi mátrix (Tenzorjelölés és matematikai alapok). A \(\boldsymbol{S}_9, \boldsymbol{F}_9\) az alábbi átrendező mátrixok, amelyekkel a geometriai merevségi tag mátrixszorzatként írható fel. Először egy \(3\times 3\) másodrendű \(\boldsymbol{A}\) tenzor 9 komponensű vektorrá rendezésére vezessük be a \([\,\cdot\,]\) jelölést:

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

Ezzel a \(\boldsymbol{F}_9\) a deformációs gradiens variációját \([\delta\boldsymbol{F}] = \boldsymbol{F}_9\, \delta\boldsymbol{u}^e\) alakban kifejező \(9\times d n_e\) mátrix. Az \(\alpha = 1, \ldots, n_e\) elemcsomóponthoz tartozó \(9\times d\) blokk

\[ [\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{ a } 3\times 3 \text{ egységmátrix}) \]

és ezeknek az elemcsomópontok sorrendjében vízszintesen egymás mellé rendezett alakja \(\boldsymbol{F}_9 = [[\boldsymbol{F}_9]_1, \ldots, [\boldsymbol{F}_9]_{n_e}]\). A \(\boldsymbol{S}_9\) úgy van megválasztva, hogy ezzel együtt a geometriai merevségi tag \(\delta\boldsymbol{u}^{eT}\, \boldsymbol{F}_9^T \boldsymbol{S}_9 \boldsymbol{F}_9\, \dot{\boldsymbol{u}}^e\) alakban legyen kifejezhető; ez az alábbi \(9\times 9\) mátrix:

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

Ez adja a mátrixot.

Updated Lagrange-módszer

Az Updated Lagrange-módszerben a relatív Kirchhoff-feszültségtenzor Jaumann-sebessége \(\hat{\boldsymbol{\sigma}}^{\nabla J}\) és a deformációsebesség-tenzor \(\boldsymbol{D}\) között a \(\hat{\boldsymbol{\sigma}}^{\nabla J} = \boldsymbol{\mathsf{C}}:\boldsymbol{D}\) lineáris kapcsolatot feltételezzük. Ez a lineárisan rugalmas, elasztoplasztikus és kúszó anyagokra közös hipoelasztikus konstitutív törvény alakja; a FrontISTR ezekhez az anyagokhoz az Updated Lagrange-módszert használja. Ekkor az aktuális konfigurációban kifejezett elem-érintőmerevség integrandusa tenzoralakban

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

alakú (\(\boldsymbol{\sigma}^{\nabla T}\) a Truesdell-sebesség, \(\boldsymbol{A}_{(L)}\) az Almansi-alakváltozás lineáris része, \(\boldsymbol{F}_t = \partial\boldsymbol{u}/\partial\boldsymbol{x}\) az aktuális konfigurációhoz viszonyított elmozdulásgradiens, \(\boldsymbol{L}\) pedig a sebességgradiens-tenzor). A jobb oldal első tagja az anyagmerevségi tag, a második a geometriai merevségi tag.

A FrontISTR megvalósításában ezt az integrandust Voigt-jelölésű mátrixalakban számítjuk:

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

Itt \(\boldsymbol{b}\) az aktuális konfigurációban felépített B-mátrix (Belső erők virtuális munkájának diszkretizálása). A \(\boldsymbol{\sigma}_9, \boldsymbol{f}_9\) a TL-módszernél definiált \(\boldsymbol{S}_9, \boldsymbol{F}_9\) mátrixokból adódik úgy, hogy a második PK-feszültséget \(\boldsymbol{S}\) Cauchy-feszültségre \(\boldsymbol{\sigma}\), a referencia-konfigurációbeli \(\partial N_\alpha^e/\partial X_i\) gradienst pedig az aktuális konfigurációbeli \(\partial N_\alpha^e/\partial x_i\) gradiensre cseréljük.

A \(\boldsymbol{G}\) egy Cauchy-feszültségtől függő korrekciós mátrix, amely ahhoz szükséges, hogy a \(\hat{\boldsymbol{\sigma}}^{\nabla J} = \boldsymbol{\mathsf{C}}:\boldsymbol{D}\) hipoelasztikus konstitutív törvény összhangba kerüljön a Truesdell-sebességen alapuló konstitutív törvényként kezelt érintőmerevségi formulációval. A negyedrendű \(G_{ijkl} = \delta_{il}\sigma_{kj} + \delta_{jl}\sigma_{ik}\) tenzorkomponenseket \(6\times 6\) Voigt-alakba rendezve

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

adódik.

A globális merevségi mátrix összeállítása

A globális \(\boldsymbol{K}\) érintőmerevséget úgy kapjuk, hogy az egyes \(\boldsymbol{K}^e\) elemmerevségeket csomópontpáronként \(d\times d\) méretű \(\boldsymbol{K}^e_{\alpha\beta}\) blokkokra bontjuk, majd a Csomóponti fizikai mennyiségek összeállítása szakaszban bevezetett másodrendű tenzoros \(\mathcal{E}^2(i_g, i_h)\) összeállítási halmaz segítségével

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

értékeit az \(i_g\) sorokba és \(i_h\) oszlopokba rendezzük. A megvalósítás nem hozza létre explicit módon az \(\mathcal{E}^2\) halmazt; ehelyett az elemciklusban közvetlenül hozzáadja az értékeket a megfelelő blokkokhoz. A mátrix négyzetes, mérete a csomópontonkénti szabadságfokok száma \(\times\) a globális csomópontok \(n_g\) száma, de mivel az elemekkel össze nem kapcsolt csomópontok közötti komponensek \(0\) értékűek, ritka mátrixos formában tárolható.

A TL- és UL-módszer elemmerevségi mátrixai a referencia-konfiguráció váltásától (csomóponti koordináták és a B-mátrix felépítésének alapja), valamint a \(\boldsymbol{G}\) mátrix jelenlététől vagy hiányától eltekintve azonos alakúak. A FrontISTR ezért mindkét formulációt közös szubrutinban valósítja meg.

Iterációs algoritmus

A fentieket összefoglalva az iteráció kezdetén \(\Delta\boldsymbol{u} = \boldsymbol{0}\) értéket állítunk be, és kiszámítjuk a kezdeti \(\boldsymbol{R}_0 = \boldsymbol{F}_{n+1} - \boldsymbol{Q}(\boldsymbol{u}_n)\) maradékot. Ezután az \(i\)-edik iterációban a következő lépéseket hajtjuk végre.

  1. Az aktuális \(\boldsymbol{u}_n + \Delta\boldsymbol{u}\) elmozdulásnál a Az érintőmerevségi mátrix felépítése eljárásával kiszámítjuk a \(\boldsymbol{K}_i = \partial\boldsymbol{Q}/\partial\boldsymbol{u}|_{\boldsymbol{u}_n + \Delta\boldsymbol{u}}\) érintőmerevséget.
  2. A geometriai peremfeltételek érvényesítéséhez módosítjuk az elmozdulási kényszerekkel terhelt szabadságfokokhoz tartozó érintőmerevségi mátrixot és maradékvektort, így \(\tilde{\boldsymbol{K}}_i, \tilde{\boldsymbol{R}}_{i-1}\) adódik (Geometriai peremfeltételek kezelése).
  3. Megoldjuk a \(\tilde{\boldsymbol{K}}_i\, d\boldsymbol{u}_i = \tilde{\boldsymbol{R}}_{i-1}\) lineáris egyenletet, és meghatározzuk a \(d\boldsymbol{u}_i\) korrekciót. Ez a művelet gyakran az iterációs számítás számítási költségének legnagyobb részét teszi ki.
  4. Az elmozdulásnövekményt \(\Delta\boldsymbol{u} \leftarrow \Delta\boldsymbol{u} + d\boldsymbol{u}_i\) szerint frissítjük, majd ennek megfelelően kiszámítjuk a \(\boldsymbol{Q}(\boldsymbol{u}_n + \Delta\boldsymbol{u})\) belsőerő-vektort és a \(\boldsymbol{R}_i = \boldsymbol{F}_{n+1} - \boldsymbol{Q}(\boldsymbol{u}_n + \Delta\boldsymbol{u})\) maradékot.
  5. Ellenőrizzük a konvergenciát, és ha teljesült, befejezzük az iterációt. A \(\boldsymbol{R}_i\) maradék azon komponensei, amelyek geometriai peremfeltétellel kényszerezett szabadságfokokhoz tartoznak, a kényszerreakcióknak megfelelő értékeket tartalmaznak, ezért ezeket kizárva a \(\tilde{\boldsymbol{R}}_i\) komponenseiből képezzük a konvergenciamutatót. A konkrét kritériumokat és küszöbértékeket a Konvergenciakritériumok fejezet tárgyalja. Ha nincs konvergencia és elérjük az iterációs korlátot, az iteráció sikertelennek minősül.

A konvergencia elérésekor a konvergált \(\Delta\boldsymbol{u}\) értéket hozzáadjuk \(\boldsymbol{u}_n\)-hez, így megkapjuk a \(t_{n+1}\) időpontbeli kumulatív \(\boldsymbol{u}_{n+1} = \boldsymbol{u}_n + \Delta\boldsymbol{u}\) elmozdulást, majd továbblépünk a következő időlépésre.

Kapcsolódó elemek

AI-assisted translation May contain errors Official docs Status