Salta el contingut

Mètode de Newton-Raphson

Linealització i recurrència iterativa

Treball virtual de les forces externes i assemblatge de l’equació global dona una equació no lineal per al desplaçament nodal en el temps \(t_{n+1}\), \(\boldsymbol{u}_{n+1}\), que es resol mitjançant el mètode de Newton-Raphson. S’assumeix conegut el desplaçament nodal fins al temps \(t_n\), \(\boldsymbol{u}_n\), i es pren l’increment de desplaçament \(\Delta\boldsymbol{u}\) com la variable desconeguda que cal determinar

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

D’ara endavant es negligeix la dependència del vector de forces externes respecte del desplaçament nodal i, amb \(\boldsymbol{F}(\boldsymbol{u}_{n+1}) = \boldsymbol{F}_{n+1}\),

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

es resol.

En la solució actual \(\Delta\boldsymbol{u}\), es defineix la rigidesa tangent

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

Utilitzant-la, la linealització de l’equació no lineal dona

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

Sigui \(d\boldsymbol{u}_i\) la correcció en la iteració \(i\)-èsima, i sigui el vector residual a l’inici de la iteració

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

Aleshores, la recurrència iterativa és

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

Així, el residual \(\boldsymbol{R}_i\) és una magnitud corresponent al desequilibri de forces respecte de l’equilibri.

Construcció de la matriu de rigidesa tangent

La rigidesa tangent \(\boldsymbol{K} = \partial\boldsymbol{Q}/\partial\boldsymbol{u}\) es construeix derivant parcialment respecte del desplaçament nodal el vector de forces internes de l’element obtingut a Discretització del treball virtual de les forces internes, integrant els integrands resultants a nivell d’element sobre cada domini de l’element i assemblant-los. Denotant l’integrand a nivell d’element per \(\boldsymbol{K}^e_X\) (notació de la configuració de referència, formulació TL) o \(\boldsymbol{K}^e_x\) (notació de la configuració actual, formulació UL), la rigidesa tangent de l’element és

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

A continuació es donen les formes finals dels integrands TL/UL. En tots dos casos es descomponen en la suma d’un terme de rigidesa del material (terme de desplaçament inicial) i un terme de rigidesa geomètrica (terme de tensió inicial).

Formulació de Lagrange total

En la formulació de Lagrange total, s’assumeix una relació lineal entre la velocitat de la segona tensió de Piola-Kirchhoff \(\dot{\boldsymbol{S}}\) i la velocitat de deformació de Green-Lagrange \(\dot{\boldsymbol{E}}\), és a dir, \(\dot{\boldsymbol{S}} = \boldsymbol{\mathsf{C}}:\dot{\boldsymbol{E}}\). Això correspon a lleis constitutives per a materials elàstics lineals (materials de St. Venant-Kirchhoff) i materials hiperelàstics, i FrontISTR utilitza la formulació de Lagrange total per a aquests materials. L’integrand de la rigidesa tangent de l’element s’escriu aleshores en forma tensorial com

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

El primer terme del membre dret és el terme de rigidesa del material (terme de desplaçament inicial), i el segon terme és el terme de rigidesa geomètrica (terme de tensió inicial).

En la implementació de FrontISTR, aquest integrand s’avalua en forma matricial utilitzant la notació de 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 \]

Cada matriu és la següent. \(\boldsymbol{B}_L, \boldsymbol{B}_{NL}\) són les matrius B introduïdes a Discretització del treball virtual de les forces internes, i \(\tilde{\boldsymbol{C}}\) és la representació de Voigt del tensor constitutiu \(\boldsymbol{\mathsf{C}}\), és a dir, una matriu de rigidesa del material de \(6\times 6\) (Notació tensorial i fonaments matemàtics). \(\boldsymbol{S}_9, \boldsymbol{F}_9\) són les matrius de reordenació següents utilitzades per expressar el terme de rigidesa geomètrica com un producte matricial. Primer, per a un tensor de segon ordre \(\boldsymbol{A}\) de \(3\times 3\), es defineix la notació \([\,\cdot\,]\) que el reordena en un vector de 9 components com

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

Amb aquesta definició, \(\boldsymbol{F}_9\) expressa la variació del gradient de deformació en la forma \([\delta\boldsymbol{F}] = \boldsymbol{F}_9\, \delta\boldsymbol{u}^e\) i és una matriu de \(9\times d n_e\). Per al node de l’element \(\alpha = 1, \ldots, n_e\), el bloc corresponent de \(9\times d\) és

\[ [\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{ és la matriu identitat de } 3\times 3) \]

i es dona per \(\boldsymbol{F}_9 = [[\boldsymbol{F}_9]_1, \ldots, [\boldsymbol{F}_9]_{n_e}]\), amb els blocs disposats horitzontalment en ordre de nodes de l’element. \(\boldsymbol{S}_9\) es tria de manera que, en combinar-la amb aquesta matriu, el terme de rigidesa geomètrica s’expressi com \(\delta\boldsymbol{u}^{eT}\, \boldsymbol{F}_9^T \boldsymbol{S}_9 \boldsymbol{F}_9\, \dot{\boldsymbol{u}}^e\); és la matriu de \(9\times 9\) següent

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

Aquesta és la matriu resultant.

Formulació de Lagrange actualitzat

En la formulació de Lagrange actualitzat, s’assumeix una relació lineal entre la velocitat de Jaumann del tensor de tensió relativa de Kirchhoff \(\hat{\boldsymbol{\sigma}}^{\nabla J}\) i el tensor de velocitat de deformació \(\boldsymbol{D}\), és a dir, \(\hat{\boldsymbol{\sigma}}^{\nabla J} = \boldsymbol{\mathsf{C}}:\boldsymbol{D}\). Aquesta és la forma d’una llei constitutiva hipoelàstica comuna als materials elàstics lineals, elastoplàstics i de fluència, i FrontISTR utilitza la formulació de Lagrange actualitzat per a aquests materials. L’integrand de la rigidesa tangent de l’element expressat en la configuració actual s’escriu aleshores en forma tensorial com

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

on \(\boldsymbol{\sigma}^{\nabla T}\) és la velocitat de Truesdell, \(\boldsymbol{A}_{(L)}\) és la part lineal de la deformació d’Almansi, \(\boldsymbol{F}_t = \partial\boldsymbol{u}/\partial\boldsymbol{x}\) és el gradient de desplaçament respecte de la configuració actual, i \(\boldsymbol{L}\) és el tensor gradient de velocitat. El primer terme del membre dret és el terme de rigidesa del material, i el segon terme és el terme de rigidesa geomètrica.

En la implementació de FrontISTR, aquest integrand s’avalua en forma matricial utilitzant la notació de 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 \]

Aquí, \(\boldsymbol{b}\) és la matriu B construïda en la configuració actual (Discretització del treball virtual de les forces internes). \(\boldsymbol{\sigma}_9, \boldsymbol{f}_9\) s’obtenen de \(\boldsymbol{S}_9, \boldsymbol{F}_9\) definides per a la formulació TL substituint la segona tensió PK \(\boldsymbol{S}\) per la tensió de Cauchy \(\boldsymbol{\sigma}\) i el gradient de la configuració de referència \(\partial N_\alpha^e/\partial X_i\) pel gradient de la configuració actual \(\partial N_\alpha^e/\partial x_i\).

\(\boldsymbol{G}\) és una matriu de correcció dependent de la tensió de Cauchy necessària perquè la llei constitutiva hipoelàstica \(\hat{\boldsymbol{\sigma}}^{\nabla J} = \boldsymbol{\mathsf{C}}:\boldsymbol{D}\) sigui coherent amb el marc de rigidesa tangent com una llei constitutiva basada en la velocitat de Truesdell. S’obté disposant els components del tensor de quart ordre \(G_{ijkl} = \delta_{il}\sigma_{kj} + \delta_{jl}\sigma_{ik}\) en forma de Voigt de \(6\times 6\) com

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

Aquesta és la matriu resultant.

Assemblatge de la matriu global de rigidesa

La rigidesa tangent global \(\boldsymbol{K}\) s’obté dividint cada rigidesa d’element \(\boldsymbol{K}^e\) en blocs \(\boldsymbol{K}^e_{\alpha\beta}\) de \(d\times d\) per a cada parell de nodes i utilitzant el conjunt d’assemblatge de tensors de segon ordre \(\mathcal{E}^2(i_g, i_h)\) introduït a Assemblatge de magnituds físiques nodals dels elements:

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

Els valors resultants es disposen com una matriu amb fila \(i_g\) i columna \(i_h\). En la implementació, el conjunt \(\mathcal{E}^2\) no es construeix explícitament; en canvi, els blocs corresponents s’afegeixen directament dins del bucle dels elements. La matriu és quadrada, amb una dimensió igual als graus de llibertat per node \(\times\) el nombre total de nodes \(n_g\), però com que els components diferents dels corresponents a nodes connectats mitjançant elements són \(0\), s’emmagatzema en forma de matriu dispersa.

Les matrius de rigidesa dels elements per a les formulacions TL i UL tenen la mateixa forma, excepte pel canvi de configuració de referència (coordenades nodals i font utilitzada per construir la matriu B) i per la presència o absència de la matriu \(\boldsymbol{G}\). Per tant, FrontISTR implementa ambdues formulacions en una subrutina comuna.

Algorisme d’iteració

Resumint l’anterior, a l’inici de la iteració s’estableix \(\Delta\boldsymbol{u} = \boldsymbol{0}\) i es calcula el residual inicial \(\boldsymbol{R}_0 = \boldsymbol{F}_{n+1} - \boldsymbol{Q}(\boldsymbol{u}_n)\). A continuació, en la iteració \(i\)-èsima, es realitza el procediment següent.

  1. En el desplaçament actual \(\boldsymbol{u}_n + \Delta\boldsymbol{u}\), calculeu la rigidesa tangent \(\boldsymbol{K}_i = \partial\boldsymbol{Q}/\partial\boldsymbol{u}|_{\boldsymbol{u}_n + \Delta\boldsymbol{u}}\) mitjançant el procediment de Construcció de la matriu de rigidesa tangent.
  2. Per imposar les condicions de contorn geomètriques, modifiqueu la matriu de rigidesa tangent i el vector residual per als graus de llibertat sotmesos a restriccions de desplaçament, obtenint \(\tilde{\boldsymbol{K}}_i, \tilde{\boldsymbol{R}}_{i-1}\) (Tractament de les condicions de contorn geomètriques).
  3. Resoleu l’equació lineal \(\tilde{\boldsymbol{K}}_i\, d\boldsymbol{u}_i = \tilde{\boldsymbol{R}}_{i-1}\) per obtenir la correcció \(d\boldsymbol{u}_i\). Aquest procediment sovint representa la major part del cost computacional del càlcul iteratiu.
  4. Actualitzeu l’increment de desplaçament com \(\Delta\boldsymbol{u} \leftarrow \Delta\boldsymbol{u} + d\boldsymbol{u}_i\) i, en conseqüència, calculeu el vector de forces internes \(\boldsymbol{Q}(\boldsymbol{u}_n + \Delta\boldsymbol{u})\) i el residual \(\boldsymbol{R}_i = \boldsymbol{F}_{n+1} - \boldsymbol{Q}(\boldsymbol{u}_n + \Delta\boldsymbol{u})\).
  5. Comproveu la convergència i finalitzeu la iteració si s’assoleix. Els components corresponents a reaccions restringides apareixen al residual \(\boldsymbol{R}_i\) en els graus de llibertat subjectes a condicions de contorn geomètriques, de manera que l’indicador de convergència es construeix a partir de \(\tilde{\boldsymbol{R}}_i\) després d’excloure aquests components. Els indicadors i llindars de convergència concrets es descriuen a Criteris de convergència. Si no s’assoleix la convergència i s’arriba al límit d’iteracions, la iteració es considera fallida.

Quan la iteració convergeix, afegiu el \(\Delta\boldsymbol{u}\) convergit a \(\boldsymbol{u}_n\) per obtenir el desplaçament acumulat en el temps \(t_{n+1}\), \(\boldsymbol{u}_{n+1} = \boldsymbol{u}_n + \Delta\boldsymbol{u}\), i continueu amb el pas de temps següent.

Temes relacionats

AI-assisted translation May contain errors Official docs Status