Ņūtona–Rafsona metode¶

Linearizācija un iteratīvā rekurence¶

Sadaļā Ārējo spēku virtuālais darbs un globālā vienādojuma montāža iegūts nelineārs vienādojums mezglu pārvietojumam laikā \(t_{n+1}\), \(\boldsymbol{u}_{n+1}\), ko risina ar Ņūtona–Rafsona metodi. Tiek pieņemts, ka mezglu pārvietojums līdz laikam \(t_n\), \(\boldsymbol{u}_n\), ir zināms, un par nosakāmo nezināmo mainīgo izvēlas pārvietojuma pieaugumu \(\Delta\boldsymbol{u}\)

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

Turpmāk ārējo spēku vektora atkarība no mezglu pārvietojuma tiek ignorēta, un, pieņemot \(\boldsymbol{F}(\boldsymbol{u}_{n+1}) = \boldsymbol{F}_{n+1}\),

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

tiek atrisināts.

Pie pašreizējā risinājuma \(\Delta\boldsymbol{u}\) definē pieskares stingumu

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

Izmantojot to, nelineārā vienādojuma linearizācija dod

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

Lai korekcija \(i\)-tajā iterācijā būtu \(d\boldsymbol{u}_i\), bet atlikuma vektors iterācijas sākumā būtu

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

Tad iteratīvā rekurence ir

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

Tādējādi atlikums \(\boldsymbol{R}_i\) ir lielums, kas atbilst spēku nelīdzsvarotībai attiecībā pret līdzsvara stāvokli.

Pieskares stinguma matricas izveide¶

Pieskares stingumu \(\boldsymbol{K} = \partial\boldsymbol{Q}/\partial\boldsymbol{u}\) izveido, daļēji diferencējot sadaļā Iekšējo spēku virtuālā darba diskretizācija iegūto elementa iekšējo spēku vektoru pēc mezglu pārvietojuma, integrējot iegūtos elementa līmeņa integrandus katra elementa apgabalā un tos samontējot. Apzīmējot elementa līmeņa integrandu ar \(\boldsymbol{K}^e_X\) (atskaites konfigurācijas pieraksts, TL formulējums) vai \(\boldsymbol{K}^e_x\) (pašreizējās konfigurācijas pieraksts, UL formulējums), elementa pieskares stingums ir

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

Tālāk dotas TL/UL integrandu galīgās formas. Abos gadījumos tie tiek sadalīti materiāla stinguma locekļa (sākotnējā pārvietojuma locekļa) un ģeometriskā stinguma locekļa (sākotnējā sprieguma locekļa) summā.

Pilnais Lagranža formulējums¶

Pilnajā Lagranža formulējumā pieņem lineāru sakarību starp otrā Piolas–Kirhofa sprieguma ātrumu \(\dot{\boldsymbol{S}}\) un Grīna–Lagranža deformācijas ātrumu \(\dot{\boldsymbol{E}}\), proti, \(\dot{\boldsymbol{S}} = \boldsymbol{\mathsf{C}}:\dot{\boldsymbol{E}}\). Tas atbilst lineāri elastīgu materiālu (Senvenāna–Kirhofa materiālu) un hiperelastīgu materiālu konstitutīvajiem likumiem, un FrontISTR šiem materiāliem izmanto pilno Lagranža formulējumu. Tad elementa pieskares stinguma integrandu tenzoru formā raksta kā

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

Pirmais loceklis labajā pusē ir materiāla stinguma loceklis (sākotnējā pārvietojuma loceklis), bet otrais — ģeometriskā stinguma loceklis (sākotnējā sprieguma loceklis).

FrontISTR realizācijā šo integrandu novērtē matricas formā, izmantojot Foigta pierakstu:

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

Katra matrica ir šāda. \(\boldsymbol{B}_L, \boldsymbol{B}_{NL}\) ir B matricas, kas ieviestas sadaļā Iekšējo spēku virtuālā darba diskretizācija, bet \(\tilde{\boldsymbol{C}}\) ir konstitutīvā tenzora \(\boldsymbol{\mathsf{C}}\) Foigta attēlojums, proti, \(6\times 6\) materiāla stinguma matrica (Tenzoru pieraksts un matemātiskie pamati). \(\boldsymbol{S}_9, \boldsymbol{F}_9\) ir tālāk dotās pārkārtošanas matricas, ko izmanto, lai ģeometriskā stinguma locekli izteiktu kā matricu reizinājumu. Vispirms \(3\times 3\) otrās kārtas tenzoram \(\boldsymbol{A}\) definē pierakstu \([\,\cdot\,]\), kas to pārkārto 9 komponentu vektorā:

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

Ar šo definīciju \(\boldsymbol{F}_9\) izsaka deformācijas gradienta variāciju formā \([\delta\boldsymbol{F}] = \boldsymbol{F}_9\, \delta\boldsymbol{u}^e\) un ir \(9\times d n_e\) matrica. Elementa mezglam \(\alpha = 1, \ldots, n_e\) atbilstošais \(9\times d\) bloks ir

\[ [\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{ ir } 3\times 3 \text{ vienības matrica}) \]

un \(\boldsymbol{F}_9 = [[\boldsymbol{F}_9]_1, \ldots, [\boldsymbol{F}_9]_{n_e}]\), blokus horizontāli sakārtojot elementa mezglu secībā. \(\boldsymbol{S}_9\) izvēlas tā, lai kopā ar šo matricu ģeometriskā stinguma locekli izteiktu kā \(\delta\boldsymbol{u}^{eT}\, \boldsymbol{F}_9^T \boldsymbol{S}_9 \boldsymbol{F}_9\, \dot{\boldsymbol{u}}^e\); tā ir šāda \(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} \]

Šī ir iegūtā matrica.

Atjauninātais Lagranža formulējums¶

Atjauninātajā Lagranža formulējumā pieņem lineāru sakarību starp relatīvā Kirhofa sprieguma tenzora Jaumana atvasinājumu \(\hat{\boldsymbol{\sigma}}^{\nabla J}\) un deformācijas ātruma tenzoru \(\boldsymbol{D}\), proti, \(\hat{\boldsymbol{\sigma}}^{\nabla J} = \boldsymbol{\mathsf{C}}:\boldsymbol{D}\). Tā ir hipoelastīga konstitutīvā likuma forma, kas kopīga lineāri elastīgiem, elastoplastiskiem un šļūdes materiāliem, un FrontISTR šiem materiāliem izmanto atjaunināto Lagranža formulējumu. Tad elementa pieskares stinguma integrandu, izteiktu pašreizējā konfigurācijā, tenzoru formā raksta kā

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

kur \(\boldsymbol{\sigma}^{\nabla T}\) ir Trūsdela ātrums, \(\boldsymbol{A}_{(L)}\) ir Almansi deformācijas lineārā daļa, \(\boldsymbol{F}_t = \partial\boldsymbol{u}/\partial\boldsymbol{x}\) ir pārvietojuma gradients attiecībā pret pašreizējo konfigurāciju, bet \(\boldsymbol{L}\) ir ātruma gradienta tenzors. Pirmais loceklis labajā pusē ir materiāla stinguma loceklis, bet otrais — ģeometriskā stinguma loceklis.

FrontISTR realizācijā šo integrandu novērtē matricas formā, izmantojot Foigta pierakstu:

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

Šeit \(\boldsymbol{b}\) ir B matrica, kas izveidota pašreizējā konfigurācijā (Iekšējo spēku virtuālā darba diskretizācija). \(\boldsymbol{\sigma}_9, \boldsymbol{f}_9\) iegūst no TL formulējumam definētajām \(\boldsymbol{S}_9, \boldsymbol{F}_9\), aizstājot otro PK spriegumu \(\boldsymbol{S}\) ar Košī spriegumu \(\boldsymbol{\sigma}\) un atskaites konfigurācijas gradientu \(\partial N_\alpha^e/\partial X_i\) ar pašreizējās konfigurācijas gradientu \(\partial N_\alpha^e/\partial x_i\).

\(\boldsymbol{G}\) ir no Košī sprieguma atkarīga korekcijas matrica, kas nepieciešama, lai hipoelastīgais konstitutīvais likums \(\hat{\boldsymbol{\sigma}}^{\nabla J} = \boldsymbol{\mathsf{C}}:\boldsymbol{D}\) pieskares stinguma shēmā būtu saderīgs ar uz Trūsdela ātrumu balstītu konstitutīvo likumu. To iegūst, sakārtojot ceturtās kārtas tenzora komponentes \(G_{ijkl} = \delta_{il}\sigma_{kj} + \delta_{jl}\sigma_{ik}\) \(6\times 6\) Foigta formā:

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

Šī ir iegūtā matrica.

Globālās stinguma matricas montāža¶

Globālo pieskares stingumu \(\boldsymbol{K}\) iegūst, sadalot katra elementa stingumu \(\boldsymbol{K}^e\) \(d\times d\) blokos \(\boldsymbol{K}^e_{\alpha\beta}\) katram mezglu pārim un izmantojot otrās kārtas tenzora montāžas kopu \(\mathcal{E}^2(i_g, i_h)\), kas ieviesta sadaļā Elementu mezglu fizikālo lielumu montāža:

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

Iegūtās vērtības sakārto matricā ar rindu \(i_g\) un kolonnu \(i_h\). Realizācijā kopa \(\mathcal{E}^2\) netiek izveidota tieši; tā vietā atbilstošos blokus pievieno tieši elementu ciklā. Matrica ir kvadrātiska, un tās dimensija ir brīvības pakāpju skaits uz mezglu \(\times\) kopējais mezglu skaits \(n_g\), taču, tā kā komponentes starp mezgliem, kas nav savienoti ar elementiem, ir \(0\), to glabā retas matricas formā.

TL un UL formulējumu elementu stinguma matricām ir vienāda forma, izņemot atskaites konfigurācijas pārslēgšanu (mezglu koordinātas un B matricas izveidei izmantotais avots) un \(\boldsymbol{G}\) matricas esamību vai neesamību. Tādēļ FrontISTR abus formulējumus realizē kopīgā apakšprogrammā.

Iterācijas algoritms¶

Apkopojot iepriekš minēto, iterācijas sākumā iestata \(\Delta\boldsymbol{u} = \boldsymbol{0}\) un aprēķina sākotnējo atlikumu \(\boldsymbol{R}_0 = \boldsymbol{F}_{n+1} - \boldsymbol{Q}(\boldsymbol{u}_n)\). Tad \(i\)-tajā iterācijā izpilda šādu procedūru.

Pie pašreizējā pārvietojuma \(\boldsymbol{u}_n + \Delta\boldsymbol{u}\) aprēķina pieskares stingumu \(\boldsymbol{K}_i = \partial\boldsymbol{Q}/\partial\boldsymbol{u}|_{\boldsymbol{u}_n + \Delta\boldsymbol{u}}\), izmantojot sadaļā Pieskares stinguma matricas izveide aprakstīto procedūru.
Lai uzliktu ģeometriskos robežnosacījumus, modificē pieskares stinguma matricu un atlikuma vektoru brīvības pakāpēm, uz kurām attiecas pārvietojuma ierobežojumi, iegūstot \(\tilde{\boldsymbol{K}}_i, \tilde{\boldsymbol{R}}_{i-1}\) (Ģeometrisko robežnosacījumu apstrāde).
Atrisina lineāro vienādojumu \(\tilde{\boldsymbol{K}}_i\, d\boldsymbol{u}_i = \tilde{\boldsymbol{R}}_{i-1}\), lai iegūtu korekciju \(d\boldsymbol{u}_i\). Šī procedūra bieži veido lielāko daļu iteratīvā aprēķina skaitļošanas izmaksu.
Atjaunina pārvietojuma pieaugumu kā \(\Delta\boldsymbol{u} \leftarrow \Delta\boldsymbol{u} + d\boldsymbol{u}_i\) un attiecīgi aprēķina iekšējo spēku vektoru \(\boldsymbol{Q}(\boldsymbol{u}_n + \Delta\boldsymbol{u})\) un atlikumu \(\boldsymbol{R}_i = \boldsymbol{F}_{n+1} - \boldsymbol{Q}(\boldsymbol{u}_n + \Delta\boldsymbol{u})\).
Pārbauda konverģenci un, ja tā sasniegta, iterāciju pabeidz. Komponentes, kas atbilst ierobežojumu reakcijām, parādās atlikumā \(\boldsymbol{R}_i\) brīvības pakāpēm ar ģeometriskajiem robežnosacījumiem, tādēļ konverģences indikatoru veido no \(\tilde{\boldsymbol{R}}_i\), izslēdzot šīs komponentes. Konkrēti konverģences indikatori un sliekšņi aprakstīti sadaļā Konverģences kritēriji. Ja konverģence nav sasniegta un tiek sasniegts iterāciju ierobežojums, iterāciju uzskata par neveiksmīgu.

Kad iterācija konverģē, konverģējušo \(\Delta\boldsymbol{u}\) pieskaita \(\boldsymbol{u}_n\), iegūstot uzkrāto pārvietojumu laikā \(t_{n+1}\), \(\boldsymbol{u}_{n+1} = \boldsymbol{u}_n + \Delta\boldsymbol{u}\), un pāriet uz nākamo laika soli.

Saistītās tēmas¶

Ārējo spēku virtuālais darbs un globālo vienādojumu montāža — risināmā nelineārā vienādojuma sākumpunkts
Iekšējā virtuālā darba diskretizācija — \(\boldsymbol{B}_L, \boldsymbol{B}_{NL}, \boldsymbol{b}\) izveide
Ģeometrisko robežnosacījumu apstrāde — pieskares stinguma matricas un atlikuma vektora modificēšana, lai uzliktu pārvietojuma ierobežojumus
Konverģences kritēriji — apstāšanās kritēriji, kas balstīti uz atlikuma normu
Tenzoru pieraksts un matemātiskie pamati — materiāla matricas \(\tilde{\boldsymbol{C}}\) Foigta attēlojums
Nelineārās iterācijas un laika integrēšana (funkcijas) — lietošana un izvēle funkciju atsaucē