Saltar a contenido

Método de Newton-Raphson

Linealización y recurrencia iterativa

La formulación del trabajo virtual de las fuerzas externas y el ensamblaje de la ecuación global proporciona una ecuación no lineal para el desplazamiento nodal \(\boldsymbol{u}_{n+1}\) en el instante \(t_{n+1}\), que se resuelve mediante el método de Newton-Raphson. Se supone conocido el desplazamiento nodal \(\boldsymbol{u}_n\) hasta el instante \(t_n\), y se toma el incremento de desplazamiento \(\Delta\boldsymbol{u}\) como variable desconocida para determinar

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

En lo sucesivo se desprecia la dependencia del vector de fuerzas externas respecto del desplazamiento nodal y, tomando \(\boldsymbol{F}(\boldsymbol{u}_{n+1}) = \boldsymbol{F}_{n+1}\), se resuelve

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

En la solución actual \(\Delta\boldsymbol{u}\) se define la rigidez tangente

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

Con ella, la linealización de la ecuación no lineal da

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

Si la corrección de la iteración \(i\)-ésima se denota por \(d\boldsymbol{u}_i\) y el vector residual al inicio de la iteración por

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

la recurrencia iterativa es

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

El residual \(\boldsymbol{R}_i\) es, por tanto, una magnitud que representa el desequilibrio de fuerzas respecto del equilibrio.

Construcción de la matriz de rigidez tangente

La rigidez tangente \(\boldsymbol{K} = \partial\boldsymbol{Q}/\partial\boldsymbol{u}\) se construye derivando parcialmente respecto del desplazamiento nodal el vector de fuerzas internas del elemento obtenido en la discretización del trabajo virtual de las fuerzas internas, integrando los integrandos resultantes a nivel de elemento sobre cada dominio elemental y ensamblándolos. Si se denota el integrando a nivel de elemento por \(\boldsymbol{K}^e_X\) (notación en la configuración de referencia, formulación TL) o por \(\boldsymbol{K}^e_x\) (notación en la configuración actual, formulación UL), la rigidez tangente del elemento es

\[ \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ón se muestran las formas finales de los integrandos TL/UL. En ambos casos se descomponen como suma de un término de rigidez material (término de desplazamiento inicial) y un término de rigidez geométrica (término de tensión inicial).

Formulación Total Lagrange

En la formulación Total Lagrange se supone una relación lineal entre la tasa de la segunda tensión de Piola-Kirchhoff \(\dot{\boldsymbol{S}}\) y la tasa de deformación de Green-Lagrange \(\dot{\boldsymbol{E}}\), es decir, \(\dot{\boldsymbol{S}} = \boldsymbol{\mathsf{C}}:\dot{\boldsymbol{E}}\). Esto corresponde a leyes constitutivas de materiales elásticos lineales (materiales de St. Venant-Kirchhoff) y materiales hiperelásticos; FrontISTR utiliza la formulación Total Lagrange para estos materiales. El integrando de la rigidez tangente del elemento se escribe entonces, en forma tensorial, como

\[ \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 término del lado derecho es el término de rigidez material (término de desplazamiento inicial), y el segundo es el término de rigidez geométrica (término de tensión inicial).

En la implementación de FrontISTR, este integrando se calcula en forma matricial utilizando la notación 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 \]

Las matrices son las siguientes. \(\boldsymbol{B}_L, \boldsymbol{B}_{NL}\) son las matrices B introducidas en la discretización del trabajo virtual de las fuerzas internas, y \(\tilde{\boldsymbol{C}}\) es la representación de Voigt del tensor constitutivo \(\boldsymbol{\mathsf{C}}\), es decir, una matriz de rigidez material de \(6\times 6\) (notación tensorial y fundamentos matemáticos). \(\boldsymbol{S}_9, \boldsymbol{F}_9\) son las siguientes matrices de reordenación utilizadas para expresar el término de rigidez geométrica como un producto matricial. Primero, para un tensor de segundo orden \(\boldsymbol{A}\) de \(3\times 3\), se define la notación \([\,\cdot\,]\) que lo reordena en un vector de 9 componentes como

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

Con esta definición, \(\boldsymbol{F}_9\) expresa la variación del gradiente de deformación en la forma \([\delta\boldsymbol{F}] = \boldsymbol{F}_9\, \delta\boldsymbol{u}^e\) y es una matriz de \(9\times d n_e\). Para el nodo de elemento \(\alpha = 1, \ldots, n_e\), el bloque correspondiente de \(9\times d\) es

\[ [\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{ es la matriz identidad de } 3\times 3) \]

y \(\boldsymbol{F}_9 = [[\boldsymbol{F}_9]_1, \ldots, [\boldsymbol{F}_9]_{n_e}]\) se obtiene disponiendo horizontalmente los bloques en el orden de los nodos del elemento. \(\boldsymbol{S}_9\) se elige de modo que, al combinarla con esta matriz, el término de rigidez geométrica se exprese como \(\delta\boldsymbol{u}^{eT}\, \boldsymbol{F}_9^T \boldsymbol{S}_9 \boldsymbol{F}_9\, \dot{\boldsymbol{u}}^e\); es la matriz de \(9\times 9\)

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

Formulación Updated Lagrange

En la formulación Updated Lagrange se supone una relación lineal entre la tasa de Jaumann del tensor de tensión relativa de Kirchhoff \(\hat{\boldsymbol{\sigma}}^{\nabla J}\) y el tensor de velocidad de deformación \(\boldsymbol{D}\), es decir, \(\hat{\boldsymbol{\sigma}}^{\nabla J} = \boldsymbol{\mathsf{C}}:\boldsymbol{D}\). Esta es la forma de una ley constitutiva hipoelástica común a materiales elásticos lineales, elastoplásticos y de creep; FrontISTR utiliza la formulación Updated Lagrange para estos materiales. El integrando de la rigidez tangente del elemento expresado en la configuración actual se escribe entonces, en forma tensorial, como

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

donde \(\boldsymbol{\sigma}^{\nabla T}\) es la tasa de Truesdell, \(\boldsymbol{A}_{(L)}\) es la parte lineal de la deformación de Almansi, \(\boldsymbol{F}_t = \partial\boldsymbol{u}/\partial\boldsymbol{x}\) es el gradiente de desplazamiento respecto de la configuración actual y \(\boldsymbol{L}\) es el tensor gradiente de velocidad. El primer término del lado derecho es el término de rigidez material, y el segundo es el término de rigidez geométrica.

En la implementación de FrontISTR, este integrando se calcula en forma matricial utilizando la notación 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}\) es la matriz B construida en la configuración actual (discretización del trabajo virtual de las fuerzas internas). \(\boldsymbol{\sigma}_9, \boldsymbol{f}_9\) se obtienen a partir de \(\boldsymbol{S}_9, \boldsymbol{F}_9\), definidas para la formulación TL, sustituyendo la segunda tensión PK \(\boldsymbol{S}\) por la tensión de Cauchy \(\boldsymbol{\sigma}\) y el gradiente de la configuración de referencia \(\partial N_\alpha^e/\partial X_i\) por el gradiente de la configuración actual \(\partial N_\alpha^e/\partial x_i\).

\(\boldsymbol{G}\) es una matriz de corrección dependiente de la tensión de Cauchy necesaria para hacer compatible la ley constitutiva hipoelástica \(\hat{\boldsymbol{\sigma}}^{\nabla J} = \boldsymbol{\mathsf{C}}:\boldsymbol{D}\) con el marco de rigidez tangente como ley constitutiva basada en la tasa de Truesdell. Se obtiene disponiendo los componentes del tensor de cuarto orden \(G_{ijkl} = \delta_{il}\sigma_{kj} + \delta_{jl}\sigma_{ik}\) en la representación de Voigt de \(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} \]

Ensamblaje de la matriz de rigidez global

La rigidez tangente global \(\boldsymbol{K}\) se obtiene dividiendo cada rigidez de elemento \(\boldsymbol{K}^e\) en bloques de \(d\times d\), \(\boldsymbol{K}^e_{\alpha\beta}\), para cada par de nodos y utilizando el conjunto de ensamblaje de tensores de segundo orden \(\mathcal{E}^2(i_g, i_h)\) introducido en el ensamblaje de magnitudes físicas nodales de elemento:

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

Los valores resultantes se disponen como una matriz con fila \(i_g\) y columna \(i_h\). En la implementación, el conjunto \(\mathcal{E}^2\) no se construye explícitamente; en su lugar, los bloques correspondientes se suman directamente dentro del bucle de elementos. La matriz es cuadrada, con dimensión igual a grados de libertad por nodo \(\times\) número total de nodos \(n_g\); como los componentes distintos de los correspondientes a nodos conectados mediante elementos son \(0\), se almacena en formato de matriz dispersa.

Las matrices de rigidez de elemento de las formulaciones TL y UL tienen la misma forma salvo por el cambio de configuración de referencia (coordenadas nodales y origen empleado para construir la matriz B) y por la presencia o ausencia de la matriz \(\boldsymbol{G}\). Por ello, FrontISTR implementa ambas formulaciones en una subrutina común.

Algoritmo de iteración

En resumen, al iniciar la iteración se fija \(\Delta\boldsymbol{u} = \boldsymbol{0}\) y se calcula el residual inicial \(\boldsymbol{R}_0 = \boldsymbol{F}_{n+1} - \boldsymbol{Q}(\boldsymbol{u}_n)\). A continuación, en la iteración \(i\)-ésima se ejecuta el procedimiento siguiente.

  1. Para el desplazamiento actual \(\boldsymbol{u}_n + \Delta\boldsymbol{u}\), se calcula la rigidez tangente \(\boldsymbol{K}_i = \partial\boldsymbol{Q}/\partial\boldsymbol{u}|_{\boldsymbol{u}_n + \Delta\boldsymbol{u}}\) mediante el procedimiento de construcción de la matriz de rigidez tangente.
  2. Para imponer las condiciones de contorno geométricas, se modifican la matriz de rigidez tangente y el vector residual en los grados de libertad sujetos a restricciones de desplazamiento, obteniéndose \(\tilde{\boldsymbol{K}}_i, \tilde{\boldsymbol{R}}_{i-1}\) (tratamiento de las condiciones de contorno geométricas).
  3. Se resuelve la ecuación lineal \(\tilde{\boldsymbol{K}}_i\, d\boldsymbol{u}_i = \tilde{\boldsymbol{R}}_{i-1}\) para obtener la corrección \(d\boldsymbol{u}_i\). Este procedimiento suele representar la mayor parte del coste computacional del cálculo iterativo.
  4. Se actualiza el incremento de desplazamiento como \(\Delta\boldsymbol{u} \leftarrow \Delta\boldsymbol{u} + d\boldsymbol{u}_i\) y, en consecuencia, se calculan el vector de fuerzas internas \(\boldsymbol{Q}(\boldsymbol{u}_n + \Delta\boldsymbol{u})\) y el residual \(\boldsymbol{R}_i = \boldsymbol{F}_{n+1} - \boldsymbol{Q}(\boldsymbol{u}_n + \Delta\boldsymbol{u})\).
  5. Se comprueba la convergencia y se termina la iteración si se ha alcanzado. En los grados de libertad sujetos a condiciones de contorno geométricas aparecen en el residual \(\boldsymbol{R}_i\) componentes correspondientes a reacciones de restricción; por ello, el indicador de convergencia se construye a partir de \(\tilde{\boldsymbol{R}}_i\), excluyendo dichos componentes. Los indicadores y umbrales de convergencia concretos se describen en criterios de convergencia. Si no se alcanza la convergencia y se llega al límite de iteraciones, la iteración se considera fallida.

Cuando la iteración converge, el valor convergido de \(\Delta\boldsymbol{u}\) se suma a \(\boldsymbol{u}_n\) para obtener el desplazamiento acumulado en el instante \(t_{n+1}\), \(\boldsymbol{u}_{n+1} = \boldsymbol{u}_n + \Delta\boldsymbol{u}\), y se avanza al siguiente paso de tiempo.

Temas relacionados