Hyppää sisältöön

Newton–Raphsonin menetelmä

Linearisointi ja iteratiivinen toistokaava

Ulkoisten voimien virtuaalinen työ ja globaalien yhtälöiden 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äisen 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 2. 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 1. termi on materiaalijäykkyystermi (alkusiirtymätermi) ja 2. 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äisen 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\) 2. 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 1. termi on materiaalijäykkyystermi ja 2. 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äisen virtuaalityön diskretointi). \(\boldsymbol{\sigma}_9, \boldsymbol{f}_9\) saadaan TL-formulointia varten määritellyistä \(\boldsymbol{S}_9, \boldsymbol{F}_9\) korvaamalla 2. 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ä 4. 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 solmusuureiden kokoamisessa esiteltyä 2. 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