Uma Abordagem Analítica e Numérica para Detecção de Pontos Limites e de Bifurcação

RESUMO

Neste trabalho descreve-se de maneira detalhada, tanto analíticamente quanto numericamente, a detecção e a classificação de pontos críticos na trajetória primária de equilíbrio de sistemas estruturais. Utiliza-se a Formulação Lagrangiana Total para descrever a cinemática de um elemento de barra bi-articulado 3D. Através desta formulação obtém-se o vetor de forças internas e a matriz de rigidez tangente que levam em conta os efeitos da não linearidade geométrica. Assume-se um modelo constitutivo linear elástico para o estado uniaxial de tensão-deformação, usando a deformação de Green-Lagrange e a tensão axial do segundo tensor de Piola-Kirchhoff que são energéticamente conjugados. Como estudo de caso apresenta-se um sistema físico simples com três graus de liberdades composto por duas barras bi-articuladas 3D e uma mola linear. Por fim, determinam-se as condições geométricas e físicas para a coalescência entre os pontos limites e de bifurcação.

Palavras-chave:
Descrição Lagrangiana total; análise não-linear geométrica; pontos críticos

ABSTRACT

Using both analytical and numerical approaches, this paper describes in detail the detection and classification of critical points in the primary equilibrium path of structural systems. The Total Lagrangian formulation is employed to describe the kinematics of a biarticulated prismatic 3D bar element. With this formulation, the internal force vector and the tangent stiffness matrix including the geometric nonlinearity effects are obtained. An elastic linear constitutive model is assumed for the uniaxial stress-strain state. Such model uses the Green-Lagrange strain tensor and the second Piola-Kirchhoff axial stress tensor which are energetically conjugate tensors. As a study case, the article presents a simple structural system with three degrees of freedom made up of two bi-articulated prismatic 3D bars and a linear spring. Finally, the geometrical and physical conditions for the coalescence between limit and bifurcation points are determined.

Keywords:
Total Lagrangian description; geometrical nonlinearity; critical points

1 INTRODUÇÃO

Tornou-se quase imperativo, para o avanço da ciência, compreender e simular os fenômenos não-lineares em diversas áreas do conhecimento, tais como, biomecânica, mecânica dos fluídos, geotecnia, mecânica dos sólidos, engenharia de tecidos humanos, etc. Por exemplo nas industrias aeronáutica, aeroespacial e petrolífera a análise não linear é imprescindível no projeto de diferentes tipologias estruturais aplicadas nesses setores. Por outro lado, para modelar esses fenômenos de maneira consistente é necessário o conhecimento de física, matemática aplicada e computacional. Nas últimas décadas muitos autores têm publicados livros texto abordando diferentes tópicos da análise não linear na área dos métodos numéricos aplicados a engenharia, por exemplo, recomenda-se a leitura das seguintes referências: Belytschko et al. ¹, Bonet e Wood ², Borst et al. ³, Crisfield ^4),(5, Doyle ⁶, Hashiguchi e Yamakawa ⁷, Kojic´ e Bathe ⁸, Krenk ⁹, Neto et al. ¹⁰, Simo e Hughes ¹¹, Voyiadjis e Woelke ¹² e Wriggers ¹³. Decorre, dáı, a necessidade de difundir os conceitos básicos e fundamentais da análise não-linear através do método dos elementos finitos, imbuido deste espírito, o trabalho, aqui apresentado, utilizando um sistema físico com apenas três graus de liberdade, tem como objetivo apresentar uma análise teórica e numérica da instabilidade estrutural que consiste na detecção e classificação de pontos singulares na trajetória de equilíbrio primária. Propõe-se utilizar elementos de barra bi-articulados por sua simplicidade teórica, o que permite descrever, facilmente, a cinemática do movimento do elemento e obter o vetor de forças internas e a matriz de rigidez tangente analíticamente, que são elementos imprescindíveis para uma análise não-linear em mecânica dos sólidos. Para simular deformações elásticas finitas, assume-se um modelo constitutivo hiperelástico para o estado uniaxial de tensão-deformação, utilizando a tensão axial do segundo tensor de Piola-Kirchhoff e a deformação de Green-Lagrange, que formam um par energeticamente conjugado. Nos itens 2, 3 e 4 descrevem-se a cinemática do elemento bi-articulado 3D adotando a Formulação Lagrangiana Total (FLT), a obtenção do vetor de forças internas e a dedução da matriz de rigidez tangente, respectivamente. No item 5 faz-se uma abordagem analítica detalhada da detecção e classificação dos pontos singulares presentes na trajeto´ ria primária de equilíbrio utilizando um sistema físico bastante simples composto por duas barras bi-articuladas 3D e uma mola linear. No item 6 faz-se uma abordagem numérica da instabilidade estrutural do sistema físico descrito anteriormente, para isso utiliza-se um programa de elementos finitos que realiza análises não lineares geométricas escrito em Fortran90 pelos autores deste artigo. Por último, apresentam-se as conclusões finais e as referências bibliográficas utilizadas neste trabalho.

2 DESCRIÇÃO CINEMÁTICA

Seja um sistema de coordenadas globais cuja base é ortonormal conforme mostra-se na Figura 1. Para expressar as variáveis cinemáticas na configuração indeformada utilizam-se as coordenadas materiais (X , Y, Z ), enquanto que, na configuração deformada as coordenadas espaciais (x, y, z). Na configuração indeformada as coordenadas nodais do elemento bi-articulado 3D são dadas por X _A = (X _A , Y _A , Z _A ) e X _B = (X _B , Y _B , Z _B ), respectivamente. Sua posição e comprimento iniciais são dados por

$$\begin{gathered} \overrightarrow{AB}={\mathit{X}}_0={\mathit{X}}_B-{\mathit{X}}_A \\ {l}_0^2={\mathit{X}}_{\mathbf{0}}^{\mathbf{T}}{\mathit{X}}_{\mathbf{0}} \end{gathered}$$ (2.1)

Figura 1:
Movimento do elemento bi-articulado 3D.

Enquanto que na configuração deformada as coordenadas nodais do elemento bi-articulado 3D são dadas por x _A = (x _A , y _A , z _A ) e x _B = (x _B , y _B , z _B ), respectivamente. Conforme mostra-se na Figura 1 a coordenada atual do nó A é dada por x _A = X _A + u _A , enquanto que a coordenada atual do nó B se expressa como x _B = X _B + u _B . Onde u _A = (u _A , v _A , w _A ) é o deslocamento do nó A e u _B = (u _B , v _B , w _B ) é o deslocamento do nó B. Portanto, a posição atual do elemento é dada por

$$\begin{gathered} \overrightarrow{ab}=\mathit{x}={\mathit{x}}_B-{\mathit{x}}_A \\ \overrightarrow{ab}=\mathit{x}=({\mathit{X}}_B+{\mathit{u}}_B)-({\mathit{X}}_A-{\mathit{u}}_A) \\ \overrightarrow{ab}=\mathit{x}=({\mathit{X}}_B-{\mathit{X}}_A)+({\mathit{u}}_B-{\mathit{u}}_A) \\ \overrightarrow{ab}=\mathit{x}={\mathit{X}}_0+{\mathit{u}}_{BA} \end{gathered}$$ (2.2)

Onde u _BA = u _B − u _A é o vetor de deslocamentos nodais relativos conforme mostra-se na Figura 2. Por outro lado, o comprimento atual se expressa como

$l^2={\mathit{x}}^T\mathit{x}=({\mathit{X}}_0+{\mathit{u}}_{BA}{)}^T({\mathit{X}}_0+{\mathit{u}}_{BA})$ (2.3)

Figura 2:
Vetor de deslocamentos nodais relativos.

Neste trabalho adota-se a formulação Lagrangiana Total para descrever o movimento do elemento bi-articulado 3D, portanto serão usadas as coordenadas materiais (X , Y, Z ) e a configuração indeformada para definir a medida de deformação do elemento. Dentre algumas familias de deformação descritas em coordenadas materiais existentes na literatura técnica, adota-se neste trabalho, a medida de deformação de Green-Lagrange que compara os quadrados dos comprimentos atual l e inicial l ₀ do elemento da seguinte maneira

${\epsilon }_G=\frac{l^2-l_0^2}{2l_0^2}=\frac{1}{l_0^2}({\mathit{X}}_0^T{\mathit{u}}_{BA}+\frac{1}{2}{\mathit{u}}_{BA}^T{\mathit{u}}_{BA})$ (2.4)

Note-se que esta medida de deformação possui termos quadráticos em relação aos deslocamentos nodais relativos.

Para obter o vetor de forças internas utiliza-se o Princípio do Trabalhos Virtuais (PTV), portanto é necessário aplicar uma variação virtual no campo de deslocamentos na configuração de equilíbrio atual, conforme mostra-se na Figura 3, o que implica em uma variação virtual da deformação de Green-Lagrange que se escreve como

$$\begin{gathered} \delta {\epsilon }_G=\frac{1}{l_0^2}({\mathit{X}}_0^T\delta {\mathit{u}}_{BA}+\frac{1}{2}\delta {\mathit{u}}_{BA}^T{\mathit{u}}_{BA}+\frac{1}{2}{\mathit{u}}_{BA}^T\delta {\mathit{u}}_{BA}) \\ \delta {\epsilon }_G=\frac{1}{l_0^2}({\mathit{X}}_0^T\delta {\mathit{u}}_{BA}+{\mathit{u}}_{BA}^T\delta {\mathit{u}}_{BA}) \\ \delta {\epsilon }_G=\frac{1}{l_0^2}({\mathit{X}}_0\mathbf{}+{\mathit{u}}_{BA}{)}^T\delta {\mathit{u}}_{BA}=\frac{1}{l_0^2}{\mathit{x}}^T\delta {\mathit{u}}_{BA} \end{gathered}$$ (2.5)

onde δu _BA = δu _B − δu _A é a variação virtual do vetor de deslocamentos relativos. Note-se que a variação virtual da deformação de Green-Lagrange consiste na projeção da variação virtual do vetor dos deslocamentos relativos sobre a posição atual do elemento definida pelo vetor x e escalado por $l_0^2$.

Figura 3:
Variação virtual dos deslocamentos nodais.

3 VETOR DE FORÇAS NODAIS

Como mostra-se na Figura 4, seja f _A =( f _Ax , f _Ay , f _Az ) o vetor de forças do nó A, e f _B = ( f _Bx , f _By , f _Bz ) o vetor de forças do nó B, respectivamente. Para obter estes vetores de forças nodais aplica-se o PTV na configuração indeformada pois se está utilizando a medida de deformação de Green-Lagrange, desta maneira expressa-se este princípio como

$\delta V=\underset{0}{\overset{l_0}{\int }}N\delta {\epsilon }_{G}ds-{\mathit{f}}_A^T\delta {\mathit{u}}_A-{\mathit{f}}_B^T\delta {\mathit{u}}_B=0,$ (3.1)

Figura 4:
Vetor de forças nodais e a variação virtual dos deslocamentos nodais.

onde N é o esforço axial que atua no elemento e é dado por N = σ _G A. Lembrando que a tensão axial σ _G é energéticamente conjugada com a medida de deformação de Green-Lagrange e é uma das tensões normais do segundo tensor de tensões de Piola-Kirchoff. Substituindo a equação (2.5c) na equação (3.1), tem-se que

$$\begin{gathered} \delta V=\underset{0}{\overset{l_0}{\int }}\frac{N}{l_0^2}\mathit{x}{}^T\delta {\mathit{u}}_{BA}ds-{\mathit{f}}_A^T\delta {\mathit{u}}_A-{\mathit{f}}_B^T\delta {\mathit{u}}_B=0 \\ \delta V=\underset{0}{\overset{l_0}{\int }}\frac{N}{l_0^2}\mathit{x}{}^T(\delta {\mathit{u}}_B-\delta {\mathit{u}}_{\mathbf{A}}\mathbf{)}ds-{\mathit{f}}_A^T\delta {\mathit{u}}_A-{\mathit{f}}_B^T\delta {\mathit{u}}_B=0 \\ \delta V=\delta {\mathit{u}}_A^T\mathbf{}\left(\underset{0}{\overset{l_0}{\int }}\frac{N}{l_0^2}\mathit{x}ds-{\mathit{f}}_A\right)+\delta {\mathit{u}}_B^T\mathbf{}\left(\underset{0}{\overset{l_0}{\int }}\frac{N}{l_0^2}\mathit{x}ds-{\mathit{f}}_B\right)=0 \\ {\mathit{f}}_A=\frac{N}{l_0}\mathit{x},{\mathit{f}}_B=\frac{N}{l_0}\mathit{x} \end{gathered}$$ (3.2)

Neste trabalho assume-se que σ _G = Eε _G , onde E é o módulo de elasticidade longitudinal do material, assim o esforço axial pode ser definido como N = E Aε _G e a equação (3.2d) pode ser reescrita da seguinte maneira

${\mathit{f}}_A=\frac{EA}{l_0}{\epsilon }_G\mathit{x},{\mathit{f}}_B=\frac{EA}{l_0}{\epsilon }_G\mathit{x}$ (3.3)

4 MATRIZ DE RIGIDEZ TANGENTE

Ao aplicar um incremento infinitesimal nos vetores de deslocamentos nodais u _A e u _B na configuração deformada, obtém-se um incremento infinitesimal dos vetores de forças internas f _A e f _B respectivamente, através da matriz de rigidez tangente. Portanto, define-se a relação entre os incrementos infinitesimais dos vetores de forças internas e dos vetores de deslocamentos nodais como

$\left\{\begin{array}{c}d{f}_A\\ d{f}_B\end{array}\right\}=K_T\left\{\begin{array}{c}d{u}_A\\ d{u}_B\end{array}\right\}.$ (4.1)

Onde K _T é a matriz de rigidez tangente de ordem 6×6. Portanto, levando em conta a equação (2.2d) e diferenciando-se a equação (3.2d) em relação ao vetor de deslocamentos relativos, obtém-se que

$$\begin{gathered} d{\mathit{q}}_A=-\mathit{x}\frac{dN}{l_0}-\frac{N}{l_0}d\mathit{x} \\ d{\mathit{q}}_A=-\left(\frac{\mathbf{x}}{{\mathbf{l}}_{\mathbf{0}}}\frac{\mathbf{d}\mathbf{N}}{\mathbf{d}{\mathbf{u}}_{\mathbf{B}\mathbf{A}}}\mathbf{+}\frac{\mathbf{N}}{{\mathbf{l}}_{\mathbf{0}}}\mathbf{I}\right)d({\mathit{u}}_B-{\mathit{u}}_A) \\ d{\mathit{q}}_B=-d{\mathit{q}}_A, \end{gathered}$$ (4.2)

onde I é a matriz identidade de ordem 3×3. Por outro lado, diferenciando-se o esforço axial em relação ao vetor dos deslocamentos relativos, e levando em conta as equações (2.2d), (2.4) e (2.5c), chega-se a

$\frac{dN}{d{\mathit{u}}_{BA}}=EA\frac{d{\epsilon }_G}{d{\mathit{u}}_{BA}}=\frac{EA}{l_0^2}({\mathit{X}}_0^T+{\mathit{u}}_{BA}^T)=\frac{EA}{l_0^2}{\mathit{x}}^T$ (4.3)

Por último, substituindo a equação (4.3) na equação (4.2b) e levando em conta a equação (4.3c), obtém-se que

$\left\{\frac{d{\mathit{q}}_A}{d{\mathit{q}}_{\mathbf{B}}}\right\}=\left(\frac{EA}{l_0^3}\left[\begin{array}{cc}\mathit{x}\otimes \mathit{x}& -\mathit{x}\otimes \mathit{x}\\ -\mathit{x}\otimes \mathit{x}& \mathit{x}\otimes \mathit{x}\end{array}\right]+\frac{N}{l_0}\left[\begin{array}{cc}\mathit{I}& -\mathit{I}\\ -\mathit{I}& \mathit{I}\end{array}\right]\right)\left\{\frac{d{\mathit{u}}_A}{d{\mathit{u}}_{\mathbf{B}}}\right\}$ (4.4)

com

$$\begin{gathered} {\mathit{K}}_M=\frac{EA}{l_0^3}\left[\begin{array}{cc}\mathit{x}\otimes \mathit{x}& -\mathit{x}\otimes \mathit{x}\\ -\mathit{x}\otimes \mathit{x}& \mathit{x}\otimes \mathit{x}\end{array}\right],{\mathit{K}}_{\sigma }=\frac{N}{l_0}\left[\begin{array}{cc}\mathit{I}& -\mathit{I}\\ -\mathit{I}& \mathit{I}\end{array}\right] \\ \mathit{x}\otimes \mathit{x}=\left[\begin{array}{ccc}{x}_{BA}^2& x_{BA}{y}_{BA}& x_{BA}{z}_{BA}\\ x_{BA}{y}_{BA}& y_{BA}^2& y_{BA}{z}_{BA}\\ x_{BA}{z}_{BA}& y_{BA}{z}_{BA}& z_{BA}^2\end{array}\right], \end{gathered}$$ (4.5)

onde K _M é a matriz de rigidez material, de ordem 6×6, que depende do vetor posição atual do elemento x cujas componentes são: x _BA = x _B − x _A , y _BA = y _B − y _A e z _BA = z _B − z _A . K _σ é a matriz de rigidez geométrica que depende do esforço axial N e é de ordem 6×6 e o símbolo ⊗é o produto tensorial ou aberto. Portanto, a matriz de rigidez tangente se expressa como K _T = K _M + K _σ .

5 FORMULAÇÃO ANALÍTICA

Para detectar os pontos críticos na trajetória primaria de equilíbrio analíticamente adota-se um sistema estrutural composto por 2 elementos bi-articulados 3D. Na Figura 5 mostram-se as condições de contorno, de carregamento, as propriedades mecânicas e geométricas deste sistema estrutural. Note-se que no vértice como condição de contorno há uma mola linear de constante k

Figura 5:
Treliça espacial com uma mola lateral.

na direção do eixo z. Neste nó aplicam-se as cargas f ₂ e f ₃ nas direções dos eixos y e z, respectivamente. Assume-se que não haverá deslocamento na direção do eixo x devido a ausência de carga aplicada nesta direção. Além disso, adota-se o comprimento c na direção do eixo z como uma imperfeição do sistema estrutural. De acordo com a Figura 5 o comprimento inicial de cada barra bi-articulada é dado por l ₀ = $\sqrt{a^2+b^2+c^2}$. Após o deslocamento do vértice dado por (0, u ₂ , u ₃) o comprimento atual de cada barra será

$l=\sqrt{(a-u_2{)}^2+b^2+(c+u_3{)}^2}.$

Em primeiro lugar, será feita a análise considerando o sistema estrutural perfeito o que implica na seguinte condição c = 0 e f ₃ = 0. Observe-se que isto não implica necessariamente que u ₃ será igual a zero como se demonstrará mais adiante. Portanto, para o sistema perfeito a trajetória primária de equilíbrio, a curva f ₂ ×u ₂ , ocorrerá no plano (x, y).

Para a facilidade do desenvolvimento algébrico a seguir adotam-se os seguintes parâmetros adimensionais

$$\begin{gathered} \alpha =\frac{a}{l_0},\beta =\frac{b}{l_0},\gamma =\frac{c}{l_0} \\ \kappa =\frac{k{l}_0}{EA},{\mu }_2=\frac{u_2}{l_0},{\mu }_3=\frac{u_3}{l_0} \\ {\lambda }_2=\frac{f_2}{EA},{\lambda }_3=\frac{f_3}{EA} \end{gathered}$$ (5.1)

Levando em consideração as definições de l0 e l e os parâmetros dados pela equação (5.1), a deformação de Green-Lagrange de cada barra bi-articulada é dado por

$$\begin{gathered} {\epsilon }_G=\frac{l^2-l_0^2}{2l_0^2}=\frac{(a-u_2{)}^2+b^2+(c+u_3{)}^2-(a^2+b^2+c^2)}{2l_0^2} \\ {\epsilon }_G=\frac{-2a{u}_2+u_2^2+2c{u}_3+u_3^2}{2l_0^2}=-{\mu }_2(\alpha -\frac{1}{2}{\mu }_2)+{\mu }_3(\gamma +\frac{1}{2}{\mu }_3) \end{gathered}$$ (5.2)

A variação virtual da deformação de Green-Lagrange em relação à variação virtual dos deslocamentos (0, δu ₂ , δu ₃ ) escreve-se como

${\delta }_{{\epsilon }_G}=\delta {\mu }_2(\alpha -{\mu }_2)+\delta {\mu }_3(\gamma +{\mu }_3).$ (5.3)

Por outro lado, levando em conta os parâmetros da equação (5.1), a energia potencial total do sistema estrutural é dada pela seguinte expressão

$$\begin{gathered} \pi =2(\frac{1}{2}EA{l}_0{\epsilon }_G^2)+\frac{1}{2}k{u}_3^2-f_2{u}_2-f_3{u}_3 \\ \pi =EA{l}_0({\epsilon }_G^2+\frac{1}{2}\frac{k{u}_3^2}{EA{l}_0}-\frac{f_2{u}_2}{EA{l}_0}-\frac{f_3{u}_3}{EA{l}_0} \\ \pi =EA{l}_0({\epsilon }_G^2+\frac{1}{2}\kappa {\mu }_3^2-{\lambda }_1{\mu }_2-{\lambda }_2{\mu }_3) \end{gathered}$$ (5.4)

Como se está considerando o sistema estrutural com 2 graus de liberdade, sua condição de equilíbrio é representada por um sistema não linear de 2 equações e 2 incógnitas. Para obter este sistema de equações aplica-se o princípio da estacionaridade no funcional da energia potencial total que se escreve como

$\delta \pi =\frac{\partial \pi }{\partial u_2}\delta u_2+\frac{\partial \pi }{\partial u_3}\delta u_3=0$ (5.5)

Levando em consideração os parâmetros definidos na equação (5.1) e aplicando a condição de estacionaridade expressa na equação (5.5) na equação (5.4c), chega-se a

$$\begin{gathered} \delta \pi =EA{l}_0(2_{{\epsilon }_G}{\delta }_{{\epsilon }_G}+\kappa {\mu }_3\delta {\mu }_3-{\lambda }_2\delta {\mu }_2-{\lambda }_3\delta {\mu }_3)=0 \\ {2}_{{\epsilon }_G}=(-\delta {\mu }_2(\alpha -{\mu }_2)+\delta {\mu }_3(\gamma +{\mu }_3\left)\right)+\kappa {\mu }_3\delta {\mu }_3-{\lambda }_2\delta {\mu }_2-{\lambda }_3\delta {\mu }_3=0, \end{gathered}$$ (5.6)

o que resulta no seguinte sistema de equações

$$\begin{gathered} -2_{{\epsilon }_G}(\alpha -{\mu }_2)={\lambda }_2 \\ {2}_{{\epsilon }_G}(\gamma +{\mu }_3)+\kappa {\mu }_3={\lambda }_3. \end{gathered}$$ (5.7)

Esta últimas equações representam as condições de equilíbrio na configuração deformada da treliça espacial nas direções dos eixos y e z, respectivamente. Note-se que ambas equações dependem de termos cúbicos em relação aos parâmetros μ ₂ e μ ₃ pois a deformação ε _G , de acordo com a equação (5.2b), depende desses parâmetros quadráticamente. Para obter a matriz de rigidez tangente diferencia-se os parâmetros de carga em relação aos parâmetros de deslocamentos em ambas equações de equilíbrio tal que

$\left\{\begin{array}{c}d{\lambda }_2\\ d{\lambda }_3\end{array}\right\}=\left[\begin{array}{cc}\frac{\partial {\lambda }_2}{\partial {\mu }_2}& \frac{\partial {\lambda }_2}{\partial {\mu }_3}\\ \frac{\partial {\lambda }_3}{\partial {\mu }_2}& \frac{\partial {\lambda }_3}{\partial {\mu }_3}\end{array}\right]\left\{\begin{array}{c}d{\mu }_2\\ d{\mu }_3\end{array}\right\}$ (5.8)

Determinam-se os coeficientes da matriz de rigidez tangente diferenciando-se as equações (5.2b) e (5.7) de acordo com a equação (5.8). Após alguns desenvolvimentos algébricos chega-se a

$\left\{\begin{array}{c}d{\lambda }_2\\ d{\lambda }_3\end{array}\right\}=\left[\begin{array}{cc}2{\left(\alpha -{\mu }_2\right)}^2+2{\epsilon }_G& -2\left(\alpha -{\mu }_2\right)\left(\gamma +{\mu }_3\right)\\ -2\left(\alpha -{\mu }_2\right)\left(\gamma +{\mu }_3\right)& 2{\left(\gamma -{\mu }_3\right)}^2+2{\epsilon }_G+\kappa \end{array}\right]\left\{\begin{array}{c}d{\lambda }_1\\ d{\lambda }_2\end{array}\right\}$ (5.9)

A matriz de rigidez tangente pode ser decomposta na matriz de rigidez material K _M e na matriz de rigidez geométrica K _σ que se expressam como

${\mathit{K}}_m=\left[\begin{array}{cc}2(\alpha -{\mu }_2{)}^2& -2(\alpha -{\mu }_2)(\gamma +{\mu }_3)\\ -2(\alpha -{\mu }_2)(\gamma +{\mu }_3)& 2(\gamma -{\mu }_3{)}^2+\kappa \end{array}\right],{\mathit{K}}_{\mathbf{\sigma }}=\left[\begin{array}{cc}{2}_{{\epsilon }_G}& 0\\ 0& 2_{{\epsilon }_G}\end{array}\right]$ (5.10)

Agora analisa-se a trajetoria de equilíbrio do sistema perfeito, isto é, com c = 0 ⇒ γ = 0 e com f ₃ = 0 ⇒ λ ₃ = 0, o que implica que a equação de equilíbrio na direção do eixo z deve cumprir que (2ε _G + κ)μ ₃ = 0. Esta equação pode ser satisfeita de duas maneiras: para μ ₃ = 0 ou 2ε _G + κ = 0. Para a primeira condição a trajetória de equilíbrio primária ocorre no plano (x, y). Esta trajetória é simétrica pois somente há o deslocamento vertical do nó do vértice da treliça espacial. De acordo com a equação (5.2b) para esta condição a deformação é igual a ${\epsilon }_G={\mu }_2(\alpha -\frac{1}{2}{\mu }_2)$ Substituindo esta expressão na equação (5.7a), obtém-se a trajetória de equilíbrio primária que se expressa como

${\lambda }_2=2{\alpha }^2{\mu }_2-3\alpha {\mu }_2^2+{\mu }_2^3$ (5.11)

Trata-se de um polinômio de grau 3 em função do parâmetro μ ₂ cujo gráfico mostra-se na Figura 6. Esta curva foi obtida para os seguintes valores: a = 1, b = 2 e c = 0. Como pode-se observar nesta figura, a trajetória de equilíbrio apresenta dois extremos, os pontos A e B, que são obtidos através da condição dλ ₂ / dμ ₂ = 0. Desta maneira, obtêm-se os pontos de máximo e mínimo com os seguintes valores

$$\begin{gathered} {\mu }_2^a=\alpha \left(1-\frac{1}{\sqrt{3}}\right)\Rightarrow {\lambda }_2^a=\frac{2}{3\sqrt{3}}{\alpha }^3 \\ {\mu }_2^b=\alpha \left(1-\frac{1}{\sqrt{3}}\right)\Rightarrow {\lambda }_2^b=-\frac{2}{3\sqrt{3}}{\alpha }^3 \end{gathered}$$ (5.12)

Figura 6:
Trajetória primária de equilíbrio.

Os pontos A e B são denominados na literatura técnica de pontos limites. Nestes pontos a matriz de rigidez tangente definida na equação (5.9) torna-se singular, portanto a estrutura atinge sua capacidade portante máxima no ponto A, tornando-se instável a partir deste ponto. Conforme mostra-se na Figura 6, a estrutura dá um salto adiante do ponto A ao ponto A’ onde a configuração de equilíbrio torna-se estável. Na literatura técnica este fenômeno é denominado de snap-through.

Por outro lado, a segunda condição ocorre quando a deformação atinge o valor de ${\epsilon }_G=-\frac{1}{2}\kappa \text{e}{\mu }_3\ne 0$. Neste caso, inicialmente a treliça espacial encontra-se na trajetória primária de equilíbrio, entretanto quando este equilíbrio torna-se instável, a treliça move-se para fora do

plano (x, y) buscando uma trajetória de equilíbrio secundária estável, desta maneira ocorrendo o fenômeno da flambagem. Para determinar o deslocamento na direção do eixo z utiliza-se a equação (5.2b), que para este caso expressa-se como ${\epsilon }_G=-{\mu }_2(\alpha -\frac{1}{2}{\mu }_2)+\frac{1}{2}{\mu }_3^2$, igualando essa equação à condição ε = $-\frac{1}{2}\kappa$, obtém-se que

${\mu }_3=\pm \sqrt{-\kappa +2{\mu }_2(\alpha -\frac{1}{2}{\mu }_2},$ (5.13)

que é válida para o intervalo: $\alpha -\sqrt{{\alpha }^2-\kappa }\le \alpha +\sqrt{{\alpha }^2-\kappa }$. A relação entre os parâmetros de deslocamentos μ₂ e μ₃ é mostrado na Figura 7. Para obter o gráfico desta figura adotaram-se os seguintes valores para as variáveis a = 1, b = 2, c = 0 e κ = 0.1, consequentemente, 0.131 ≤ μ ₂ ≤ 0.763. Observe-se que na Figura 8 há dois sentidos possíveis, do ponto vista analítico, da treliça flambar fora do plano (x, y) na direção do eixo z. A treliça pode flambar na direção positiva do eixo z, isto e´, para u ₃ > 0, ou na direcção negativa do eixo z com u ₃ < 0. O deslocamento máximo na direção z, pontos E e F, é dada por

$\frac{d{\mu }_3}{d{\mu }_2}=0\Rightarrow {\mu }_3=\pm \sqrt{{\alpha }^2-\kappa }.$

Figura 7:
Curva no espaço de deslocamentos adimensionais.

Figura 8:
Trajetória secundária de equilíbrio.

De acordo com a Figura 8 após a bifurcação representada pelo ponto C o deslocamento u ₃

do vértice da treliça vai aumentando até alcançar um máximo de ±0.316. Logo após atingir o máximo o deslocamento vai diminuindo até anular-se no ponto D, que é quando a treliça retorna ao plano (x, y). É importante destacar que o que possibilita, do ponto de vista físico, que a treliça saia e retorne ao plano (x, y) é a presença da mola na direção do eixo z colocada no vértice da treliça.

Na Figura 8 mostra-se a trajetória de equilíbrio secundária que é obtida levando em conta a equação (5.7a) e a restrição εG = $-\frac{1}{2}\kappa$, cuja expressão é dada por

${\lambda }_2=-2{\epsilon }_G(\alpha -{\mu }_2)=\kappa (\alpha -{\mu }_2)$ (5.14)

Conforme mostra-se na Figura 8, os pontos C e D onde há a intersecção entre as trajetórias primária e secundária de equilíbrio denominam-se pontos de bifurcação. Nestes pontos a matriz de rigidez tangente dada pela equação (5.8) também torna-se singular. No ponto de bifurcação C a estrutura muda para a trajetória secundária de equilíbrio enquanto que no ponto de bifurcação D ela retorna à trajetória primária de equílibrio. Substituindo os valores de μ ₂ , obtidos ao impor que μ ₃ = 0 na equação (5.13), na equação (5.14) determinam-se os pontos de bifurcação da seguinte maneira

$$\begin{gathered} {\mu }_2^c=\alpha -\sqrt{{\alpha }^2-\kappa }\Rightarrow {\lambda }_2^c=\kappa \sqrt{{\alpha }^2-\kappa } \\ {\mu }_2^d=\alpha -\sqrt{{\alpha }^2-\kappa }\Rightarrow {\lambda }_2^d=-\kappa \sqrt{{\alpha }^2-\kappa } \end{gathered}$$ (5.15)

Note-se para que haja flambagem deve verificar-se α ² − κ > 0. Por fim, considerando as trajetórias de equilíbrio mostradas nas Figuras 6 e 8, respectivamente, pode-se concluir que a trajetória de equilíbrio primária da treliça analisada possui 4 pontos críticos, 2 pontos limites e 2 pontos de bifurcação. Como foi dito anteriormente, a matriz de rigidez tangente dada pela equação (5.9) para as condições γ = 0 e μ ₃ = 0 torna-se singular nos pontos críticos, isto é, seu determinante anula-se. É importante destacar que a sequência de ocorrência desses pontos críticos depende da relação entre os parâmetros adimensionais α e κ, respectivamente. Assim para que ocorra o ponto de bifurcação C antes do ponto limite A é necessário que ${\mu }_2^c<{\mu }_2^a$. Para que o ponto de bifurcação C ocorra depois do ponto limite A é necessário que ${\mu }_2^c>{\mu }_2^a\text{e}\kappa <{\alpha }^2$ e para que haja bifurcação é necessário que κ < α ² Por último, para que o ponto de bifurcação C coincida com o ponto limite A é necessário que ${\mu }_2^c={\mu }_2^a$. Essas condições estão sumarizadas na equação abaixo

$$\begin{gathered} 0<\kappa <\frac{2}{3}{\alpha }^2\text{o ponto de bifurcação ocorre antes do ponto limite}; \\ \frac{2}{3}{\mathrm{\alpha }}^2<\mathrm{\kappa }<{\mathrm{\alpha }}^2\text{o ponto de bifurcação ocorre depois do ponto limite}; \\ {\mathrm{\alpha }}^2<\mathrm{\kappa }\text{não ocorre bifurcação}; \\ \mathrm{\kappa }<\frac{2}{3}{\mathrm{\alpha }}^2\text{o ponto de bifurcação coincide com o ponto limite}. \end{gathered}$$ (5.16)

Mostram-se na Figura 9 as trajetórias de equilíbrio apresentadas nas Figuras 6 e 7 em um sistema de referência tridimensional. Nota-se nesta figura que o caminho de equilíbrio secundário projeta-se como um círculo no plano (u ₂ , u ₃ ) mostrado na Figura 7, e como uma reta que liga os pontos C e D no plano (u ₂ , f ₂ ) mostrado na Figura 8. Como pode ser observado nesta figura, nos pontos de bifurcação C e D há a interseção entre as trajetórias de equilíbrio primária e secundária, respectivamente. Além disso, nos pontos E e F o deslocamento u ₃ atinge o valores máximos absolutos de mesma magnitude. Os pontos limites A e B são extremos em relação à carga e pertencem à trajetória primária de equilíbrio.

Figura 9:
Trajetórias de equilíbrio em 3D.

O próximo passo é o cálculo do determinante da matriz de rigidez tangente para a trajetória de equilíbrio primária, ou seja, considerando-se a treliça espacial contida no plano (x, y). Levando em consideração a equação (5.9), define-se a matriz de rigidez tangente para as condições γ = 0 e μ ₃ = 0 como

$\mathit{K}=\left[\begin{array}{cc}2(\alpha -{\mu }_2{)}^2+2{\epsilon }_G& 0\\ 0& 2{\epsilon }_G+\kappa \end{array}\right],\text{com}{\epsilon }_G=-{\mu }_2(\alpha -\frac{1}{2}{\mu }_2),$ (5.17)

cujo determinante em função dos parâmetros α, κ e μ ₂ e´ dado por

$$\begin{gathered} det\mathit{K}=\left(2\right(\alpha -{\mu }_2{)}^2+2{\epsilon }_G\left)\right(2{\epsilon }_G+\kappa ) \\ det\mathit{K}=3{\mu }_2^4-12\alpha {\mu }_2^3+(14{\alpha }^2+3\kappa ){\mu }_2^2-(4{\alpha }^3+6\alpha \kappa ){\mu }_2+2{\alpha }^2\kappa =0. \end{gathered}$$ (5.18)

Trata-se de um polinômio de grau 4 com 4 raízes reais como mostra-se na Figura 10. Essas raízes representam os quatros pontos críticos definidos nas equações (5.12) e (5.15), sendo 2 pontos limites (A, B) e 2 pontos de bifurcação (C, D).

Figura 10:
Determinante da matriz de rigidez tangente.

6 ANÁLISE NUMÉRICA

Para realizar a análise numérica do exemplo apresentado no item anterior utilizou-se um programa escrito em linguagem, Fortran90 denominado gnla_truss.f90 escrito pelos autores deste artigo. O programa faz a análise incremental-iterativa utilizando o método de Newton-Raphson conjuntamente com o método de comprimento de arco. Utiliza-se a formulação Lagrangeana Total do elemento de barra bi-articulado 3D descrito nos itens 2, 3 e 4. Na implementação computacional utilizou-se a equação (3.2d) para o cálculo do vetor de forças internas expresso em coordenadas globais, e a equação (4.4) para o cálculo da matriz de rigidez tangente cujos coeficientes também estão expressos em coordenadas globais. Para capturar a trajetória secundária de equilíbrio considerou-se a treliça espacial com uma imperfeição na direção e sentido positivo do eixo z da ordem de grandeza de 10⁻³, isto é, c = 0.001. Discretizou-se a treliça espacial com 2 elementos de barra bi-articulado 3D. Adotou-se o módulo de elasticidade de E = 100 e a área da seção transversal de A = 1. Novamente, adotaram-se os seguintes valores para as variáveis a = 1, b = 2 e κ = 0.1. Portanto, de acordo com a equação (5.1) a rigidez da mola é dada por k = $\kappa =\frac{EA}{l_0}\kappa \Rightarrow k=2\sqrt{5}$. Para o cálculo da rigidez da mola não se levou em consideração a imperfeição c no cálculo do comprimento indeformado l ₀

Nas Figuras 11a e 11b mostram-se a evolução dos deslocamentos do vértice da treliça na direção dos eixos y e z, respectivamente. Na Figura 11a comparam-se os resultados analíticos com os numéricos para a trajetória secundária de equilíbrio expressa pela curva μ ₂ ×λ ₂ , como pode ser observado há uma boa concordância entre ambos. Como já foi dito anteriormente os pontos C e D são pontos de bifurcação e a reta que os une é a projeção da trajetória secundária de equilíbrio projetada no plano (μ ₂ ×λ ₂ ), veja a Figura 9. Na Figura 11b mostra-se a curva entre a carga aplicada no vértice da treliça na direção do eixo y e o seu deslocamento na direção do eixo z quando ocorre a flambagem. Nesta figura o ponto E representa o máximo deslocamento que o vértice alcança fora do plano (x, y) na direção da imperfeição c. Também, observa-se nesta figura que após o deslocamento u ₃ alcançar o máximo em E seu valor diminui até zero. Por fim, na Figura 11c mostra-se a relação entre os deslocamentos do vértice na direção dos eixos y e z, respectivamente. Observe que inicialmente o deslocamento u ₃ é zero, quando ocorre a bifurcação em C há uma transição suave para valores diferentes de zero devido á magnitude adotada para imperfeição, c = 0.001. Novamente, mostra-se nesta figura, que este deslocamento atinge um máximo em E . Após este ponto o deslocamento u ₃ vai diminuindo até o valor zero no ponto D, que é o segundo ponto de bifurcação. No ponto de bifurcação D a treliça volta para o caminho de equilíbrio primário. Também na Figura 11c comparam-se os resultados numéricos com os analíticos onde pode-se observar a boa coincidência entre ambos resultados. Para obter as trajetórias secundárias de equilíbrio mostradas nas Figuras 11a, 11b e 11c, respectivamente, utilizou-se um processo incremental-iterativo com comprimento de arco constante de 0.025 para 120 passos de carga. Adotou-se uma tolerância para a convergência de 10⁻⁵. O número médio de iterações foi de 2.1.

Figura 11:
Trajetórias secundárias de equilíbrio. a) μ ₂ ×λ ₂ b) μ ₂ × λ ₂ c) μ ₂ × μ ₃.

Para detectar e classificar os pontos singulares no caminho de equilíbrio primário serão adotadas neste trabalho somente duas funções de prova dentre várias que estão descritas em detalhes nas referências ⁴, ⁵. A primeira função de prova adotada será o número de pivôs positivos/negativos da matriz de rigidez tangente. Quando o equilíbrio da estrutura passa da fase estável para a fase instável a matriz de rigidez tangente deixa de ser positiva definida e passa a ser indefinida, isto é, esta matriz que possuia somente pivôs positivos anteriormente ao ponto crítico passa a ter também pivôs negativos após esse ponto crítico. No programa de elementos finitos desenvolvido pelos autores para realizar a análise numérica aquí apresentada utiliza para a resolução do sistema de equações lineares a seguinte decomposição da matriz de rigidez tangente

$$\begin{gathered} \left\{\begin{array}{l}\mathit{K}=\mathit{L}\mathit{D}{\mathit{L}}^T \\ \text{onde}\mathit{L}\text{é uma matriz triangular inferior}\\ \text{com}{L}_{ii}=1e{\mathit{L}}^T={\mathit{L}}^{-1}\end{array}\right. \end{gathered}$$ (6.1)

e D é uma matriz diagonal

$\mathit{D}=diag\left[D_{ii}\right]=\left[\begin{array}{cccc}{D}_{11}& 0& \cdots & 0\\ 0& D_{22}& \cdots & 0\\ ⋮& ⋮& \ddots & ⋮\\ 0& 0& \cdots & D_{mm}\end{array}\right]$ (6.2)

Desta maneira define-se esta função de prova como

$$\begin{gathered} \tau =\left\{\begin{array}{l}+1\text{se o número de pivôs positivos (ou negativos) de}\mathit{D} \\ \text{é igual entre os passos de carga}n-1en\\ -1\text{se o número de pivôs positivos (ou negativos) de}\mathit{D} \\ \text{é diferente entre os passos de carga}n-1en\end{array}\right. \end{gathered}$$ (6.3)

A segunda função de prova é específica para determinar pontos limites. Quando a estrutura atinge um ponto limite ela perde a sua capacidade portante. Assim, além da matriz de rigidez tangente tornar-se singular no ponto limite a estrutura não suporta acréscimo de carregamento após esse ponto crítico. Esta função de prova será denominada de parâmetro de rigidez do sistema estrutural e é definida como

$k=\frac{\mathrm{\Delta }{\mathit{u}}_p^t\mathit{p}}{\mathrm{\Delta }{\mathit{u}}_p^t\mathrm{\Delta }{\mathit{u}}_{\mathbf{p}}},\mathrm{\Delta }{\mathit{u}}_{\mathrm{p}}={\mathit{K}}^{-1}\mathit{p},$ (6.4)

onde ∆u _p é o incremento do vetor de deslocamentos nodais da fase preditora e p é o vetor de forças nodais externas. Normaliza-se esta função de prova da seguinte maneira: S _p =$S_p=\frac{k^n}{k^0}$ onde k ⁰ é o parâmetro de rigidez avaliado na primeira iteração do primeiro passo de carga e k” é o parâmetro de rigidez avaliado na fase preditora do passo de carga atual. O cálculo de S _p é realizado na fase preditora de cada passo de carga. Quando se alcança o ponto limite o parâmetro de rigidez k e consequentemente S _p tendem a zero, enquanto que para o ponto de bifurcação tanto k quanto S _p atingem um valor arbitrário diferente de zero. Portanto, utilizam-se o parâmetro de rigidez S _p e a função de prova τ para detectar e classificar os pontos singulares da seguinte maneira

$$\begin{gathered} \text{ponto limite:}\tau =0S_p=0 \\ \text{ponto de bifurcação:}\tau =0S_p\ne 0 \end{gathered}$$ (6.5)

Além das funções de provas dadas na equação (6.5), no programa gnla_truss.f90 estão implementadas as demais funções de provas descritas nas referências ⁴, ⁵. Para detectar os quatros pontos singulares no caminho primário de equilíbrio analisou-se a treliça espacial sem a imperfeição, isto é, c = 0. Os valores assumidos para as demais variáveis foram E = 100 A=1, a = 1, b = 2, k = 0.1 e k = $2\sqrt{5}$. Novamente, a treliça espacial foi discretizada com 2 elementos de barra bi-articulados 3D. Neste caso, utilizou-se um processo incremental-iterativo com comprimento de arco constante de 0.025 para 90 passos de carga. Adotou-se uma tolerância para a convergência de 10⁻⁵. O número médio de iterações foi de 2.0. Na Figura 12a mostram-se as curvas do determinante da matriz de rigidez tangente obtidas analíticamente e numericamente, respectivamente. Nesta figura pode-se observar a boa concordância entre esses resultados. Por outro lado, observa-se também que há 4 pontos singulares, os pontos A, B, C e D, onde o determinante anula-se. Somente com a condição de detK = 0, como mostrado na Figura 12a, não há como classificar esses pontos. Na Figura 12 mostram-se os resultados numéricos obtidos para as funções de prova definidas na equações (6.3) e (6.4). Observe-se que para os pontos A e B ambas funções de prova tornam-se nulas, que de acordo com as condições expressas na equação (6.5), tratam-se de pontos limites. Entretanto para os pontos C e D somente a função de prova sinal do pivô anula-se enquanto que o parâmetro de rigidez é diferente de zero, portanto tratam-se de pontos de bifurcação.

Figura 12:
a) Determinante da matriz de rigidez tangente. b) Funções de prova.

7 CONCLUSÕES

De maneira sucinta pode-se dizer que neste trabalho foi descrito de forma objetiva, concisa e detalhada a abordagem analítica para a detecção, classificação e sequência de occorrência de pontos críticos na trajetória primária de equilíbrio de um sistema físico simples. Determinaram-se as condições geométricas e físicas para as quais ocorrem a coalescência entre os pontos limites e os pontos de bifurcação. A condição de coalescência de pontos críticos deve ser um requisito em estudos de otimização relativos à estabilidade estrutural. A presença de pontos críticos na trajetória de equilíbrio traduz-se na singularidade da matriz de rigidez tangente. Neste quesito as funções de provas apresentadas neste trabalho, além da simplicidade conceitual e da facilidade de implementação computacional, mostraram-se capazes de detectar e classificar os pontos críticos. Portanto, é necessário levar em consideração a não linearidade geométrica na análise da estabilidade de equilíbrio de sistemas estruturais. É importante destacar que o principal objetivo deste trabalho foi demonstrar a necessidade de compreender melhor os fenômenos não lineares para projetar sistemas estruturais mais seguros. Por fim, outro ponto importante a comentar é que para uma melhor simulação dos fenômenos não lineares observados nos mais diferentes sistemas físicos é necessário um bom domínio da matemática aplicada e computacional.

REFERÊNCIAS

• ¹
  T. Belytschko, W.K. Liu & B. Moran. “Nonlinear finite elements for continua and strucutures”, First edition, John Wiley, 1 (2000), 650 p.
• ²
  J. Bonet & R.D. Wood. “Nonlinear continuum mechanics for finite element analysis”, 2nd Edition, Cambridge University Press, 1 (2008), 318 p.
• ³
  R. Borst, M.A. Crisfield, J.J.C. Remmers & C.V. Verhoosel. “Non-linear finite element analysis of solids and structures”, 2nd Edition, Wiley, 1 (2012), 516 p.
• ⁴
  M.A. Crisfield. “Non-linear finite element analysis of solids and structures”, Essentials, John Wiley, 1 (1991), 345 p.
• ⁵
  M.A. Crisfield. “Non-linear finite element analysis of solids and structures”, Advanced Topics, John Wiley, 2 (1997), 494 p.
• ⁶
  J.F. Doyle. “Nonlinear analysis of thin-walled structures. Statics, Dynamics and Stability”, Springer, 1 (2001), 511 p.
• ⁷
  K. Hashiguchi & Y. Yamakawa. “Introduction to finite strain theory for continuum elasto-plasticity”, Wiley, 1 (2013), 417 p.
• ⁸
  M. Kojić & K.J. Bathe. “Inelastic analysis of solids and structures”, Springer, 1 (2005), 414 p.
• ⁹
  S. Krenk. “Non-linear modeling and analysis of solids and structures”, Cambridge University Press, 1 (2009), 349 p.
• ¹⁰
  E.A.S. Neto, D. Perić & D.R.J. Owen. “Computational methods for plasticity. Theory and aplications”, John Wiley, 1 (2008), 791 p.
• ¹¹
  J.C. Simo & T.J.R. Hughes. “Computational Inelasticity”, Springer, 1 (1998), 392 p.
• ¹²
  G.Z. Voyiadjis & P. Woelke. “Elasto-plastic and damage analysis of plates and shells”, Springer, 1 (2010), 208 p.
• ¹³
  P. Wriggers. “Nonlinear finite element methods”, Springer, 1 (2008), 559 p.

Datas de Publicação

• Publicação nesta coleção
  Dez 2017

Histórico

• Recebido
  19 Jan 2016
• Aceito
  04 Set 2017