Pereiti prie turinio

Newton-Raphson metodas

Tiesinimas ir iteracinė rekurencija

Virtualaus išorinių jėgų darbo ir globaliosios lygties surinkimo skyriuje gaunama netiesinė lygtis mazgų poslinkiui laiku \(t_{n+1}\), \(\boldsymbol{u}_{n+1}\), kuri sprendžiama Newton-Raphson metodu. Laikoma, kad mazgų poslinkis iki laiko \(t_n\), \(\boldsymbol{u}_n\), yra žinomas, o nustatytinas nežinomasis yra poslinkio prieaugis \(\Delta\boldsymbol{u}\):

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

Toliau nepaisoma išorinės jėgos vektoriaus priklausomybės nuo mazgų poslinkio ir, laikant \(\boldsymbol{F}(\boldsymbol{u}_{n+1}) = \boldsymbol{F}_{n+1}\), sprendžiama lygtis

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

Esamam sprendiniui \(\Delta\boldsymbol{u}\) apibrėžiamas tangentinis standumas

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

Naudojant šį dydį, netiesinės lygties tiesinimas duoda

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

Tegu \(i\)-osios iteracijos pataisa yra \(d\boldsymbol{u}_i\), o iteracijos pradžios liekamosios jėgos vektorius

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

Tuomet iteracinė rekurencija yra

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

Taigi liekamoji jėga \(\boldsymbol{R}_i\) yra dydis, atitinkantis jėgų disbalansą nuo pusiausvyros.

Tangentinio standumo matricos sudarymas

Tangentinis standumas \(\boldsymbol{K} = \partial\boldsymbol{Q}/\partial\boldsymbol{u}\) sudaromas dalinai diferencijuojant Vidinių jėgų virtualaus darbo diskretizavimo skyriuje gautą elemento vidinių jėgų vektorių pagal mazgų poslinkius, integruojant gautus elemento lygio pointegralinius reiškinius kiekvieno elemento srityje ir juos surenkant. Pažymėjus elemento lygio pointegralinį reiškinį \(\boldsymbol{K}^e_X\) (atskaitos konfigūracijos žymėjimas, TL formuluotė) arba \(\boldsymbol{K}^e_x\) (dabartinės konfigūracijos žymėjimas, UL formuluotė), elemento tangentinis standumas yra

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

Toliau pateikiamos galutinės TL ir UL pointegralinių reiškinių formos. Abiem atvejais jos išskaidomos į medžiagos standumo dėmens (pradinio poslinkio dėmens) ir geometrinio standumo dėmens (pradinio įtempio dėmens) sumą.

Total Lagrange formuluotė

Total Lagrange formuluotėje daroma prielaida apie tiesinį ryšį tarp antrojo Piola-Kirchhoff įtempių kitimo greičio \(\dot{\boldsymbol{S}}\) ir Green-Lagrange deformacijų greičio \(\dot{\boldsymbol{E}}\), t. y. \(\dot{\boldsymbol{S}} = \boldsymbol{\mathsf{C}}:\dot{\boldsymbol{E}}\). Tai atitinka tiesiškai elastinių medžiagų (St. Venant-Kirchhoff medžiagų) ir hiperelastinių medžiagų konstitucinius dėsnius; šioms medžiagoms FrontISTR naudoja Total Lagrange formuluotę. Tada elemento tangentinio standumo pointegralinis reiškinys tenzorine forma užrašomas taip:

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

Pirmasis dešiniosios pusės dėmuo yra medžiagos standumo dėmuo (pradinio poslinkio dėmuo), o antrasis – geometrinio standumo dėmuo (pradinio įtempio dėmuo).

FrontISTR realizacijoje šis pointegralinis reiškinys skaičiuojamas matricine forma naudojant Voigt žymėjimą:

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

Matricos apibrėžiamos taip. \(\boldsymbol{B}_L, \boldsymbol{B}_{NL}\) yra B matricos, įvestos Vidinių jėgų virtualaus darbo diskretizavimo skyriuje, o \(\tilde{\boldsymbol{C}}\) yra konstitucinio tenzoriaus \(\boldsymbol{\mathsf{C}}\) Voigt atvaizdas, t. y. \(6\times 6\) medžiagos standumo matrica (Tenzorinis žymėjimas ir matematiniai pagrindai). \(\boldsymbol{S}_9, \boldsymbol{F}_9\) yra toliau pateiktos pertvarkymo matricos, naudojamos geometrinio standumo dėmeniui užrašyti matricos sandauga. Pirmiausia \(3\times 3\) antrojo rango tenzoriui \(\boldsymbol{A}\) apibrėžiamas žymėjimas \([\,\cdot\,]\), pertvarkantis jį į 9 komponentų vektorių:

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

Pagal šį apibrėžimą \(\boldsymbol{F}_9\) išreiškia deformacijų gradiento variaciją forma \([\delta\boldsymbol{F}] = \boldsymbol{F}_9\, \delta\boldsymbol{u}^e\) ir yra \(9\times d n_e\) matrica. Elemento mazgui \(\alpha = 1, \ldots, n_e\) atitinkamas \(9\times d\) blokas yra

\[ [\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{ yra } 3\times 3 \text{ tapatybės matrica}) \]

ir \(\boldsymbol{F}_9 = [[\boldsymbol{F}_9]_1, \ldots, [\boldsymbol{F}_9]_{n_e}]\), kur blokai horizontaliai išdėstyti elementų mazgų tvarka. \(\boldsymbol{S}_9\) parenkama taip, kad kartu su šia matrica geometrinio standumo dėmuo būtų išreikštas kaip \(\delta\boldsymbol{u}^{eT}\, \boldsymbol{F}_9^T \boldsymbol{S}_9 \boldsymbol{F}_9\, \dot{\boldsymbol{u}}^e\); tai yra ši \(9\times 9\) matrica:

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

Tai yra galutinė matrica.

Updated Lagrange formuluotė

Updated Lagrange formuluotėje daroma prielaida apie tiesinį ryšį tarp santykinio Kirchhoff įtempių tenzoriaus Jaumann greičio \(\hat{\boldsymbol{\sigma}}^{\nabla J}\) ir deformacijos greičio tenzoriaus \(\boldsymbol{D}\), t. y. \(\hat{\boldsymbol{\sigma}}^{\nabla J} = \boldsymbol{\mathsf{C}}:\boldsymbol{D}\). Tai yra hipoelastinio konstitucinio dėsnio forma, bendra tiesiškai elastinėms, elastoplastinėms ir valkšnumo medžiagoms; šioms medžiagoms FrontISTR naudoja Updated Lagrange formuluotę. Dabartinėje konfigūracijoje išreikštas elemento tangentinio standumo pointegralinis reiškinys užrašomas tenzorine forma:

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

Čia \(\boldsymbol{\sigma}^{\nabla T}\) yra Truesdell greitis, \(\boldsymbol{A}_{(L)}\) – tiesinė Almansi deformacijos dalis, \(\boldsymbol{F}_t = \partial\boldsymbol{u}/\partial\boldsymbol{x}\) – poslinkio gradientas dabartinės konfigūracijos atžvilgiu, o \(\boldsymbol{L}\) – greičio gradiento tenzorius. Pirmasis dešiniosios pusės dėmuo yra medžiagos standumo dėmuo, o antrasis – geometrinio standumo dėmuo.

FrontISTR realizacijoje šis pointegralinis reiškinys skaičiuojamas matricine forma naudojant Voigt žymėjimą:

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

Čia \(\boldsymbol{b}\) yra B matrica, sudaryta dabartinėje konfigūracijoje (Vidinių jėgų virtualaus darbo diskretizavimas). \(\boldsymbol{\sigma}_9, \boldsymbol{f}_9\) gaunamos iš TL formuluotei apibrėžtų \(\boldsymbol{S}_9, \boldsymbol{F}_9\), pakeičiant antrąjį PK įtempį \(\boldsymbol{S}\) į Cauchy įtempį \(\boldsymbol{\sigma}\) ir atskaitos konfigūracijos gradientą \(\partial N_\alpha^e/\partial X_i\) į dabartinės konfigūracijos gradientą \(\partial N_\alpha^e/\partial x_i\).

\(\boldsymbol{G}\) yra nuo Cauchy įtempių priklausanti korekcijos matrica, reikalinga tam, kad hipoelastinis konstitucinis dėsnis \(\hat{\boldsymbol{\sigma}}^{\nabla J} = \boldsymbol{\mathsf{C}}:\boldsymbol{D}\) būtų suderintas su tangentinio standumo sistema, kurioje konstitucinis dėsnis grindžiamas Truesdell greičiu. Ji gaunama ketvirtojo rango tenzoriaus komponentus \(G_{ijkl} = \delta_{il}\sigma_{kj} + \delta_{jl}\sigma_{ik}\) išdėsčius \(6\times 6\) Voigt forma:

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

Tai yra galutinė matrica.

Globaliosios standumo matricos surinkimas

Globalusis tangentinis standumas \(\boldsymbol{K}\) gaunamas kiekvieno elemento standumą \(\boldsymbol{K}^e\) padalijus į \(d\times d\) blokus \(\boldsymbol{K}^e_{\alpha\beta}\) kiekvienai mazgų porai ir naudojant antrojo rango tenzorinio surinkimo aibę \(\mathcal{E}^2(i_g, i_h)\), įvestą Elemento mazgų fizikinių dydžių surinkimo skyriuje:

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

Gautos reikšmės išdėstomos matricoje, kurios eilutė yra \(i_g\), o stulpelis – \(i_h\). Realizacijoje aibė \(\mathcal{E}^2\) aiškiai nesudaroma; vietoj to atitinkami blokai tiesiogiai pridedami elementų cikle. Matrica yra kvadratinė, jos matmuo yra laisvės laipsnių skaičius viename mazge \(\times\) bendras mazgų skaičius \(n_g\), tačiau kadangi komponentai tarp mazgų, nesusietų elementais, yra \(0\), ji saugoma retosios matricos forma.

TL ir UL formuluočių elementų standumo matricos yra tos pačios formos, išskyrus atskaitos konfigūracijos pakeitimą (mazgų koordinates ir B matricos sudarymo šaltinį) bei \(\boldsymbol{G}\) matricos buvimą arba nebuvimą. Todėl FrontISTR abi formuluotes realizuoja bendroje paprogramėje.

Iteracijų algoritmas

Apibendrinant, iteracijų pradžioje nustatoma \(\Delta\boldsymbol{u} = \boldsymbol{0}\) ir apskaičiuojama pradinė liekamoji jėga \(\boldsymbol{R}_0 = \boldsymbol{F}_{n+1} - \boldsymbol{Q}(\boldsymbol{u}_n)\). Tada \(i\)-ojoje iteracijoje atliekama ši procedūra.

  1. Esant dabartiniam poslinkiui \(\boldsymbol{u}_n + \Delta\boldsymbol{u}\), pagal Tangentinio standumo matricos sudarymo procedūrą apskaičiuojamas tangentinis standumas \(\boldsymbol{K}_i = \partial\boldsymbol{Q}/\partial\boldsymbol{u}|_{\boldsymbol{u}_n + \Delta\boldsymbol{u}}\).
  2. Geometrinėms kraštinėms sąlygoms taikyti tangentinio standumo matrica ir liekamosios jėgos vektorius modifikuojami poslinkio apribojimus turintiems laisvės laipsniams; gaunama \(\tilde{\boldsymbol{K}}_i, \tilde{\boldsymbol{R}}_{i-1}\) (Geometrinių kraštinių sąlygų apdorojimas).
  3. Išsprendžiama tiesinė lygtis \(\tilde{\boldsymbol{K}}_i\, d\boldsymbol{u}_i = \tilde{\boldsymbol{R}}_{i-1}\) ir gaunama pataisa \(d\boldsymbol{u}_i\). Ši procedūra dažnai sudaro didžiausią iteracinio skaičiavimo sąnaudų dalį.
  4. Poslinkio prieaugis atnaujinamas pagal \(\Delta\boldsymbol{u} \leftarrow \Delta\boldsymbol{u} + d\boldsymbol{u}_i\); pagal jį apskaičiuojamas vidinių jėgų vektorius \(\boldsymbol{Q}(\boldsymbol{u}_n + \Delta\boldsymbol{u})\) ir liekamoji jėga \(\boldsymbol{R}_i = \boldsymbol{F}_{n+1} - \boldsymbol{Q}(\boldsymbol{u}_n + \Delta\boldsymbol{u})\).
  5. Patikrinama konvergencija ir, jei ji pasiekta, iteracijos baigiamos. Geometrines kraštines sąlygas turinčiuose laisvės laipsniuose liekamojoje jėgoje \(\boldsymbol{R}_i\) atsiranda suvaržytas reakcijas atitinkantys komponentai, todėl konvergencijos rodiklis sudaromas iš \(\tilde{\boldsymbol{R}}_i\), pašalinus šiuos komponentus. Konkretūs konvergencijos rodikliai ir slenksčiai aprašyti Konvergencijos kriterijų skyriuje. Jei konvergencija nepasiekiama ir pasiekiama iteracijų riba, iteracija laikoma nesėkminga.

Iteracijai konvergavus, konverguotas \(\Delta\boldsymbol{u}\) pridedamas prie \(\boldsymbol{u}_n\), gaunant sukauptą poslinkį laiku \(t_{n+1}\), \(\boldsymbol{u}_{n+1} = \boldsymbol{u}_n + \Delta\boldsymbol{u}\), ir pereinama prie kito laiko žingsnio.

Susijusios temos

AI-assisted translation May contain errors Official docs Status