Flambagem elástica.md

Uso de autovetores de flambagem como imperfeição inicial

Autor: Prof. Dr. Eduardo Vicente Wolf Trentini

Ideia central

Em uma análise não linear geométrica, uma estrutura perfeitamente reta e carregada de forma perfeitamente simétrica pode seguir uma trajetória idealizada. Para representar imperfeições geométricas iniciais, pode-se perturbar a geometria inicial usando uma forma modal de flambagem:

$$\mathbf{x}_0^{\mathrm{imp}} = \mathbf{x}_0 + e_0\,\boldsymbol{\phi}_n .$$

O vetor φ_n é uma versão normalizada do modo de flambagem, e e₀ é a amplitude física escolhida para a imperfeição.

Sumário

• Objetivo
• Por que não basta usar o autovetor de K?
• Análise linear de referência
• Matriz geométrica do elemento de pórtico 2D
• Problema de flambagem linear
• Redução para graus de liberdade livres
• Conversão para um problema padrão
• Normalização do modo
• Alteração da geometria inicial
• Algoritmo completo
• Esboço de implementação com Eigen
• Como escolher o modo
• Cuidados numéricos e físicos
• Exemplo conceitual: coluna comprimida
• Resumo

Objetivo

O objetivo deste texto é explicar como obter uma forma de imperfeição inicial para uma estrutura a partir de um problema de autovalor. A sequência usual é:

1. montar a matriz de rigidez elástica K;
2. aplicar uma carga de referência F_ref;
3. resolver uma análise linear para obter u_ref;
4. calcular os esforços normais de referência N_ref nos elementos;
5. montar a matriz de rigidez geométrica K_G;
6. resolver o problema de flambagem linear;
7. usar o primeiro modo de flambagem como imperfeição inicial.

A formulação apresentada é direcionada a estruturas reticuladas planas, especialmente elementos de pórtico 2D de Euler-Bernoulli. A mesma lógica geral aparece em outras formulações, mas a matriz geométrica do elemento muda.

Por que não basta usar o autovetor de K?

Se for resolvido o problema

$$\mathbf{K}\boldsymbol{\phi}_i = \lambda_i \boldsymbol{\phi}_i,$$

obtêm-se direções principais de rigidez da matriz elástica. O autovetor associado ao menor autovalor positivo pode indicar uma forma flexível da estrutura, mas não necessariamente representa o modo de instabilidade provocado pelo carregamento axial.

Na flambagem, o efeito importante é a redução da rigidez lateral causada por esforços normais de compressão. Esse efeito não está contido apenas em K. Ele entra pela matriz de rigidez geométrica K_G, montada a partir dos esforços normais existentes na estrutura.

Conclusão prática

Para imperfeição inicial associada à instabilidade, o mais adequado é usar o autovetor do problema de flambagem linear, e não simplesmente um autovetor de K.

Análise linear de referência

Primeiro, monta-se a rigidez elástica global da estrutura:

$$\mathbf{K}\mathbf{u}=\mathbf{F}.$$

Escolhe-se uma carga de referência F_ref. Ela deve ter o mesmo padrão do carregamento que pode causar instabilidade. Por exemplo, se a estrutura será comprimida por cargas verticais nos nós superiores, essas cargas devem aparecer em F_ref. Pode ser inclusive o carregamento que se deseja aplicar à estrutura.

A análise linear de referência é:

$$\mathbf{K}\mathbf{u}_{\mathrm{ref}}=\mathbf{F}_{\mathrm{ref}}.$$

A partir do vetor u_ref, calcula-se o esforço normal de referência em cada elemento. Para um elemento de pórtico 2D, no sistema local, uma forma simples é:

$$N_{\mathrm{ref}}^{(e)} = \frac{EA}{L}\left(u_{2,\ell}^{(e)} - u_{1,\ell}^{(e)}\right),$$

em que:

• E é o módulo de elasticidade;
• A é a área da seção;
• L é o comprimento do elemento;
• u_1,ℓ^(e) e u_2,ℓ^(e) são os deslocamentos axiais locais dos nós do elemento.

Neste texto será adotada a convenção:

$$N > 0 \quad \text{tração}, \qquad N < 0 \quad \text{compressão}.$$

Com essa convenção, elementos comprimidos têm N < 0 e reduzem a rigidez lateral da estrutura.

Matriz geométrica do elemento de pórtico 2D

Considere um elemento de pórtico plano de Euler-Bernoulli com graus de liberdade locais ordenados como:

$$\mathbf{d}_{\ell}^{(e)} = \begin{bmatrix} u_1 & v_1 & \theta_1 & u_2 & v_2 & \theta_2 \end{bmatrix}^{T}.$$

Uma forma usual da matriz geométrica local é:

$$\mathbf{k}_{G,\ell}^{(e)} = \frac{N^{(e)}}{30L} \begin{bmatrix} 0 & 0 & 0 & 0 & 0 & 0 \\\ 0 & 36 & 3L & 0 & -36 & 3L \\\ 0 & 3L & 4L^2 & 0 & -3L & -L^2 \\\ 0 & 0 & 0 & 0 & 0 & 0 \\\ 0 & -36 & -3L & 0 & 36 & -3L \\\ 0 & 3L & -L^2 & 0 & -3L & 4L^2 \end{bmatrix}.$$

Como a matriz acima está no sistema local do elemento, ela deve ser transformada para o sistema global:

$$\mathbf{k}_{G}^{(e)} = \mathbf{T}_e^{T}\mathbf{k}_{G,\ell}^{(e)}\mathbf{T}_e.$$

Depois, cada matriz k_G^(e) é assemblada na matriz global K_G, do mesmo modo que se faz com a matriz de rigidez elástica.

Sinal da matriz geométrica

A matriz acima foi escrita com N > 0 em tração e N < 0 em compressão. Assim, para elementos comprimidos, k_G^(e) contribui negativamente para a rigidez tangente lateral. Se outra convenção for usada, o sinal do problema de autovalor deve ser ajustado.

Problema de flambagem linear

A rigidez tangente linearizada em torno do estado de referência pode ser escrita como:

$$\mathbf{K}_T(\lambda) = \mathbf{K} + \lambda \mathbf{K}_G,$$

em que λ é um fator multiplicador da carga de referência.

A perda de estabilidade linear ocorre quando existe um vetor não nulo φ tal que:

$$\left(\mathbf{K}+\lambda\mathbf{K}_G\right)\boldsymbol{\phi}=\mathbf{0}.$$

Essa equação pode ser reescrita como o problema de autovalor generalizado:

$$\mathbf{K}\boldsymbol{\phi} = -\lambda\mathbf{K}_G\boldsymbol{\phi}.$$

O menor autovalor positivo λ₁ é o fator crítico de flambagem linear. Se F_ref for a carga de referência, a carga crítica aproximada é:

$$\mathbf{F}_{\mathrm{cr}} \approx \lambda_1\mathbf{F}_{\mathrm{ref}}.$$

O autovetor associado φ₁ é o primeiro modo de flambagem.

Redução para graus de liberdade livres

Para resolver o problema de autovalor, é recomendável trabalhar apenas com os graus de liberdade livres. Separando os graus de liberdade em livres f e restringidos r, escreve-se:

$$\mathbf{u}=\begin{bmatrix}\mathbf{u}_f\\\mathbf{u}_r\end{bmatrix}.$$

Como os apoios impõem u_r = 0, o problema de flambagem deve ser resolvido com as submatrizes livres:

$$\mathbf{K}_{ff}\boldsymbol{\phi}_f = -\lambda \mathbf{K}_{Gff}\boldsymbol{\phi}_f.$$

Depois de calculado φ_f, reconstrói-se o vetor global completo φ, preenchendo com zero os graus de liberdade restringidos.

Cuidado de implementação

Para análise estática linear, é comum impor apoios zerando linhas e colunas da matriz global. Para problema de autovalor, o procedimento mais limpo é montar diretamente as matrizes reduzidas K_ff e K_Gff, contendo apenas os graus de liberdade livres.

Conversão para um problema padrão

Algumas bibliotecas resolvem diretamente o problema generalizado:

$$\mathbf{A}\boldsymbol{\phi}=\lambda\mathbf{B}\boldsymbol{\phi}.$$

Nesse caso, pode-se usar:

$$\mathbf{A}=\mathbf{K}_{ff}, \qquad \mathbf{B}=-\mathbf{K}_{Gff}.$$

Se a biblioteca disponível só resolver problemas padrão do tipo

$$\mathbf{C}\boldsymbol{\phi}=\mu\boldsymbol{\phi},$$

então pode-se multiplicar a equação anterior por K_ff⁻¹:

$$\mathbf{K}_{ff}^{-1}\left(-\mathbf{K}_{Gff}\right)\boldsymbol{\phi}_f = \frac{1}{\lambda}\boldsymbol{\phi}_f.$$

Definindo:

$$\mathbf{C}=\mathbf{K}_{ff}^{-1}\left(-\mathbf{K}_{Gff}\right),$$

tem-se:

$$\mathbf{C}\boldsymbol{\phi}_f = \mu\boldsymbol{\phi}_f, \qquad \mu=\frac{1}{\lambda}.$$

Logo:

$$\lambda=\frac{1}{\mu}.$$

Nesse formato, o menor λ positivo corresponde ao maior μ positivo.

Normalização do modo

O autovetor não tem escala física. Se φ é autovetor, então 10φ, 0,01φ ou −φ também são autovetores.

Para usar o autovetor como imperfeição geométrica, deve-se escolher uma normalização. Para estruturas planas, uma opção simples é normalizar pelo maior deslocamento translacional:

$$\boldsymbol{\phi}_n = \frac{\boldsymbol{\phi}} {\displaystyle \max_i \sqrt{\phi_{x_i}^{2}+\phi_{y_i}^{2}}}.$$

Assim, o maior deslocamento translacional nodal de φ_n passa a valer 1.

A amplitude física da imperfeição pode ser definida por um valor e₀. Por exemplo:

$$e_0=\frac{L_{\mathrm{ref}}}{1000},$$

em que L_ref pode ser o comprimento de uma barra, a altura de um pilar, o vão de uma estrutura, ou outra dimensão de referência adequada ao problema.

Alteração da geometria inicial

Depois da normalização, a geometria imperfeita é obtida por:

$$x_i^{\mathrm{imp}} = x_i + e_0\phi_{x_i,n},$$ $$y_i^{\mathrm{imp}} = y_i + e_0\phi_{y_i,n}.$$

Em um pórtico plano, o autovetor costuma conter, para cada nó:

$$\boldsymbol{\phi}_i = \begin{bmatrix} \phi_{x_i} & \phi_{y_i} & \phi_{\theta_i} \end{bmatrix}^{T}.$$

Para modificar as coordenadas dos nós, usa-se normalmente apenas:

$$\phi_{x_i}, \qquad \phi_{y_i}.$$

A componente rotacional φ_θi pertence ao modo, mas não é uma coordenada geométrica nodal. Portanto, ela não precisa ser aplicada diretamente como alteração de geometria.

Algoritmo completo

A sequência de implementação pode ser resumida assim:

1. Montar K elástica global.

2. Montar o vetor de carga de referência F_ref.

3. Identificar graus de liberdade livres.

4. Resolver K_ffu_ref,f = F_ref,f.

5. Reconstruir u_ref global.

6. Para cada elemento:

   1. transformar deslocamentos globais para locais;
   2. calcular N_ref do elemento;
   3. montar k_G,ℓ^(e) local com N_ref;
   4. transformar k_G,ℓ^(e) para coordenadas globais;
   5. assemblar k_G^(e) na matriz global K_G.

7. Reduzir K e K_G para os graus de liberdade livres.

8. Resolver o problema de flambagem linear:

   $$\mathbf{K}_{ff}\boldsymbol{\phi} = -\lambda\mathbf{K}_{Gff}\boldsymbol{\phi}.$$

9. Escolher o menor λ positivo. Alternativamente, permitir a visualização dos modos e deixar o usuário escolher.

10. Reconstruir φ global.

11. Normalizar φ pelo maior deslocamento translacional.

12. Escolher e₀.

13. Criar geometria imperfeita: x^imp = x + e₀φ_x; y^imp = y + e₀φ_y.

14. Usar essa geometria como ponto inicial da análise não linear geométrica.

Esboço de implementação com Eigen

O trecho abaixo é apenas um esqueleto. Ele mostra a lógica, não uma implementação completa.

// Matrizes globais completas
Eigen::MatrixXd K  = AssembleElasticStiffness(model);
Eigen::VectorXd Fr = AssembleReferenceLoad(model);
Eigen::MatrixXd KG = Eigen::MatrixXd::Zero(K.rows(), K.cols());

// Lista de graus de liberdade livres
std::vector<int> freeDofs = BuildFreeDofList(model);

// Redução para graus de liberdade livres
Eigen::MatrixXd Kff  = ExtractSubmatrix(K, freeDofs, freeDofs);
Eigen::VectorXd Frf  = ExtractSubvector(Fr, freeDofs);

// Análise linear de referência
Eigen::VectorXd ur_f = Kff.ldlt().solve(Frf);
Eigen::VectorXd ur   = ExpandToGlobalVector(ur_f, freeDofs, K.rows());

// Cálculo dos esforços normais e montagem de KG
for (Element& e : model.elements)
{
    Eigen::VectorXd ue_global = ExtractElementDisplacements(ur, e);
    Eigen::VectorXd ue_local  = e.T * ue_global;

    double N = e.E * e.A / e.L * (ue_local(3) - ue_local(0));

    Eigen::Matrix<double, 6, 6> kG_local =
        GeometricStiffness2DFrame(e.L, N);

    Eigen::Matrix<double, 6, 6> kG_global =
        e.T.transpose() * kG_local * e.T;

    AssembleElementMatrix(KG, kG_global, e.globalDofs);
}

Eigen::MatrixXd KGff = ExtractSubmatrix(KG, freeDofs, freeDofs);

// Problema padrão equivalente: C * phi = mu * phi
// C = Kff^{-1} * (-KGff), mu = 1/lambda
Eigen::MatrixXd C = Kff.ldlt().solve(-KGff);
Eigen::EigenSolver<Eigen::MatrixXd> solver(C);

// Selecionar maior mu positivo real; lambda = 1/mu.
// O autovetor associado e o modo de flambagem.

Para uma implementação mais robusta, é preferível usar um resolvedor de autovalor generalizado simétrico quando disponível, evitando formar explicitamente K⁻¹K_G. Em muitos casos educacionais e protótipos pequenos, a forma padrão acima é suficiente para estudo, desde que os resultados sejam verificados.

Como escolher o modo

Depois de resolver o problema de autovalor, deve-se filtrar os resultados:

1. descartar autovalores complexos com parte imaginária significativa;
2. descartar autovalores negativos ou muito próximos de zero, conforme a convenção adotada;
3. ordenar os autovalores positivos;
4. escolher o menor λ positivo;
5. verificar visualmente se o autovetor corresponde a uma forma de flambagem plausível.

Em estruturas simétricas, pode haver modos próximos, ou até modos repetidos. Nesses casos, pequenas diferenças numéricas podem alterar a combinação linear retornada pelo resolvedor. Isso não significa necessariamente erro.

Cuidados numéricos e físicos

Apoios insuficientes

Se a estrutura tiver modos de corpo rígido, K será singular ou quase singular. Antes de fazer análise de flambagem, verifique se a análise linear está bem condicionada.

Sinal dos esforços normais

A matriz geométrica depende diretamente de N. Se a convenção de sinal do elemento for alterada, o problema de autovalor também muda. Sempre teste com um caso conhecido, como uma coluna biarticulada.

Escala da imperfeição

O autovetor só fornece a forma. A amplitude e₀ deve ser escolhida separadamente. Para estudo numérico, valores como L/1000, L/500 ou L/300 são comuns como hipóteses de imperfeição. Para projeto, deve-se seguir a norma ou critério adotado.

Não aplicar rotações como deslocamento

Em pórticos 2D, o autovetor inclui rotações nodais. Para alterar a geometria, use apenas os deslocamentos translacionais. As rotações podem ser usadas para visualização do modo, mas não devem ser somadas às coordenadas x e y. Existem funções de visualização que usam as rotações para desenhar a forma deslocada entre os nós. Em uma implementação inicial, isso pode ser omitido.

Magnitude da carga de referência

Se F_ref for multiplicada por um fator c, os esforços normais de referência também serão multiplicados por c, e o autovalor crítico será aproximadamente dividido por c. O produto λF_ref, porém, representa a carga crítica aproximada. Por isso, a forma do carregamento é mais importante do que sua magnitude absoluta.

Exemplo conceitual: coluna comprimida

Considere uma coluna reta com comprimento L, comprimida por uma carga axial de referência P_ref. O primeiro modo de flambagem de uma coluna biarticulada tem forma semelhante a:

$$v_0(x)=e_0\sin\left(\frac{\pi x}{L}\right).$$

No modelo discreto, o autovetor de flambagem fornece uma aproximação nodal dessa forma. A geometria imperfeita é criada deslocando os nós transversalmente:

$$y_i^{\mathrm{imp}} = y_i + e_0\phi_{y_i,n}.$$

Essa imperfeição faz com que a análise não linear já comece com uma pequena curvatura inicial, permitindo capturar os efeitos de segunda ordem associados ao carregamento compressivo.

Resumo

Etapa	O que fazer
Análise linear	Resolver Kuref = Fref.
Esforços normais	Calcular Nref(e) em cada elemento.
Matriz geométrica	Montar KG usando os esforços normais de referência.
Flambagem linear	Resolver Kφ = −λKGφ.
Modo crítico	Escolher o menor λ positivo.
Normalização	Escalar φ para que o maior deslocamento translacional seja 1.
Imperfeição inicial	Aplicar x0imp = x0 + e0φn.
Análise não linear	Usar a geometria imperfeita como geometria inicial.

Mensagem final

O autovetor fornece a forma da imperfeição inicial, mas não define sua amplitude. A matriz geométrica K_G representa a influência dos esforços normais sobre a rigidez tangente e, consequentemente, sobre a estabilidade da estrutura. Portanto, a sequência de análise consiste em: realizar uma análise linear de referência, calcular os esforços normais nos elementos, montar a matriz geométrica K_G, resolver o problema de flambagem e aplicar o modo normalizado como imperfeição geométrica inicial.

Ressalta-se que, para que o modo de flambagem seja representado adequadamente, as barras devem estar discretizadas em segmentos menores. Caso contrário, uma barra longa modelada por um único elemento retilíneo pode não representar satisfatoriamente a forma deformada associada ao modo crítico.