Gå till innehållet

Newton-Raphson-metoden

Linjärisering och iterativ rekursion

I Yttre krafters virtuella arbete och assemblering av det globala ekvationssystemet erhålls en icke-linjär ekvation för nodförskjutningen \(\boldsymbol{u}_{n+1}\) vid tiden \(t_{n+1}\), vilken löses med Newton-Raphson-metoden. Nodförskjutningen \(\boldsymbol{u}_n\) fram till tiden \(t_n\) antas vara känd, och förskjutningsinkrementet \(\Delta\boldsymbol{u}\) tas som den okända variabel som ska bestämmas

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

I fortsättningen försummas den yttre kraftvektorns beroende av nodförskjutningen, och med \(\boldsymbol{F}(\boldsymbol{u}_{n+1}) = \boldsymbol{F}_{n+1}\)

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

ska lösas.

Vid den aktuella lösningen \(\Delta\boldsymbol{u}\) definieras tangentstyvheten

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

Med denna ger linjärisering av den icke-linjära ekvationen

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

Låt korrektionen i den \(i\):te iterationen vara \(d\boldsymbol{u}_i\), och definiera residualvektorn vid iterationens början som

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

Då blir den iterativa rekursionen

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

Residualen \(\boldsymbol{R}_i\) är alltså en storhet som motsvarar kraftobalansen från jämvikt.

Konstruktion av tangentstyvhetsmatrisen

Tangentstyvheten \(\boldsymbol{K} = \partial\boldsymbol{Q}/\partial\boldsymbol{u}\) konstrueras genom att partiellt derivera elementets inre kraftvektor, som erhölls i Diskretisering av inre virtuellt arbete, med avseende på nodförskjutningen, integrera de resulterande integranderna på elementnivå över respektive elementdomän och därefter sammanställa dem. Om integranden på elementnivå betecknas \(\boldsymbol{K}^e_X\) (notation i referenskonfigurationen, TL-formulering) eller \(\boldsymbol{K}^e_x\) (notation i aktuell konfiguration, UL-formulering), är elementets tangentstyvhet

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

Nedan ges slutformerna för TL-/UL-integranderna. I båda fallen delas de upp i summan av en materialstyvhetsterm (initialförskjutningsterm) och en geometrisk styvhetsterm (initialspänningsterm).

Total-Lagrange-formulering

I total-Lagrange-formuleringen antas ett linjärt samband mellan hastigheten för den 2:a Piola-Kirchhoff-spänningen \(\dot{\boldsymbol{S}}\) och Green-Lagrange-töjningshastigheten \(\dot{\boldsymbol{E}}\), nämligen \(\dot{\boldsymbol{S}} = \boldsymbol{\mathsf{C}}:\dot{\boldsymbol{E}}\). Detta motsvarar konstitutiva lagar för linjärelastiska material (St. Venant-Kirchhoff-material) och hyperelastiska material, och FrontISTR använder total-Lagrange-formuleringen för dessa material. Integranden för elementets tangentstyvhet skrivs då i tensorform som

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

Den 1:a termen i högerledet är materialstyvhetstermen (initialförskjutningstermen), och den 2:a termen är den geometriska styvhetstermen (initialspänningstermen).

I FrontISTR-implementationen utvärderas denna integrand i matrisform med Voigt-notation:

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

Matriserna är följande. \(\boldsymbol{B}_L, \boldsymbol{B}_{NL}\) är B-matriserna som infördes i Diskretisering av inre virtuellt arbete, och \(\tilde{\boldsymbol{C}}\) är Voigt-representationen av den konstitutiva tensorn \(\boldsymbol{\mathsf{C}}\), det vill säga en \(6\times 6\) materialstyvhetsmatris (Tensornotation och matematiska grunder). \(\boldsymbol{S}_9, \boldsymbol{F}_9\) är följande omordningsmatriser som används för att uttrycka den geometriska styvhetstermen som en matrisprodukt. För en 2:a ordningens \(3\times 3\)-tensor \(\boldsymbol{A}\) definieras först notation \([\,\cdot\,]\), som ordnar om den till en vektor med 9 komponenter, som

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

Med denna definition uttrycker \(\boldsymbol{F}_9\) variationen av deformationsgradienten i formen \([\delta\boldsymbol{F}] = \boldsymbol{F}_9\, \delta\boldsymbol{u}^e\) och är en \(9\times d n_e\)-matris. För elementnod \(\alpha = 1, \ldots, n_e\) är motsvarande \(9\times d\)-block

\[ [\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{ är } 3\times 3 \text{ identitetsmatrisen}) \]

och \(\boldsymbol{F}_9 = [[\boldsymbol{F}_9]_1, \ldots, [\boldsymbol{F}_9]_{n_e}]\) erhålls genom att blocken ordnas horisontellt i elementnodernas ordning. \(\boldsymbol{S}_9\) väljs så att den geometriska styvhetstermen tillsammans med denna matris uttrycks som \(\delta\boldsymbol{u}^{eT}\, \boldsymbol{F}_9^T \boldsymbol{S}_9 \boldsymbol{F}_9\, \dot{\boldsymbol{u}}^e\); det är följande \(9\times 9\)-matris

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

Detta är den resulterande matrisen.

Updated-Lagrange-formulering

I updated-Lagrange-formuleringen antas ett linjärt samband mellan Jaumann-hastigheten för den relativa Kirchhoff-spänningstensorn \(\hat{\boldsymbol{\sigma}}^{\nabla J}\) och deformationshastighetstensor \(\boldsymbol{D}\), nämligen \(\hat{\boldsymbol{\sigma}}^{\nabla J} = \boldsymbol{\mathsf{C}}:\boldsymbol{D}\). Detta är formen för en hypoelastisk konstitutiv lag som är gemensam för linjärelastiska, elastoplastiska och krypmaterial, och FrontISTR använder updated-Lagrange-formuleringen för dessa material. Integranden för elementets tangentstyvhet uttryckt i den aktuella konfigurationen skrivs då i tensorform som

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

där \(\boldsymbol{\sigma}^{\nabla T}\) är Truesdell-hastigheten, \(\boldsymbol{A}_{(L)}\) är den linjära delen av Almansi-töjningen, \(\boldsymbol{F}_t = \partial\boldsymbol{u}/\partial\boldsymbol{x}\) är förskjutningsgradienten med avseende på den aktuella konfigurationen och \(\boldsymbol{L}\) är hastighetsgradienttensorn. Den 1:a termen i högerledet är materialstyvhetstermen och den 2:a termen är den geometriska styvhetstermen.

I FrontISTR-implementationen utvärderas denna integrand i matrisform med Voigt-notation:

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

Här är \(\boldsymbol{b}\) B-matrisen som konstrueras i den aktuella konfigurationen (Diskretisering av inre virtuellt arbete). \(\boldsymbol{\sigma}_9, \boldsymbol{f}_9\) erhålls från \(\boldsymbol{S}_9, \boldsymbol{F}_9\), definierade för TL-formuleringen, genom att ersätta den 2:a PK-spänningen \(\boldsymbol{S}\) med Cauchy-spänningen \(\boldsymbol{\sigma}\) och gradienten i referenskonfigurationen \(\partial N_\alpha^e/\partial X_i\) med gradienten i den aktuella konfigurationen \(\partial N_\alpha^e/\partial x_i\).

\(\boldsymbol{G}\) är en Cauchy-spänningsberoende korrektionsmatris som krävs för att göra den hypoelastiska konstitutiva lagen \(\hat{\boldsymbol{\sigma}}^{\nabla J} = \boldsymbol{\mathsf{C}}:\boldsymbol{D}\) förenlig med tangentstyvhetsramverket som en Truesdell-hastighetsbaserad konstitutiv lag. Den erhålls genom att ordna 4:e ordningens tensorkomponenter \(G_{ijkl} = \delta_{il}\sigma_{kj} + \delta_{jl}\sigma_{ik}\) i \(6\times 6\) Voigt-form som

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

Detta är den resulterande matrisen.

Assemblering av den globala styvhetsmatrisen

Den globala tangentstyvheten \(\boldsymbol{K}\) erhålls genom att dela varje elementstyvhet \(\boldsymbol{K}^e\) i \(d\times d\)-block \(\boldsymbol{K}^e_{\alpha\beta}\) för varje nodpar och använda mängden för assemblering av 2:a ordningens tensorer \(\mathcal{E}^2(i_g, i_h)\) som infördes i Assemblering av elementens nodala fysikaliska storheter:

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

De erhållna värdena ordnas som en matris med rad \(i_g\) och kolumn \(i_h\). I implementationen konstrueras inte mängden \(\mathcal{E}^2\) explicit; i stället adderas motsvarande block direkt inne i elementloopen. Matrisen är kvadratisk med en dimension lika med antalet frihetsgrader per nod \(\times\) det totala antalet noder \(n_g\), men eftersom komponenter utöver dem mellan noder som är kopplade via element är \(0\) lagras den i gles matrisform.

Elementstyvhetsmatriserna för TL- och UL-formuleringarna har samma form förutom byte av referenskonfiguration (nodkoordinater och den källa som används för att konstruera B-matrisen) samt förekomst eller frånvaro av \(\boldsymbol{G}\)-matrisen. FrontISTR implementerar därför båda formuleringarna i en gemensam subrutin.

Iterationsalgoritm

Sammanfattningsvis sätts \(\Delta\boldsymbol{u} = \boldsymbol{0}\) vid iterationens början och den initiala residualen \(\boldsymbol{R}_0 = \boldsymbol{F}_{n+1} - \boldsymbol{Q}(\boldsymbol{u}_n)\) beräknas. I den \(i\):te iterationen utförs därefter följande procedur.

  1. Vid den aktuella förskjutningen \(\boldsymbol{u}_n + \Delta\boldsymbol{u}\) beräknas tangentstyvheten \(\boldsymbol{K}_i = \partial\boldsymbol{Q}/\partial\boldsymbol{u}|_{\boldsymbol{u}_n + \Delta\boldsymbol{u}}\) med proceduren i Konstruktion av tangentstyvhetsmatrisen.
  2. För att införa geometriska randvillkor modifieras tangentstyvhetsmatrisen och residualvektorn för de frihetsgrader som omfattas av förskjutningsvillkor, vilket ger \(\tilde{\boldsymbol{K}}_i, \tilde{\boldsymbol{R}}_{i-1}\) (Behandling av geometriska randvillkor).
  3. Lös den linjära ekvationen \(\tilde{\boldsymbol{K}}_i\, d\boldsymbol{u}_i = \tilde{\boldsymbol{R}}_{i-1}\) för att erhålla korrektionen \(d\boldsymbol{u}_i\). Denna procedur står ofta för större delen av beräkningskostnaden i iterationsberäkningen.
  4. Uppdatera förskjutningsinkrementet som \(\Delta\boldsymbol{u} \leftarrow \Delta\boldsymbol{u} + d\boldsymbol{u}_i\) och beräkna därefter den inre kraftvektorn \(\boldsymbol{Q}(\boldsymbol{u}_n + \Delta\boldsymbol{u})\) samt residualen \(\boldsymbol{R}_i = \boldsymbol{F}_{n+1} - \boldsymbol{Q}(\boldsymbol{u}_n + \Delta\boldsymbol{u})\).
  5. Kontrollera konvergensen och avsluta iterationen om konvergens har uppnåtts. Komponenter som motsvarar inspänningsreaktioner uppträder i residualen \(\boldsymbol{R}_i\) för frihetsgrader som omfattas av geometriska randvillkor, så konvergensindikatorn konstrueras från \(\tilde{\boldsymbol{R}}_i\) efter att dessa komponenter har uteslutits. Konkreta konvergensindikatorer och tröskelvärden beskrivs i Konvergenskriterier. Om konvergens inte uppnås och iterationsgränsen nås behandlas iterationen som misslyckad.

När iterationen konvergerar adderas det konvergerade \(\Delta\boldsymbol{u}\) till \(\boldsymbol{u}_n\) för att erhålla den ackumulerade förskjutningen vid tiden \(t_{n+1}\), \(\boldsymbol{u}_{n+1} = \boldsymbol{u}_n + \Delta\boldsymbol{u}\), och beräkningen går vidare till nästa tidssteg.

Relaterade ämnen