Ir para o conteúdo

Método de Newton-Raphson

Linearização e recorrência iterativa

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

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

Daqui em diante, despreza-se a dependência do vetor de forças externas relativamente ao deslocamento nodal e, com \(\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}} \]

Utilizando-a, a linearização da equação não linear resulta em

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

Seja \(d\boldsymbol{u}_i\) a correção na \(i\)-ésima iteração 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 \]

Assim, o residual \(\boldsymbol{R}_i\) é uma grandeza correspondente ao desequilíbrio de forças relativamente ao equilíbrio.

Construção da matriz de rigidez tangente

A rigidez tangente \(\boldsymbol{K} = \partial\boldsymbol{Q}/\partial\boldsymbol{u}\) é construída derivando 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 ao nível do elemento sobre cada domínio de elemento e efetuando a respetiva montagem. Denotando o integrando ao 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}) \]

Seguem-se as formas finais dos integrandos TL/UL. Em ambos os casos, estes são decompostos na soma de um termo de rigidez do material (termo de deslocamento inicial) e de 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}}\). Isto corresponde a leis constitutivas para materiais elásticos lineares (materiais de St. Venant-Kirchhoff) e materiais hiperelásticos, e o FrontISTR utiliza a formulação Total Lagrange para estes materiais. O integrando da rigidez tangente do elemento é então 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 membro direito é o termo de rigidez do material (termo de deslocamento inicial), e o segundo é o termo de rigidez geométrica (termo de tensão inicial).

Na implementação do FrontISTR, este integrando é avaliado em forma matricial utilizando 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 \]

As matrizes são as seguintes. \(\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 de \(6\times 6\) (Notação tensorial e fundamentos matemáticos). \(\boldsymbol{S}_9, \boldsymbol{F}_9\) são as seguintes matrizes de rearranjo utilizadas para exprimir o termo de rigidez geométrica como um produto matricial. Primeiro, para um tensor de segunda ordem \(\boldsymbol{A}\) de \(3\times 3\), define-se a notação \([\,\cdot\,]\) que o rearranja num 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 esta definição, \(\boldsymbol{F}_9\) exprime a variação do gradiente de deformação na forma \([\delta\boldsymbol{F}] = \boldsymbol{F}_9\, \delta\boldsymbol{u}^e\) e é uma matriz de \(9\times d n_e\). Para o nó de elemento \(\alpha = 1, \ldots, n_e\), o bloco correspondente de \(9\times d\) é

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

e é dada por \(\boldsymbol{F}_9 = [[\boldsymbol{F}_9]_1, \ldots, [\boldsymbol{F}_9]_{n_e}]\), com os blocos dispostos horizontalmente pela ordem dos nós do elemento. \(\boldsymbol{S}_9\) é escolhida de modo que, em combinação com esta 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 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} \]

Esta é 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 relativa de Kirchhoff \(\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}\). Esta é a forma de uma lei constitutiva hipoelástica comum a materiais elásticos lineares, elastoplásticos e com fluência, e o FrontISTR utiliza a formulação Updated Lagrange para estes materiais. O integrando da rigidez tangente do elemento expresso na configuração atual é então 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 relativamente à configuração atual e \(\boldsymbol{L}\) é o tensor gradiente de velocidade. O primeiro termo do membro direito é o termo de rigidez do material e o segundo é o termo de rigidez geométrica.

Na implementação do FrontISTR, este integrando é avaliado em forma matricial utilizando 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\) definidas 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 a estrutura de rigidez tangente enquanto lei constitutiva baseada na taxa de Truesdell. É obtida dispondo os componentes do tensor de quarta ordem \(G_{ijkl} = \delta_{il}\sigma_{kj} + \delta_{jl}\sigma_{ik}\) na forma de Voigt de \(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} \]

Esta é a matriz resultante.

Montagem da matriz de rigidez global

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

\[ \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 dispostos numa 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 adicionados diretamente no ciclo dos elementos. A matriz é quadrada, com dimensão igual ao número de graus de liberdade por nó \(\times\) número total de nós \(n_g\), mas, como as componentes que não correspondem a pares de nós ligados por elementos são \(0\), é armazenada na forma de matriz esparsa.

As matrizes de rigidez dos elementos para as formulações TL e UL têm a mesma forma, exceto pela mudança da configuração de referência (coordenadas nodais e origem utilizada 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 numa sub-rotina comum.

Algoritmo de iteração

Resumindo o acima exposto, 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)\). Depois, na \(i\)-ésima iteração, executa-se o procedimento seguinte.

  1. Para o 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}}\) utilizando o procedimento de Construção da matriz de rigidez tangente.
  2. Para impor condições de fronteira geométricas, modifique a matriz de rigidez tangente e o vetor residual para os graus de liberdade sujeitos a restrições de deslocamento, obtendo \(\tilde{\boldsymbol{K}}_i, \tilde{\boldsymbol{R}}_{i-1}\) (Tratamento de condições de fronteira 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\). Este procedimento representa frequentemente 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, em conformidade, 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 termine a iteração se esta tiver sido atingida. As componentes correspondentes às reações de restrição aparecem no residual \(\boldsymbol{R}_i\) nos graus de liberdade sujeitos a condições de fronteira geométricas, pelo que o indicador de convergência é construído a partir de \(\tilde{\boldsymbol{R}}_i\) depois de excluir essas componentes. Os indicadores e limiares de convergência específicos são descritos em Critérios de convergência. Se a convergência não for atingida e for alcançado o limite de iterações, a iteração é considerada falhada.

Quando a iteração converge, adicione 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 avance para o passo de tempo seguinte.

Tópicos relacionados