Ir para o conteúdo

Método de Newton-Raphson

Linearização e recorrência iterativa

A equação não linear para o deslocamento nodal no instante \(t_{n+1}\), \(\boldsymbol{u}_{n+1}\), obtida em Trabalho virtual das forças externas e montagem da equação global, é resolvida pelo método de Newton-Raphson. Considera-se conhecido o deslocamento nodal até o instante \(t_n\), \(\boldsymbol{u}_n\), e toma-se o incremento de deslocamento \(\Delta\boldsymbol{u}\) como variável desconhecida para determinar

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

A partir daqui, despreza-se a dependência do vetor de forças externas em relação ao deslocamento nodal e, tomando \(\boldsymbol{F}(\boldsymbol{u}_{n+1}) = \boldsymbol{F}_{n+1}\),

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

é resolvida.

Na solução atual \(\Delta\boldsymbol{u}\), define-se a rigidez tangente

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

Usando-a, a linearização da equação não linear fornece

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

Na \(i\)-ésima iteração, seja a correção \(d\boldsymbol{u}_i\) e seja o vetor residual no início da iteração

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

Então, a recorrência iterativa é

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

Dessa forma, o residual \(\boldsymbol{R}_i\) é uma grandeza correspondente ao desequilíbrio de forças em relação ao equilíbrio.

Construção da matriz de rigidez tangente

A rigidez tangente \(\boldsymbol{K} = \partial\boldsymbol{Q}/\partial\boldsymbol{u}\) é construída diferenciando parcialmente, em relação ao deslocamento nodal, o vetor de forças internas do elemento obtido em Discretização do trabalho virtual das forças internas, integrando os integrandos resultantes de cada elemento sobre o domínio do elemento e reunindo-os. Denotando o integrando no nível do elemento por \(\boldsymbol{K}^e_X\) (notação na configuração de referência, formulação TL) ou \(\boldsymbol{K}^e_x\) (notação na configuração atual, formulação UL), a rigidez tangente do elemento é

\[ \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 seguir são apresentadas as formas finais dos integrandos TL/UL. Em ambos os casos, elas se decompõem na soma de um termo de rigidez do material (termo de deslocamento inicial) e um termo de rigidez geométrica (termo de tensão inicial).

Formulação Total Lagrange

Na formulação Total Lagrange, assume-se uma relação linear entre a taxa da segunda tensão de Piola-Kirchhoff \(\dot{\boldsymbol{S}}\) e a taxa da deformação de Green-Lagrange \(\dot{\boldsymbol{E}}\), isto é, \(\dot{\boldsymbol{S}} = \boldsymbol{\mathsf{C}}:\dot{\boldsymbol{E}}\). Isso corresponde às leis constitutivas de materiais elásticos lineares (materiais de St. Venant-Kirchhoff) e materiais hiperelásticos, e o FrontISTR utiliza a formulação Total Lagrange para esses materiais. Nesse caso, o integrando da rigidez tangente do elemento é escrito, em 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}}) \]

O primeiro termo do lado direito é o termo de rigidez do material (termo de deslocamento inicial), e o segundo termo é o termo de rigidez geométrica (termo de tensão inicial).

Na implementação do FrontISTR, esse integrando é calculado em forma matricial usando a notação 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 matriz é definida a seguir. \(\boldsymbol{B}_L, \boldsymbol{B}_{NL}\) são as matrizes B introduzidas em Discretização do trabalho virtual das forças internas, e \(\tilde{\boldsymbol{C}}\) é a representação de Voigt do tensor constitutivo \(\boldsymbol{\mathsf{C}}\), isto é, uma matriz de rigidez do material \(6\times 6\) (Notação tensorial e fundamentos matemáticos). \(\boldsymbol{S}_9, \boldsymbol{F}_9\) são as matrizes de rearranjo a seguir, usadas para representar o termo de rigidez geométrica como um produto matricial. Primeiro, para um tensor de segunda ordem \(3\times 3\) \(\boldsymbol{A}\), define-se a notação \([\,\cdot\,]\) que o rearranja em um vetor 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 \]

Com essa definição, \(\boldsymbol{F}_9\) representa a variação do gradiente de deformação na forma \([\delta\boldsymbol{F}] = \boldsymbol{F}_9\, \delta\boldsymbol{u}^e\) e é uma matriz \(9\times d n_e\). Para o nó de elemento \(\alpha = 1, \ldots, n_e\), o bloco \(9\times d\) correspondente é

\[ [\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{ é uma matriz } 3\times 3 \text{ identidade}) \]

e \(\boldsymbol{F}_9 = [[\boldsymbol{F}_9]_1, \ldots, [\boldsymbol{F}_9]_{n_e}]\) é obtida organizando esses blocos horizontalmente na ordem dos nós do elemento. \(\boldsymbol{S}_9\) é escolhida de modo que, em combinação com essa matriz, o termo de rigidez geométrica seja expresso como \(\delta\boldsymbol{u}^{eT}\, \boldsymbol{F}_9^T \boldsymbol{S}_9 \boldsymbol{F}_9\, \dot{\boldsymbol{u}}^e\); trata-se da seguinte matriz \(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} \]

Essa é a matriz resultante.

Formulação Updated Lagrange

Na formulação Updated Lagrange, assume-se uma relação linear entre a taxa de Jaumann do tensor de tensão de Kirchhoff relativa \(\hat{\boldsymbol{\sigma}}^{\nabla J}\) e o tensor taxa de deformação \(\boldsymbol{D}\), isto é, \(\hat{\boldsymbol{\sigma}}^{\nabla J} = \boldsymbol{\mathsf{C}}:\boldsymbol{D}\). Essa é a forma de uma lei constitutiva hipoelástica comum a materiais elásticos lineares, elastoplásticos e de fluência, e o FrontISTR utiliza a formulação Updated Lagrange para esses materiais. Nesse caso, o integrando da rigidez tangente do elemento expresso na configuração atual é escrito, em 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}) \]

onde \(\boldsymbol{\sigma}^{\nabla T}\) é a taxa de Truesdell, \(\boldsymbol{A}_{(L)}\) é a parte linear da deformação de Almansi, \(\boldsymbol{F}_t = \partial\boldsymbol{u}/\partial\boldsymbol{x}\) é o gradiente de deslocamento em relação à configuração atual e \(\boldsymbol{L}\) é o tensor gradiente de velocidade. O primeiro termo do lado direito é o termo de rigidez do material e o segundo termo é o termo de rigidez geométrica.

Na implementação do FrontISTR, esse integrando é calculado em forma matricial usando a notação 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 \]

Aqui, \(\boldsymbol{b}\) é a matriz B construída na configuração atual (Discretização do trabalho virtual das forças internas). \(\boldsymbol{\sigma}_9, \boldsymbol{f}_9\) são obtidas a partir de \(\boldsymbol{S}_9, \boldsymbol{F}_9\) definidos para a formulação TL, substituindo a segunda tensão PK \(\boldsymbol{S}\) pela tensão de Cauchy \(\boldsymbol{\sigma}\) e o gradiente na configuração de referência \(\partial N_\alpha^e/\partial X_i\) pelo gradiente na configuração atual \(\partial N_\alpha^e/\partial x_i\).

\(\boldsymbol{G}\) é uma matriz de correção dependente da tensão de Cauchy, necessária para tornar a lei constitutiva hipoelástica \(\hat{\boldsymbol{\sigma}}^{\nabla J} = \boldsymbol{\mathsf{C}}:\boldsymbol{D}\) consistente com o quadro da rigidez tangente como uma lei constitutiva baseada na taxa de Truesdell. Ela é obtida organizando os componentes do tensor de quarta ordem \(G_{ijkl} = \delta_{il}\sigma_{kj} + \delta_{jl}\sigma_{ik}\) na representação de Voigt \(6\times 6\) como

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

Essa é a matriz resultante.

Montagem da matriz de rigidez global

A rigidez tangente global \(\boldsymbol{K}\) é obtida dividindo a rigidez de cada elemento \(\boldsymbol{K}^e\) em blocos \(d\times d\) \(\boldsymbol{K}^e_{\alpha\beta}\) para cada par de nós e usando o conjunto de montagem de tensor de segunda ordem \(\mathcal{E}^2(i_g, i_h)\) introduzido em Montagem de grandezas físicas nodais dos elementos:

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

Os valores resultantes são organizados em uma matriz com linha \(i_g\) e coluna \(i_h\). Na implementação, o conjunto \(\mathcal{E}^2\) não é construído explicitamente; em vez disso, os blocos correspondentes são somados diretamente dentro do loop de elementos. A matriz é quadrada, com dimensão igual ao número de graus de liberdade por nó \(\times\) o número total de nós \(n_g\), mas, como os componentes entre nós que não estão conectados por elementos são \(0\), ela é armazenada em formato de matriz esparsa.

As matrizes de rigidez dos elementos para as formulações TL e UL têm a mesma forma, exceto pela troca da configuração de referência (coordenadas nodais e origem usada para construir a matriz B) e pela presença ou ausência da matriz \(\boldsymbol{G}\). Por isso, o FrontISTR implementa ambas as formulações em uma sub-rotina comum.

Algoritmo de iteração

Resumindo o que foi apresentado, no início da iteração define-se \(\Delta\boldsymbol{u} = \boldsymbol{0}\) e calcula-se o residual inicial \(\boldsymbol{R}_0 = \boldsymbol{F}_{n+1} - \boldsymbol{Q}(\boldsymbol{u}_n)\). Em seguida, na \(i\)-ésima iteração, executa-se o procedimento a seguir.

  1. No deslocamento atual \(\boldsymbol{u}_n + \Delta\boldsymbol{u}\), calcule a rigidez tangente \(\boldsymbol{K}_i = \partial\boldsymbol{Q}/\partial\boldsymbol{u}|_{\boldsymbol{u}_n + \Delta\boldsymbol{u}}\) pelo procedimento de Construção da matriz de rigidez tangente.
  2. Para impor as condições de contorno geométricas, modifique a matriz de rigidez tangente e o vetor residual nos graus de liberdade sujeitos a restrições de deslocamento, obtendo \(\tilde{\boldsymbol{K}}_i, \tilde{\boldsymbol{R}}_{i-1}\) (Tratamento das condições de contorno geométricas).
  3. Resolva a equação linear \(\tilde{\boldsymbol{K}}_i\, d\boldsymbol{u}_i = \tilde{\boldsymbol{R}}_{i-1}\) para obter a correção \(d\boldsymbol{u}_i\). Esse procedimento costuma representar a maior parte do custo computacional do cálculo iterativo.
  4. Atualize o incremento de deslocamento como \(\Delta\boldsymbol{u} \leftarrow \Delta\boldsymbol{u} + d\boldsymbol{u}_i\) e, de acordo com isso, calcule o vetor de forças internas \(\boldsymbol{Q}(\boldsymbol{u}_n + \Delta\boldsymbol{u})\) e o residual \(\boldsymbol{R}_i = \boldsymbol{F}_{n+1} - \boldsymbol{Q}(\boldsymbol{u}_n + \Delta\boldsymbol{u})\).
  5. Verifique a convergência e encerre a iteração se ela for atingida. Nos graus de liberdade sujeitos a condições de contorno geométricas, aparecem no residual \(\boldsymbol{R}_i\) componentes correspondentes às reações de restrição; portanto, o indicador de convergência é construído a partir de \(\tilde{\boldsymbol{R}}_i\), excluindo esses componentes. Os indicadores e limites específicos de convergência são descritos em Critério de convergência. Se a convergência não for atingida e o limite de iterações for alcançado, a iteração é considerada malsucedida.

Quando a iteração converge, some o \(\Delta\boldsymbol{u}\) convergido a \(\boldsymbol{u}_n\) para obter o deslocamento acumulado no instante \(t_{n+1}\), \(\boldsymbol{u}_{n+1} = \boldsymbol{u}_n + \Delta\boldsymbol{u}\), e prossiga para o próximo passo de tempo.

Tópicos relacionados