Gå till innehållet

Newton-Raphson-metoden

Linjärisering och iterativ rekursion

I Yttre krafters virtuella arbete och sammansättning 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 det virtuella arbetet av inre krafter, 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 andra 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 första termen i högerledet är materialstyvhetstermen (initialförskjutningstermen), och den andra 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 det virtuella arbetet av inre krafter, 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 andra 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 första termen i högerledet är materialstyvhetstermen och den andra 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 det virtuella arbetet av inre krafter). \(\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 andra 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 fjärde 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.

Sammansättning 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 sammansättning av andra ordningens tensorer \(\mathcal{E}^2(i_g, i_h)\) som infördes i Sammansättning 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

AI-assisted translation May contain errors Official docs Status