Sari la conținut

Metoda Newton-Raphson

Liniarizare și relația recurentă iterativă

Ecuația neliniară pentru deplasarea nodală \(\boldsymbol{u}_{n+1}\) la timpul \(t_{n+1}\), obținută în Lucrul virtual al forțelor externe și asamblarea ecuației globale, se rezolvă prin metoda Newton-Raphson. Se consideră cunoscută deplasarea nodală \(\boldsymbol{u}_n\) până la timpul \(t_n\), iar incrementul de deplasare \(\Delta\boldsymbol{u}\) este necunoscuta, astfel încât se determină

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

. În continuare se neglijează dependența vectorului forțelor externe de deplasarea nodală și se presupune \(\boldsymbol{F}(\boldsymbol{u}_{n+1}) = \boldsymbol{F}_{n+1}\), rezolvând

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

.

Rigiditatea tangentă la soluția curentă \(\Delta\boldsymbol{u}\) este

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

Folosind-o pentru liniarizarea ecuației neliniare se obține

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

. Dacă corecția la iterația \(i\) este \(d\boldsymbol{u}_i\), iar vectorul rezidual la începutul iterației este

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

atunci relația recurentă iterativă este

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

. Reziduul \(\boldsymbol{R}_i\) reprezintă dezechilibrul de forțe față de starea de echilibru.

Construirea matricei de rigiditate tangentă

Rigiditatea tangentă \(\boldsymbol{K} = \partial\boldsymbol{Q}/\partial\boldsymbol{u}\) se construiește prin derivarea parțială, în raport cu deplasările nodale, a vectorului forțelor interne elementare obținut în Discretizarea lucrului virtual al forțelor interne, integrarea termenului integrand pentru fiecare element pe domeniul elementului și asamblarea rezultatelor. Notând integrandul la nivel de element cu \(\boldsymbol{K}^e_X\) (în configurația de referință, metoda TL) sau \(\boldsymbol{K}^e_x\) (în configurația curentă, metoda UL), rigiditatea tangentă a elementului este

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

. Mai jos sunt prezentate formele finale ale integranzilor TL/UL. Ambele se descompun în suma unui termen de rigiditate materială (termen de deplasare inițială) și a unui termen de rigiditate geometrică (termen de tensiune inițială).

Metoda Total Lagrange

În metoda Total Lagrange se presupune relația liniară \(\dot{\boldsymbol{S}} = \boldsymbol{\mathsf{C}}:\dot{\boldsymbol{E}}\) dintre rata tensiunii a 2-a Piola-Kirchhoff \(\dot{\boldsymbol{S}}\) și rata deformației Green-Lagrange \(\dot{\boldsymbol{E}}\). Aceasta corespunde legilor constitutive pentru materiale liniar elastice (St.Venant-Kirchhoff) și hiperelastice, iar în FrontISTR metoda Total Lagrange este utilizată pentru aceste materiale. Integrandul rigidității tangente elementare se scrie atunci, în formă tensorială, ca

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

. Termenul 1 din membrul drept este termenul de rigiditate materială (termenul de deplasare inițială), iar termenul 2 este termenul de rigiditate geometrică (termenul de tensiune inițială).

În implementarea FrontISTR, acest integrand este calculat în forma matricială cu notație Voigt

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

. Matricele sunt definite astfel. \(\boldsymbol{B}_L, \boldsymbol{B}_{NL}\) sunt matricile B introduse în Discretizarea lucrului virtual al forțelor interne, iar \(\tilde{\boldsymbol{C}}\) este matricea de rigiditate materială \(6\times 6\) a legii constitutive \(\boldsymbol{\mathsf{C}}\) în notație Voigt (Notație tensorială și baze matematice). \(\boldsymbol{S}_9, \boldsymbol{F}_9\) sunt matricele de reordonare de mai jos, utilizate pentru a exprima termenul de rigiditate geometrică drept produs matricial. Mai întâi se definește simbolul \([\,\cdot\,]\), care rearanjează tensorul de ordin 2 \(\boldsymbol{A}\) de dimensiune \(3\times 3\) într-un vector cu 9 componente, prin

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

. Atunci \(\boldsymbol{F}_9\) este o matrice \(9\times d n_e\) care exprimă variația gradientului de deformare în forma \([\delta\boldsymbol{F}] = \boldsymbol{F}_9\, \delta\boldsymbol{u}^e\) și este alcătuită din blocurile \(9\times d\) corespunzătoare nodurilor elementului \(\alpha = 1, \ldots, n_e\):

\[ [\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{ este matricea unitate } 3\times 3) \]

dispuse orizontal în ordinea nodurilor elementului, astfel încât \(\boldsymbol{F}_9 = [[\boldsymbol{F}_9]_1, \ldots, [\boldsymbol{F}_9]_{n_e}]\). \(\boldsymbol{S}_9\) este matricea \(9\times 9\) aleasă astfel încât, împreună cu aceasta, termenul de rigiditate geometrică să poată fi exprimat ca \(\delta\boldsymbol{u}^{eT}\, \boldsymbol{F}_9^T \boldsymbol{S}_9 \boldsymbol{F}_9\, \dot{\boldsymbol{u}}^e\):

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

.

Metoda Updated Lagrange

În metoda Updated Lagrange se presupune relația liniară \(\hat{\boldsymbol{\sigma}}^{\nabla J} = \boldsymbol{\mathsf{C}}:\boldsymbol{D}\) dintre rata Jaumann a tensorului relativ al tensiunii Kirchhoff \(\hat{\boldsymbol{\sigma}}^{\nabla J}\) și tensorul ratei de deformare \(\boldsymbol{D}\). Aceasta este forma unei legi constitutive hipoelastice comune materialelor liniar elastice, elastoplastice și cu fluaj, iar în FrontISTR metoda Updated Lagrange este utilizată pentru aceste materiale. Integrandul rigidității tangente elementare exprimat în configurația curentă se scrie, în formă tensorială, ca

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

. Aici \(\boldsymbol{\sigma}^{\nabla T}\) este rata Truesdell, \(\boldsymbol{A}_{(L)}\) este partea liniară a deformației Almansi, \(\boldsymbol{F}_t = \partial\boldsymbol{u}/\partial\boldsymbol{x}\) este gradientul deplasării raportat la configurația curentă, iar \(\boldsymbol{L}\) este tensorul gradientului vitezei. Termenul 1 din membrul drept este termenul de rigiditate materială, iar termenul 2 este termenul de rigiditate geometrică.

În implementarea FrontISTR, acest integrand este calculat în forma matricială cu notație Voigt

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

. \(\boldsymbol{b}\) este matricea B construită în configurația curentă (Discretizarea lucrului virtual al forțelor interne). \(\boldsymbol{\sigma}_9, \boldsymbol{f}_9\) se obțin din \(\boldsymbol{S}_9, \boldsymbol{F}_9\) definite pentru metoda TL, înlocuind tensiunea a 2-a PK \(\boldsymbol{S}\) cu tensiunea Cauchy \(\boldsymbol{\sigma}\) și gradientul în configurația de referință \(\partial N_\alpha^e/\partial X_i\) cu gradientul în configurația curentă \(\partial N_\alpha^e/\partial x_i\).

\(\boldsymbol{G}\) este matricea de corecție dependentă de tensiunea Cauchy necesară pentru a face legea constitutivă hipoelastică \(\hat{\boldsymbol{\sigma}}^{\nabla J} = \boldsymbol{\mathsf{C}}:\boldsymbol{D}\) compatibilă cu cadrul rigidității tangente ca lege constitutivă bazată pe rata Truesdell; ea reprezintă componentele tensorului de ordin 4 \(G_{ijkl} = \delta_{il}\sigma_{kj} + \delta_{jl}\sigma_{ik}\) în notație Voigt \(6\times 6\):

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

.

Asamblarea matricei globale de rigiditate

Rigiditatea tangentă globală \(\boldsymbol{K}\) se obține împărțind rigiditatea elementară \(\boldsymbol{K}^e\) în blocuri \(d\times d\) pe noduri, \(\boldsymbol{K}^e_{\alpha\beta}\), și folosind mulțimea de asamblare pentru tensori de ordin 2 \(\mathcal{E}^2(i_g, i_h)\) introdusă în Asamblarea mărimilor fizice nodale ale elementelor:

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

plasând rezultatul pe linia \(i_g\) și coloana \(i_h\). În implementare, mulțimea \(\mathcal{E}^2\) nu este construită explicit; blocurile corespunzătoare sunt adunate direct în bucla pe elemente. Matricea este pătratică, de ordinul număr de grade de libertate pe nod \(\times\) numărul total de noduri \(n_g\), dar componentele dintre noduri care nu sunt conectate prin elemente sunt \(0\), astfel încât matricea este păstrată în format rar.

Rigiditățile elementare pentru metodele TL și UL au aceeași formă, cu excepția schimbării configurației de referință (coordonatele nodale și sursa pentru construirea matricei B) și a prezenței sau absenței matricei \(\boldsymbol{G}\); de aceea, în FrontISTR ambele sunt implementate într-o subrutină comună.

Algoritmul iterativ

Rezumând cele de mai sus, la începutul iterației se ia \(\Delta\boldsymbol{u} = \boldsymbol{0}\) și se calculează reziduul inițial \(\boldsymbol{R}_0 = \boldsymbol{F}_{n+1} - \boldsymbol{Q}(\boldsymbol{u}_n)\), iar la iterația \(i\) se execută pașii următori.

  1. Se calculează rigiditatea tangentă \(\boldsymbol{K}_i = \partial\boldsymbol{Q}/\partial\boldsymbol{u}|_{\boldsymbol{u}_n + \Delta\boldsymbol{u}}\) la deplasarea curentă \(\boldsymbol{u}_n + \Delta\boldsymbol{u}\), folosind procedura din Construirea matricei de rigiditate tangentă.
  2. Pentru a aplica condițiile la limită geometrice, matricea de rigiditate tangentă și vectorul rezidual sunt modificați pentru gradele de libertate cu deplasări constrânse, obținând \(\tilde{\boldsymbol{K}}_i, \tilde{\boldsymbol{R}}_{i-1}\) (Tratarea condițiilor la limită geometrice).
  3. Se rezolvă ecuația liniară \(\tilde{\boldsymbol{K}}_i\, d\boldsymbol{u}_i = \tilde{\boldsymbol{R}}_{i-1}\) pentru a determina corecția \(d\boldsymbol{u}_i\). Această operație reprezintă adesea cea mai mare parte a costului de calcul din procesul iterativ.
  4. Se actualizează incrementul de deplasare prin \(\Delta\boldsymbol{u} \leftarrow \Delta\boldsymbol{u} + d\boldsymbol{u}_i\), apoi se calculează corespunzător vectorul forțelor interne \(\boldsymbol{Q}(\boldsymbol{u}_n + \Delta\boldsymbol{u})\) și reziduul \(\boldsymbol{R}_i = \boldsymbol{F}_{n+1} - \boldsymbol{Q}(\boldsymbol{u}_n + \Delta\boldsymbol{u})\).
  5. Se verifică convergența și, dacă este satisfăcută, iterația se încheie. În reziduul \(\boldsymbol{R}_i\), gradele de libertate supuse condițiilor la limită geometrice conțin componente corespunzătoare reacțiunilor de constrângere; de aceea indicatorul de convergență se construiește din componentele \(\tilde{\boldsymbol{R}}_i\) după excluderea acestora. Indicatorii și pragurile concrete sunt descrise în Criteriul de convergență. Dacă nu se obține convergența și se atinge limita de iterații, iterația este considerată eșuată.

După convergență, \(\Delta\boldsymbol{u}\) este adăugat la \(\boldsymbol{u}_n\) pentru a obține deplasarea cumulată la timpul \(t_{n+1}\), \(\boldsymbol{u}_{n+1} = \boldsymbol{u}_n + \Delta\boldsymbol{u}\), apoi calculul trece la următorul pas de timp.

Subiecte conexe