Keri sisuni

Newtoni-Raphsoni meetod

Lineariseerimine ja iteratiivne rekurrents

Jaotises Välisjõudude virtuaaltöö ja globaalse võrrandi koostamine saadakse ajahetke \(t_{n+1}\) sõlmenihke \(\boldsymbol{u}_{n+1}\) jaoks mittelineaarne võrrand, mis lahendatakse Newtoni-Raphsoni meetodiga. Eeldatakse, et sõlmeniht kuni ajahetkeni \(t_n\), \(\boldsymbol{u}_n\), on teada, ning määratavaks tundmatuks võetakse nihke juurdekasv \(\Delta\boldsymbol{u}\)

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

Edaspidi jäetakse välisjõuvektori sõltuvus sõlmenihkest arvestamata ning tingimusega \(\boldsymbol{F}(\boldsymbol{u}_{n+1}) = \boldsymbol{F}_{n+1}\) lahendatakse

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

.

Praeguse lahendi \(\Delta\boldsymbol{u}\) juures defineeritakse puutujäikus

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

Seda kasutades annab mittelineaarse võrrandi lineariseerimine

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

Olgu \(i\)-nda iteratsiooni parand \(d\boldsymbol{u}_i\) ning iteratsiooni alguse jääkvektor

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

Siis on iteratiivne rekurrents

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

Seega on jääk \(\boldsymbol{R}_i\) suurus, mis väljendab tasakaalust kõrvalekalduvat jõudude tasakaalutust.

Puutujäikusmaatriksi koostamine

Puutujäikus \(\boldsymbol{K} = \partial\boldsymbol{Q}/\partial\boldsymbol{u}\) koostatakse, diferentseerides jaotises Sisejõudude virtuaaltöö diskretiseerimine saadud elemendi sisejõuvektori osaliselt sõlmenihke suhtes, integreerides saadud elemenditaseme integraalavaldised iga elemendi piirkonnas ja need kokku monteerides. Tähistades elemenditaseme integraalavaldisi \(\boldsymbol{K}^e_X\)-ga (lähtekonfiguratsiooni tähistus, TL-formulatsioon) või \(\boldsymbol{K}^e_x\)-ga (hetkekonfiguratsiooni tähistus, UL-formulatsioon), on elemendi puutujäikus

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

Allpool esitatakse TL- ja UL-integraalavaldiste lõppkujud. Mõlemal juhul lagunevad need materjalijäikuse liikme (algnihke liikme) ja geomeetrilise jäikuse liikme (algpinge liikme) summaks.

Total Lagrange'i formulatsioon

Total Lagrange'i formulatsioonis eeldatakse lineaarset seost 2. Piola-Kirchhoffi pingetensori muutumiskiiruse \(\dot{\boldsymbol{S}}\) ja Green-Lagrange'i deformatsioonitensori muutumiskiiruse \(\dot{\boldsymbol{E}}\) vahel, s.t. \(\dot{\boldsymbol{S}} = \boldsymbol{\mathsf{C}}:\dot{\boldsymbol{E}}\). See vastab lineaarsete elastsete materjalide (St. Venant-Kirchhoffi materjalide) ja hüperelastsete materjalide konstitutiivvõrranditele ning FrontISTR kasutab nende materjalide puhul Total Lagrange'i formulatsiooni. Elemendi puutujäikuse integraalavaldis kirjutatakse siis tensorkujul

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

Parema poole 1. liige on materjalijäikuse liige (algnihke liige) ja 2. liige geomeetrilise jäikuse liige (algpinge liige).

FrontISTR-i realisatsioonis hinnatakse seda integraalavaldist maatrikskujul, kasutades Voigti tähistust:

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

Maatriksid on järgmised. \(\boldsymbol{B}_L, \boldsymbol{B}_{NL}\) on jaotises Sisejõudude virtuaaltöö diskretiseerimine kasutusele võetud B-maatriksid ning \(\tilde{\boldsymbol{C}}\) on konstitutiivtensori \(\boldsymbol{\mathsf{C}}\) Voigti esitus ehk \(6\times 6\) materjalijäikusmaatriks (Tensoritähistus ja matemaatilised alused). \(\boldsymbol{S}_9, \boldsymbol{F}_9\) on järgmised ümberjärjestusmaatriksid, mida kasutatakse geomeetrilise jäikuse liikme esitamiseks maatrikskorrutisena. Kõigepealt defineeritakse \(3\times 3\) 2. järku tensori \(\boldsymbol{A}\) jaoks tähistus \([\,\cdot\,]\), mis järjestab selle 9-komponendiliseks vektoriks:

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

Selle definitsiooni korral väljendab \(\boldsymbol{F}_9\) deformatsioonigradiendi variatsiooni kujul \([\delta\boldsymbol{F}] = \boldsymbol{F}_9\, \delta\boldsymbol{u}^e\) ning on \(9\times d n_e\) maatriks. Elemendi sõlme \(\alpha = 1, \ldots, n_e\) jaoks on vastav \(9\times d\) plokk

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

ning \(\boldsymbol{F}_9 = [[\boldsymbol{F}_9]_1, \ldots, [\boldsymbol{F}_9]_{n_e}]\), kus plokid on paigutatud horisontaalselt elemendi sõlmede järjekorras. \(\boldsymbol{S}_9\) valitakse nii, et koos selle maatriksiga avalduks geomeetrilise jäikuse liige kujul \(\delta\boldsymbol{u}^{eT}\, \boldsymbol{F}_9^T \boldsymbol{S}_9 \boldsymbol{F}_9\, \dot{\boldsymbol{u}}^e\); see on järgmine \(9\times 9\) maatriks

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

See on saadud maatriks.

Updated Lagrange'i formulatsioon

Updated Lagrange'i formulatsioonis eeldatakse lineaarset seost suhtelise Kirchhoffi pingetensori Jaumanni kiiruse \(\hat{\boldsymbol{\sigma}}^{\nabla J}\) ja deformatsioonikiiruse tensori \(\boldsymbol{D}\) vahel, s.t. \(\hat{\boldsymbol{\sigma}}^{\nabla J} = \boldsymbol{\mathsf{C}}:\boldsymbol{D}\). See on hüpoelastse konstitutiivvõrrandi kuju, mis on ühine lineaarsetele elastsetele, elastoplastsetele ja roomematerjalidele, ning FrontISTR kasutab nende materjalide jaoks Updated Lagrange'i formulatsiooni. Hetkekonfiguratsioonis väljendatud elemendi puutujäikuse integraalavaldis kirjutatakse siis tensorkujul

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

kus \(\boldsymbol{\sigma}^{\nabla T}\) on Truesdelli kiirus, \(\boldsymbol{A}_{(L)}\) Almansi deformatsiooni lineaarne osa, \(\boldsymbol{F}_t = \partial\boldsymbol{u}/\partial\boldsymbol{x}\) hetkekonfiguratsiooni suhtes võetud nihkegradient ja \(\boldsymbol{L}\) kiirusgradiendi tensor. Parema poole 1. liige on materjalijäikuse liige ja 2. liige geomeetrilise jäikuse liige.

FrontISTR-i realisatsioonis hinnatakse seda integraalavaldist maatrikskujul, kasutades Voigti tähistust:

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

Siin on \(\boldsymbol{b}\) hetkekonfiguratsioonis koostatud B-maatriks (Sisejõudude virtuaaltöö diskretiseerimine). \(\boldsymbol{\sigma}_9, \boldsymbol{f}_9\) saadakse TL-formulatsiooni jaoks määratletud \(\boldsymbol{S}_9, \boldsymbol{F}_9\) maatriksitest, asendades 2. PK-pinge \(\boldsymbol{S}\) Cauchy pingega \(\boldsymbol{\sigma}\) ja lähtekonfiguratsiooni gradiendi \(\partial N_\alpha^e/\partial X_i\) hetkekonfiguratsiooni gradiendiga \(\partial N_\alpha^e/\partial x_i\).

\(\boldsymbol{G}\) on Cauchy pingest sõltuv parandusmaatriks, mis on vajalik, et hüpoelastne konstitutiivvõrrand \(\hat{\boldsymbol{\sigma}}^{\nabla J} = \boldsymbol{\mathsf{C}}:\boldsymbol{D}\) oleks kooskõlas puutujäikuse raamistikuga kui Truesdelli kiirusel põhineva konstitutiivvõrrandiga. See saadakse 4. järku tensori komponentide \(G_{ijkl} = \delta_{il}\sigma_{kj} + \delta_{jl}\sigma_{ik}\) paigutamisel \(6\times 6\) Voigti kujule:

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

See on saadud maatriks.

Globaalse jäikusmaatriksi koostamine

Globaalne puutujäikus \(\boldsymbol{K}\) saadakse, jagades iga elemendi jäikuse \(\boldsymbol{K}^e\) iga sõlmepaari jaoks \(d\times d\) plokkideks \(\boldsymbol{K}^e_{\alpha\beta}\) ja kasutades jaotises Elemendi-sõlme füüsikaliste suuruste koostamine kasutusele võetud 2. järku tensori koostehulka \(\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} \]

Saadud väärtused paigutatakse maatriksisse rea \(i_g\) ja veeru \(i_h\) järgi. Realisatsioonis hulka \(\mathcal{E}^2\) otseselt ei moodustata; selle asemel liidetakse vastavad plokid elemenditsüklis otse. Maatriks on ruutmaatriks mõõtmega vabadusastmete arv sõlme kohta \(\times\) sõlmede koguarv \(n_g\), kuid kuna elementidega ühendamata sõlmede vahelised komponendid on \(0\), salvestatakse see hõremaatriksi kujul.

TL- ja UL-formulatsiooni elemendi jäikusmaatriksid on sama kujuga, erinevuseks on üksnes lähtekonfiguratsiooni vahetamine (sõlmekoordinaadid ja B-maatriksi koostamiseks kasutatav lähteandmestik) ning maatriksi \(\boldsymbol{G}\) olemasolu või puudumine. Seetõttu realiseerib FrontISTR mõlemad formulatsioonid ühises alamprogrammis.

Iteratsioonialgoritm

Eeltoodu kokkuvõttena seatakse iteratsiooni alguses \(\Delta\boldsymbol{u} = \boldsymbol{0}\) ja arvutatakse algjääk \(\boldsymbol{R}_0 = \boldsymbol{F}_{n+1} - \boldsymbol{Q}(\boldsymbol{u}_n)\). Seejärel tehakse \(i\)-ndal iteratsioonil järgmised sammud.

  1. Praeguse nihke \(\boldsymbol{u}_n + \Delta\boldsymbol{u}\) juures arvutatakse puutujäikus \(\boldsymbol{K}_i = \partial\boldsymbol{Q}/\partial\boldsymbol{u}|_{\boldsymbol{u}_n + \Delta\boldsymbol{u}}\) jaotises Puutujäikusmaatriksi koostamine kirjeldatud protseduuriga.
  2. Geomeetriliste rajatingimuste rakendamiseks muudetakse nihkepiirangutega vabadusastmete puutujäikusmaatriksit ja jääkvektorit ning saadakse \(\tilde{\boldsymbol{K}}_i, \tilde{\boldsymbol{R}}_{i-1}\) (Geomeetriliste rajatingimuste käsitlus).
  3. Lahendatakse lineaarvõrrand \(\tilde{\boldsymbol{K}}_i\, d\boldsymbol{u}_i = \tilde{\boldsymbol{R}}_{i-1}\) ja leitakse parand \(d\boldsymbol{u}_i\). See samm moodustab sageli suurema osa iteratiivarvutuse arvutuskulust.
  4. Nihke juurdekasvu uuendatakse kujul \(\Delta\boldsymbol{u} \leftarrow \Delta\boldsymbol{u} + d\boldsymbol{u}_i\) ning vastavalt arvutatakse sisejõuvektor \(\boldsymbol{Q}(\boldsymbol{u}_n + \Delta\boldsymbol{u})\) ja jääk \(\boldsymbol{R}_i = \boldsymbol{F}_{n+1} - \boldsymbol{Q}(\boldsymbol{u}_n + \Delta\boldsymbol{u})\).
  5. Kontrollitakse koonduvust ja koonduvuse saavutamisel iteratsioon lõpetatakse. Geomeetriliste rajatingimustega seotud vabadusastmete jäägis \(\boldsymbol{R}_i\) esinevad piiratud reaktsioonidele vastavad komponendid, mistõttu koonduvusnäitaja moodustatakse pärast nende komponentide eemaldamist vektorist \(\tilde{\boldsymbol{R}}_i\). Konkreetsed koonduvusnäitajad ja läved on kirjeldatud jaotises Koonduvuskriteeriumid. Kui koonduvust ei saavutata ja iteratsioonipiir saab täis, loetakse iteratsioon ebaõnnestunuks.

Kui iteratsioon koondub, liidetakse koondunud \(\Delta\boldsymbol{u}\) väärtusele \(\boldsymbol{u}_n\), et saada ajahetke \(t_{n+1}\) kumulatiivne nihe \(\boldsymbol{u}_{n+1} = \boldsymbol{u}_n + \Delta\boldsymbol{u}\), ning liigutakse järgmise ajasammu juurde.

Seotud teemad