Newton–Raphsonin menetelmä

Linearisointi ja iteratiivinen toistokaava

Ulkoisten voimien virtuaalityö ja globaalin yhtälön kokoaminen antaa ajan \(t_{n+1}\) solmusiirtymälle \(\boldsymbol{u}_{n+1}\) epälineaarisen yhtälön, joka ratkaistaan Newton–Raphsonin menetelmällä. Oletetaan, että solmusiirtymä ajanhetkeen \(t_n\) saakka, \(\boldsymbol{u}_n\), tunnetaan, ja ratkaistavaksi tuntemattomaksi muuttujaksi otetaan siirtymäinkrementti \(\Delta\boldsymbol{u}\)

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

Seuraavassa ulkoisen voimavektorin riippuvuus solmusiirtymästä jätetään huomiotta, ja kun \(\boldsymbol{F}(\boldsymbol{u}_{n+1}) = \boldsymbol{F}_{n+1}\),

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

ratkaistaan.

Määritellään nykyisessä ratkaisussa \(\Delta\boldsymbol{u}\) tangenttijäykkyys

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

Tämän avulla epälineaarisen yhtälön linearisointi antaa

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

Olkoon \(i\):nnen iteraation korjaus \(d\boldsymbol{u}_i\) ja iteraation alun residuaalivektori

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

Tällöin iteratiivinen toistokaava on

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

Residuaali \(\boldsymbol{R}_i\) on siis tasapainotilasta poikkeavaa voimien epätasapainoa vastaava suure.

Tangenttijäykkyysmatriisin muodostaminen

Tangenttijäykkyys \(\boldsymbol{K} = \partial\boldsymbol{Q}/\partial\boldsymbol{u}\) muodostetaan derivoimalla sisäisten voimien virtuaalityön diskretoinnissa saatu elementin sisäisen voiman vektori osittain solmusiirtymän suhteen, integroimalla näin saadut elementtitason integroitavat funktiot kunkin elementin alueella ja kokoamalla ne. Kun elementtitason integroitavaa merkitään \(\boldsymbol{K}^e_X\):llä (referenssikonfiguraation merkintä, TL-formulointi) tai \(\boldsymbol{K}^e_x\):llä (nykykonfiguraation merkintä, UL-formulointi), elementin tangenttijäykkyys on

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

Seuraavassa esitetään TL- ja UL-integroitavien lopulliset muodot. Molemmissa tapauksissa ne jaetaan materiaalijäykkyystermin (alkusiirtymätermin) ja geometrisen jäykkyystermin (alkujännitystermin) summaan.

Total Lagrangen formulointi

Total Lagrangen formuloinnissa oletetaan lineaarinen yhteys toisen Piola–Kirchhoffin jännityksen nopeuden \(\dot{\boldsymbol{S}}\) ja Green–Lagrangen venymänopeuden \(\dot{\boldsymbol{E}}\) välillä, eli \(\dot{\boldsymbol{S}} = \boldsymbol{\mathsf{C}}:\dot{\boldsymbol{E}}\). Tämä vastaa lineaarisesti elastisten materiaalien (St. Venant–Kirchhoffin materiaalien) ja hyperelastisten materiaalien konstitutiivisia lakeja, ja FrontISTR käyttää näille materiaaleille Total Lagrangen formulointia. Elementin tangenttijäykkyyden integroitava kirjoitetaan tällöin tensorimuodossa

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

Oikean puolen ensimmäinen termi on materiaalijäykkyystermi (alkusiirtymätermi) ja toinen termi geometrinen jäykkyystermi (alkujännitystermi).

FrontISTRin toteutuksessa tämä integroitava arvioidaan matriisimuodossa Voigtin merkintää käyttäen:

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

Matriisit ovat seuraavat. \(\boldsymbol{B}_L, \boldsymbol{B}_{NL}\) ovat sisäisten voimien virtuaalityön diskretoinnissa esitellyt B-matriisit, ja \(\tilde{\boldsymbol{C}}\) on konstitutiivisen tensorin \(\boldsymbol{\mathsf{C}}\) Voigt-esitys eli \(6\times 6\) materiaalijäykkyysmatriisi (Tensorimerkintä ja matemaattiset perusteet). \(\boldsymbol{S}_9, \boldsymbol{F}_9\) ovat seuraavat uudelleenjärjestelymatriisit, joilla geometrinen jäykkyystermi esitetään matriisitulona. Määritellään ensin \(3\times 3\) toisen kertaluvun tensorille \(\boldsymbol{A}\) merkintä \([\,\cdot\,]\), joka järjestää sen 9-komponenttiseksi vektoriksi

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

Tällä määrittelyllä \(\boldsymbol{F}_9\) esittää muodonmuutosgradientin variaation muodossa \([\delta\boldsymbol{F}] = \boldsymbol{F}_9\, \delta\boldsymbol{u}^e\) ja on \(9\times d n_e\) -matriisi. Elementin solmulle \(\alpha = 1, \ldots, n_e\) vastaava \(9\times d\) -lohko on

\[ [\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{ on } 3\times 3 \text{ identiteettimatriisi}) \]

ja saadaan muodossa \(\boldsymbol{F}_9 = [[\boldsymbol{F}_9]_1, \ldots, [\boldsymbol{F}_9]_{n_e}]\), jossa lohkot on järjestetty vaakasuunnassa elementin solmujärjestyksessä. \(\boldsymbol{S}_9\) valitaan siten, että yhdessä tämän matriisin kanssa geometrinen jäykkyystermi voidaan esittää muodossa \(\delta\boldsymbol{u}^{eT}\, \boldsymbol{F}_9^T \boldsymbol{S}_9 \boldsymbol{F}_9\, \dot{\boldsymbol{u}}^e\); se on seuraava \(9\times 9\) -matriisi

\[ \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ästä saadaan kyseinen matriisi.

Updated Lagrangen formulointi

Updated Lagrangen formuloinnissa oletetaan lineaarinen yhteys suhteellisen Kirchhoffin jännitystensorin Jaumann-nopeuden \(\hat{\boldsymbol{\sigma}}^{\nabla J}\) ja muodonmuutosnopeustensorin \(\boldsymbol{D}\) välillä, eli \(\hat{\boldsymbol{\sigma}}^{\nabla J} = \boldsymbol{\mathsf{C}}:\boldsymbol{D}\). Tämä on lineaarisesti elastisille, elastoplastisille ja virumismateriaaleille yhteinen hypoelastisen konstitutiivisen lain muoto, ja FrontISTR käyttää näille materiaaleille Updated Lagrangen formulointia. Nykykonfiguraatiossa ilmaistu elementin tangenttijäykkyyden integroitava kirjoitetaan tällöin tensorimuodossa

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

missä \(\boldsymbol{\sigma}^{\nabla T}\) on Truesdellin nopeus, \(\boldsymbol{A}_{(L)}\) Almansi-venymän lineaarinen osa, \(\boldsymbol{F}_t = \partial\boldsymbol{u}/\partial\boldsymbol{x}\) nykykonfiguraation suhteen määritetty siirtymägradientti ja \(\boldsymbol{L}\) nopeusgradienttitensori. Oikean puolen ensimmäinen termi on materiaalijäykkyystermi ja toinen termi geometrinen jäykkyystermi.

FrontISTRin toteutuksessa tämä integroitava arvioidaan matriisimuodossa Voigtin merkintää käyttäen:

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

Tässä \(\boldsymbol{b}\) on nykykonfiguraatiossa muodostettu B-matriisi (sisäisten voimien virtuaalityön diskretointi). \(\boldsymbol{\sigma}_9, \boldsymbol{f}_9\) saadaan TL-formulointia varten määritellyistä \(\boldsymbol{S}_9, \boldsymbol{F}_9\) korvaamalla toinen PK-jännitys \(\boldsymbol{S}\) Cauchyn jännityksellä \(\boldsymbol{\sigma}\) ja referenssikonfiguraation gradientti \(\partial N_\alpha^e/\partial X_i\) nykykonfiguraation gradientilla \(\partial N_\alpha^e/\partial x_i\).

\(\boldsymbol{G}\) on Cauchyn jännityksestä riippuva korjausmatriisi, jota tarvitaan saattamaan hypoelastinen konstitutiivinen laki \(\hat{\boldsymbol{\sigma}}^{\nabla J} = \boldsymbol{\mathsf{C}}:\boldsymbol{D}\) yhteensopivaksi tangenttijäykkyyskehyksen kanssa Truesdell-nopeuteen perustuvana konstitutiivisena lakina. Se saadaan järjestämällä neljännen kertaluvun tensorikomponentit \(G_{ijkl} = \delta_{il}\sigma_{kj} + \delta_{jl}\sigma_{ik}\) \(6\times 6\) Voigt-muotoon seuraavasti

\[ \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ästä saadaan kyseinen matriisi.

Globaalin jäykkyysmatriisin kokoaminen

Globaali tangenttijäykkyys \(\boldsymbol{K}\) saadaan jakamalla kunkin elementin jäykkyys \(\boldsymbol{K}^e\) \(d\times d\) -lohkoihin \(\boldsymbol{K}^e_{\alpha\beta}\) kullekin solmuparille ja käyttämällä elementin solmufysikaalisten suureiden kokoamisessa esiteltyä toisen kertaluvun tensorien kokoamisjoukkoa \(\mathcal{E}^2(i_g, i_h)\):

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

Saadut arvot järjestetään matriisiksi, jonka rivi on \(i_g\) ja sarake \(i_h\). Toteutuksessa joukkoa \(\mathcal{E}^2\) ei muodosteta eksplisiittisesti, vaan vastaavat lohkot lisätään suoraan elementtisilmukan sisällä. Matriisi on neliömatriisi, jonka dimensio on solmukohtaisten vapausasteiden määrä \(\times\) solmujen kokonaismäärä \(n_g\), mutta koska muut kuin elementtien kautta toisiinsa kytkettyjen solmujen väliset komponentit ovat \(0\), matriisi tallennetaan harvana matriisina.

TL- ja UL-formulointien elementtijäykkyysmatriiseilla on sama muoto lukuun ottamatta referenssikonfiguraation vaihtumista (solmukoordinaatit ja B-matriisin muodostamiseen käytetty lähde) sekä \(\boldsymbol{G}\)-matriisin olemassaoloa tai puuttumista. FrontISTR toteuttaa siksi molemmat formuloinnit yhteisessä aliohjelmassa.

Iteraatioalgoritmi

Yhteenvetona asetetaan iteraation alussa \(\Delta\boldsymbol{u} = \boldsymbol{0}\) ja lasketaan alkuperäinen residuaali \(\boldsymbol{R}_0 = \boldsymbol{F}_{n+1} - \boldsymbol{Q}(\boldsymbol{u}_n)\). Tämän jälkeen \(i\):nnellä iteraatiolla suoritetaan seuraava menettely.

1. Laske nykyisellä siirtymällä \(\boldsymbol{u}_n + \Delta\boldsymbol{u}\) tangenttijäykkyys \(\boldsymbol{K}_i = \partial\boldsymbol{Q}/\partial\boldsymbol{u}|_{\boldsymbol{u}_n + \Delta\boldsymbol{u}}\) kohdassa Tangenttijäykkyysmatriisin muodostaminen kuvatulla menettelyllä.
2. Geometristen reunaehtojen asettamiseksi muokkaa tangenttijäykkyysmatriisia ja residuaalivektoria siirtymärajoitusten alaisille vapausasteille ja saa näin \(\tilde{\boldsymbol{K}}_i, \tilde{\boldsymbol{R}}_{i-1}\) (Geometristen reunaehtojen käsittely).
3. Ratkaise lineaarinen yhtälö \(\tilde{\boldsymbol{K}}_i\, d\boldsymbol{u}_i = \tilde{\boldsymbol{R}}_{i-1}\) ja määritä korjaus \(d\boldsymbol{u}_i\). Tämä vaihe muodostaa usein suurimman osan iteratiivisen laskennan laskentakustannuksesta.
4. Päivitä siirtymäinkrementti muodossa \(\Delta\boldsymbol{u} \leftarrow \Delta\boldsymbol{u} + d\boldsymbol{u}_i\) ja laske vastaavasti sisäisen voiman vektori \(\boldsymbol{Q}(\boldsymbol{u}_n + \Delta\boldsymbol{u})\) sekä residuaali \(\boldsymbol{R}_i = \boldsymbol{F}_{n+1} - \boldsymbol{Q}(\boldsymbol{u}_n + \Delta\boldsymbol{u})\).
5. Tarkista konvergenssi ja lopeta iterointi, jos konvergenssi on saavutettu. Geometristen reunaehtojen alaisiin vapausasteisiin liittyvät rajoitusreaktiot näkyvät residuaalissa \(\boldsymbol{R}_i\), joten konvergenssi-indikaattori muodostetaan \(\tilde{\boldsymbol{R}}_i\):stä sen jälkeen, kun nämä komponentit on poistettu. Erityiset konvergenssi-indikaattorit ja raja-arvot kuvataan kohdassa Konvergenssikriteerit. Jos konvergenssia ei saavuteta ja iteraatioraja täyttyy, iterointi katsotaan epäonnistuneeksi.

Kun iterointi konvergoi, lisää konvergoitunut \(\Delta\boldsymbol{u}\) arvoon \(\boldsymbol{u}_n\), jolloin saadaan ajan \(t_{n+1}\) kumulatiivinen siirtymä \(\boldsymbol{u}_{n+1} = \boldsymbol{u}_n + \Delta\boldsymbol{u}\), ja siirry seuraavaan aika-askeleeseen.

Liittyvät aiheet

• Ulkoisten voimien virtuaalityö ja globaalin yhtälön kokoaminen — Ratkaistavan epälineaarisen yhtälön lähtökohta
• Sisäisten voimien virtuaalityön diskretointi — Matriisien \(\boldsymbol{B}_L, \boldsymbol{B}_{NL}, \boldsymbol{b}\) muodostaminen
• Geometristen reunaehtojen käsittely — Tangenttijäykkyysmatriisin ja residuaalivektorin muokkaaminen siirtymärajoitusten asettamiseksi
• Konvergenssikriteerit — Residuaalin normiin perustuvat pysäytyskriteerit
• Tensorimerkintä ja matemaattiset perusteet — Materiaalimatriisin \(\tilde{\boldsymbol{C}}\) Voigt-esitys
• Epälineaarinen iterointi ja aikaintegrointi (toiminnot) — Käyttö ja valinta toimintoviitteessä