ANÁLISE DA INSTABILIDADE DINÂMICA DE ESTRUTURAS …livros01.livrosgratis.com.br/cp074018.pdf · CARVALHO, EULHER CHAVES (2008). Análise da Instabilidade Dinâmica de Estruturas

UNIVERSIDADE FEDERAL DE GOIÁS

ESCOLA DE ENGENHARIA CIVIL

EULHER CHAVES CARVALHO

ANÁLISE DA INSTABILIDADE DINÂMICA DE ESTRUTURAS ESTAIADAS

Goiânia

2008

Livros Grátis

http://www.livrosgratis.com.br

Milhares de livros grátis para download.



Dissertação apresentada ao Curso de Mestrado em Engenharia Civil da Universidade Federal de Goiás, como parte integrante dos requisitos para obtenção do título de Mestre em Engenharia Civil. Área de concentração: Estruturas. Orientador: Professor D. Sc. Zenón José Guzmán Núñez Del Prado.

Goiânia

2008

DADOS INTERNACIONAIS DE CATALOGAÇÃO NA PUBLICAÇÃO (CIP)

(GPT/BC/UFG)

Carvalho, Eulher Chaves.

C331a Análise da instabilidade dinâmica de estruturas estaiadas.

[masculino] / Eulher Chaves Carvalho – 2008.

xxi, 121 f.: il.; figs., tabs., grafs.

Orientador: Prof. Dr. Zenón José Guzmán Núñez Del Prado.

Dissertação (Mestrado) – Universidade Federal de Goiás, Escola de

Engenharia Civil, 2008.

Bibliografia: f.115.

Inclui lista de figuras, tabelas, símbolos e abreviaturas.

Apêndice.

1. Estruturas Estaiadas. 2. Imperfeições Geométricas. 3. Instabilidade

Dinâmica. I. Núñez Del Prado, Zenón José. II. Universidade Federal de Goiás,

Escola de Engenharia Civil. III. Título.

CDU: 624.072

REFERÊNCIA BIBLIOGRÁFICA

CARVALHO, EULHER CHAVES (2008). Análise da Instabilidade Dinâmica de Estruturas

Estaiadas. Publicação/2008, Curso de Mestrado em Engenharia Civil, Universidade Federal

de Goiás, 121 f.

CESSÃO DE DIREITOS

Nome do Autor: Eulher Chaves Carvalho.

Título da Dissertação: Análise da Instabilidade Dinâmica de Estruturas Estaiadas.

Grau/Ano: Mestrando/2008.

Na qualidade de titular dos direitos de autor, autorizo a Universidade Federal de Goiás – UFG

a disponibilizar gratuitamente através da Biblioteca Digital de Teses e Dissertações –

BDTD/UFG, sem ressarcimento dos direitos autorais, de acordo com a Lei nº 9610/98, este

documento, para fins de leitura, impressão e/ou download, a título de divulgação da produção

científica brasileira.

__________________________________________

Eulher Chaves Carvalho - [email protected]



Dissertação defendida no Curso de Mestrado em Engenharia Civil da Universidade Federal de Goiás, para obtenção do grau de Mestre, aprovada em 15 de agosto de 2008, pela Banca Examinadora constituída pelos seguintes professores:

_________________________________________________________ Professor D. Sc. Zenón José Guzmán Núñez Del Prado (UFG) Orientador.

_________________________________________________________ Professora D. Sc. Sylvia Regina Mesquita de Almeida (UFG) Examinador Interno.

_________________________________________________________ Professor D. Sc. Paulo Batista Gonçalves (PUC-RIO) Examinador Externo.

Termo de Ciência e de Autorização para Disponibilizar as Teses e Dissertações Eletrônicas (TEDE) na Biblioteca Digital da UFG

Na qualidade de titular dos direitos de autor, autorizo a Universidade Federal de Goiás –

UFG a disponibilizar gratuitamente através da Biblioteca Digital de Teses e Dissertações – BDTD/UFG, sem ressarcimento dos direitos autorais, de acordo com a Lei nº 9610/98, o documento conforme permissões assinaladas abaixo, para fins de leitura, impressão e/ou download, a título de divulgação da produção científica brasileira, a partir desta data.

1. Identificação do material bibliográfico: [ X ] Dissertação [ X ] Tese

2. Identificação da Tese ou Dissertação

Autor (a): Eulher Chaves Carvalho CPF: 88082555149 E-mail: [email protected] Seu e-mail pode ser disponibilizado na página? [ X ]Sim [ X ] Não

Vínculo Empregatício do autor Bolsista

Agência de fomento: Conselho Nacional de Desenvolvimento

Científico e Tecnológico Sigla: CNPq

País: Brasil UF: GO CNPJ: Título: Análise da Instabilidade Dinâmica de Estruturas Estaiadas

Palavras - chave: Estruturas Estaiadas, Instabilidade Dinâmica e Imperfeições Geométricas

Título em outra língua: Dynamic Stability Analysis of Cable Stayed Structures Palavra - chave em outra língua:

Cable Stayed Structures, Dynamic Stability and Geometric Imperfections

Área de concentração: Estruturas Data defesa: (dd/mm/aaaa) 15 de agosto de 2008 Programa de Pós-Graduação: Engenharia Civil Orientador (a): Professor D. Sc. Zenón José Guzmán Núñez Del Prado CPF: 053370937-70 E-mail: [email protected] Co-orientador (a): CPF: E-mail: 3. Informações de acesso ao documento: Liberação para disponibilização?1 [ X ] total [ X ] parcial Em caso de disponibilização parcial, assinale as permissões: [ ] Capítulos. Especifique: _____________________________________________________ [ ] Outras restrições: _________________________________________________________

Havendo concordância com a disponibilização eletrônica, torna-se imprescindível o envio

do(s) arquivo(s) em formato digital PDF ou DOC da tese ou dissertação. O Sistema da Biblioteca Digital de Teses e Dissertações garante aos autores, que os

arquivos contendo eletronicamente as teses e ou dissertações, antes de sua disponibilização, receberão procedimentos de segurança, criptografia (para não permitir cópia e extração de conteúdo, permitindo apenas impressão fraca) usando o padrão do Acrobat.

________________________________________ Data: 15 / 10 / 2008 Eulher Chaves Carvalho (Autor)

1 Em caso de restrição, esta poderá ser mantida por até um ano a partir da data de defesa. A extensão deste prazo suscita

justificativa junto à coordenação do curso. Todo resumo e metadados ficarão sempre disponibilizados.

AGRADECIMENTOS

A minha família e familiares, por tudo.

Ao professor Zenón, pela sábia orientação e ajuda amiga durante essa empreitada.

Aos professores do curso, em especial, a Sylvia e ao Fred, pela colaboração e pelos

ensinamentos.

Ao amigo Luciano, pela enriquecedora convivência.

A banca examinadora, por aceitar o convite de colaborar com este trabalho.

Aos funcionários da escola.

Ao CNPq, pelo apoio financeiro.

RESUMO

As estruturas estaiadas são freqüentemente empregadas na concepção de torres elevadas, coberturas de grandes espaços, estruturas off-shore e na sustentação de tabuleiros de pontes. Por serem leves, esbeltas e flexíveis, estão sujeitas a ter perda de estabilidade sob a ação de cargas estáticas e dinâmicas bem como podem ser sensíveis a imperfeições geométricas iniciais. Neste trabalho, analisa-se a perda de estabilidade estática e dinâmica de estruturas estaiadas planas, com e sem imperfeições geométricas iniciais. Duas formulações são utilizadas, a saber: a) Modelo de elementos finitos não-lineares usando um referencial Lagrangiano atualizado e, b) Modelo discreto não-linear considerando os cabos como molas lineares extensíveis.

No modelo de elementos finitos, a solução das equações não-lineares é realizada através do método de Newton-Raphson associado à técnica do comprimento de arco e a integração no tempo é realizada usando o método de Newmark. São apresentados exemplos de validação e é analisada a influência das imperfeições geométricas iniciais e do tensionamento dos cabos nos caminhos pós-críticos de torres estaiadas sob a ação de cargas axiais estáticas. É estudada também, a influência das imperfeições geométricas nas freqüências naturais e nos modos de vibração de torres. Utilizando o critério de Budianski, é avaliada a perda de estabilidade de torres quando submetidas à ação de cargas súbitas e harmônicas.

No modelo discreto, as equações no tempo são integradas através do métodos de Runge-Kutta e o método da força bruta é utilizado para calcular os diagramas de bifurcação e bacias de atração. São estudados os caminhos pós-críticos, as freqüências naturais bem como as fronteiras de estabilidade, mecanismos de escape e a evolução da estabilidade global através das bacias de atração quando as torres estaiadas são submetidas a cargas axiais harmônicas.

Pode-se observar a grande influencia do tensionamento dos cabos e das imperfeições geométricas iniciais no valor das cargas críticas, nos caminhos pós-críticos, nos mecanismos de escape e na estabilidade global das torres estaiadas. Os resultados obtidos podem servir de base a engenheiros e projetistas para a melhor compreensão do comportamento estático e dinâmico destas estruturas.

ABSTRACT

Cable stayed structures are widely used to build high towers, cover wide spans,

off-shore structures and to support long bridges. As light and flexible structures, cable stayed

structures, can loss stability under static and dynamic loads and could be sensitive to initial

geometric imperfections. In this work, the static and dynamic stability of both, perfect and

imperfect, plane cable stayed structures is studied. Two formulations are used a) Finite

element non-linear model using a Lagrangean actualized referential and b) Discrete non-linear

model considering extensible linear springs.

In the finite element model, the non-linear equations are solved using the Newton-Raphson method associated to the arc-length technique and the Newmark method is used to calculate the time response. Validation examples are presented and the influence of initial geometric imperfections and cable tensioning is studied when stayed towers are subjected to axial static loads. The influence of geometric imperfections on modes and natural frequencies is also studied. Using the Budianski´s criterion, the loss of stability under sudden and harmonic loads is analyzed.

In the discrete model, the non-linear equations are solved using the Runge-Kutta method and the brute force method is used to calculate the bifurcations diagrams and basins of attraction. The post-critical paths, natural frequencies, instability boundaries, escape mechanisms and the evolution of the global stability are studied when stayed towers are subjected to axial harmonic loads.

The high influence of cable tensioning and initial imperfections on critical loads, non-linear paths, escape mechanisms and basins of attraction can be observed. The results obtained can be used as a base for design criterions for engineers allowing a better understanding of the static and dynamic behavior of cable stayed structures.

vi

SUMÁRIO

SUMÁRIO VI

LISTA DE FIGURAS IX

LISTA DE TABELAS XIV

LISTA DE SÍMBOLOS E ABREVIATURAS XV

1 INTRODUÇÃO 1

1.1 JUSTIFICATIVA 1

1.2 REVISÃO BIBLIOGRÁFICA 3

1.3 OBJETIVOS 9

1.4 ESCOPO 10

2 FORMULAÇÃO MATEMÁTICA 11

2.1 INTRODUÇÃO 11

2.2 SISTEMAS DE REFERÊNCIAS 12

2.3 RELAÇÃO DEFORMAÇÃO-DESLOCAMENTO 14

2.4 RELAÇÕES CONSTITUTIVAS 16

2.5 FUNCIONAL DE ENERGIA 18

2.6 FORMULAÇÃO DO ELEMENTO FINITO 21

vii

2.6.1 FORMULAÇÃO DO ELEMENTO DE PÓRTICO PLANO 21

2.6.2 FORMULAÇÃO DO ELEMENTO DE CABO 23

2.7 MATRIZ DE RIGIDEZ DOS ELEMENTOS 25

2.8 VETOR DE FORÇAS INTERNAS 26

2.9 ACOPLAMENTO DOS ELEMENTOS DE CABO E DE PÓRTICO 29

2.10 MATRIZ DE MASSA E DE AMORTECIMENTO 29

2.11 AMORTECIMENTO 30

2.12 EQUAÇÃO DINÂMICA NÃO-LINEAR DE EQUILÍBRIO 31

2.13 FREQÜÊNCIAS NATURAIS 31

2.14 MODELO DISCRETO 32

2.14.1 FORMULAÇÃO ESTÁTICA 33

2.14.2 FORMULAÇÃO DINÂMICA 35

3 MÉTODOS NUMÉRICOS 39

3.1 INTRODUÇÃO 39

3.2 PRINCÍPIOS DA TEORIA DE ESTABILIDADE ESTRUTURAL 40

3.3 DETERMINAÇÃO DA TRAJETÓRIA DE EQUILÍBRIO 41

3.4 MÉTODO DE SOLUÇÃO PARA ANÁLISE NÃO-LINEAR 43

3.4.1 MÉTODO DE NEWTON-RAPHSON 44

3.4.2 SOLUÇÃO INCREMENTAL PREDITA 45

3.4.3 CICLO DE ITERAÇÕES 47

3.4.4 MÉTODO DE NEWTON-RAPHSON ASSOCIADO À TÉCNICA DO COMPRIMENTO

DE ARCO 49

3.4.5 CRITÉRIOS DE CONVERGÊNCIA 52

3.5 MÉTODOS DE SOLUÇÃO PARA ANÁLISE DINÂMICA 53

3.5.1 CONCEITOS BÁSICOS 53

3.5.2 EVOLUÇÃO DA RESPOSTA NO TEMPO E PLANO FASE 53

3.5.3 SOLUÇÕES DE EQUILÍBRIO 58

3.5.4 SOLUÇÕES PERIÓDICAS 58

3.5.5 SEÇÕES DE POINCARÉ 58

3.5.6 BIFURCAÇÃO DAS SOLUÇÕES DE EQUILÍBRIO 60

3.5.7 FRONTEIRAS DE INSTABILIDADE E ESCAPE DE SISTEMAS FORÇADOS 64

viii

3.5.8 ANÁLISE GLOBAL DE SISTEMAS DINÂMICOS – DETERMINAÇÃO DAS BACIAS DE

ATRAÇÃO 65

3.6 IMPLEMENTAÇÃO COMPUTACIONAL 66

4 RESULTADOS NUMÉRICOS 70

4.1 INTRODUÇÃO 70

4.2 ANÁLISE ESTÁTICA NÃO-LINEAR VIA MÉTODO DOS ELEMENTOS FINITOS 70

4.2.1 VALIDAÇÃO DO MÉTODO 70

4.2.2 EFEITO DAS IMPERFEIÇÕES GEOMÉTRICAS INICIAIS 77

4.3 ANÁLISE DINÂMICA LINEAR VIA MÉTODO DOS ELEMENTOS FINITOS 85


4.3.2 EFEITO DAS IMPERFEIÇÕES GEOMÉTRICAS INICIAIS NAS FREQÜÊNCIAS

NATURAIS 87

4.4 ANÁLISE DINÂMICA NÃO-LINEAR VIA MÉTODO DOS ELEMENTOS FINITOS 89

4.4.1 INSTABILIDADE DINÂMICA PARA CARGA SÚBITA 89

4.4.2 INSTABILIDADE DINÂMICA PARA CARGA HARMÔNICA 92

4.5 ANÁLISE ESTÁTICA LINEAR VIA MODELO DISCRETO 99


4.6 ANÁLISE DINÂMICA LINEAR VIA MODELO DISCRETO 100

4.6.1 INSTABILIDADE PARAMÉTRICA 101

4.6.2 MECANISMOS DE ESCAPE 103

4.6.3 EVOLUÇÃO DA ESTABILIDADE GLOBAL 108

5 CONCLUSÕES E SUGESTÕES 112

6 REFERÊNCIAS BIBLIOGRÁFICAS 115

7 APÊNDICE 118

ix

LISTA DE FIGURAS

Figura 1.1 – Torre estaiada Mast Cabauw. ................................................................................. 2

Figura 1.2 – Torre estaiada Zendmast Lopik. ............................................................................. 2

Figura 1.3 – Ponte estaiada Juscelino Kubitschek. ..................................................................... 2

Figura 1.4 – Ponte suspensa Tsing Ma. ...................................................................................... 4

Figura 1.5 – Ponte suspensa Great Belt. ..................................................................................... 5

Figura 2.1 – Referencial Lagrangiano Total. ............................................................................ 13

Figura 2.2 – Referencial Lagrangiano Atualizado.................................................................... 13

Figura 2.3 – Deformação na flexão: (a) Elemento indeformado; (b) Elemento

deformado; (c) Elemento infinitesimal. .................................................................................... 15

Figura 2.4 – Esforços no elemento. .......................................................................................... 16

Figura 2.5 – Movimento de um corpo tridimensional. ............................................................. 18

Figura 2.6 – Elemento de Pórtico. ............................................................................................ 21

Figura 2.7 – Elemento de treliça. .............................................................................................. 23

Figura 2.8 – Cálculo das forças internas: (a) deslocamentos incrementais; (b)

geometria deformada. ............................................................................................................... 28

Figura 2.9 – Modelo de molas: (a) configuração inicial; (b) configuração

perturbada. ................................................................................................................................ 33

Figura 2.10 – Posição do elemento de barra com comprimento infinitesimal, ds. ................... 36

Figura 3.1 – Trajetória de equilíbrio com salto dinâmico do tipo snap-through...................... 42

Figura 3.2 – Trajetória de equilíbrio com salto dinâmico do tipo snap-back. .......................... 42

Figura 3.3 – Método de Newton-Raphson: (a) Padrão; (b) Modificado. ................................. 46

Figura 3.4 – Solução incremental-iterativa em sistema com um grau de liberdade. ................ 47

x

Figura 3.5 – Técnica do comprimento de arco aplicada a um problema com um

grau de liberdade. ..................................................................................................................... 50

Figura 3.6 – Resposta no tempo (a) e plano fase (b). ............................................................... 57

Figura 3.7 – Trajetória e velocidade no espaço fase 3D. .......................................................... 57

Figura 3.8 – Planos de Poincaré. .............................................................................................. 59

Figura 3.9 - Bifurcações estáticas das soluções de equilíbrio (DEL PRADO,

2001). ........................................................................................................................................ 61

Figura 3.10 - Bifurcação de Hopf: (a) super-crítica; (b) sub-crítica (DEL PRADO,

2001). ........................................................................................................................................ 62

Figura 3.11 - Estrutura das raízes para a equação 0=++ cxxbx &&& , onde b é o

amortecimento efetivo e c a rigidez efetiva. cbD 422 −= . ..................................................... 63

Figura 3.12 – Diagrama de bifurcação típico. .......................................................................... 65

Figura 3.13 – Fluxograma do programa principal para a análise não-linear. ........................... 68

Figura 4.1 – Pórtico de Lee: (a) Sem cabo; (b) Com cabo. ...................................................... 71

Figura 4.2 – Trajetória de equilíbrio para Pórtico de Lee: (a) P x v; (b) P x u. ........................ 71

Figura 4.3 – Trajetória de equilíbrio para o Pórtico de Lee estaiado: (a) P x v; (b)

P x u. ......................................................................................................................................... 72

Figura 4.4 – Deformada do Pórtico de Lee: (a) Sem cabo; (b) Com cabo. .............................. 73

Figura 4.5 – Coluna engastada. ................................................................................................ 74

Figura 4.6 – Trajetória de equilíbrio na extremidade livre (topo) da coluna: (a) P x

v; (b) P x u. ............................................................................................................................... 74

Figura 4.7 – Deformada da coluna engastada-livre. ................................................................. 75

Figura 4.8 – Coluna estaiada com dois cabos a 60°. ................................................................ 75

Figura 4.9 – Trajetória de equilíbrio para a coluna engastada estaiada: (a) P x v;

(b) P x u. ................................................................................................................................... 76

Figura 4.10 – Deformada da coluna estaiada. .......................................................................... 76

Figura 4.11 – Pórtico de Lee com imperfeições iniciais. ......................................................... 77

Figura 4.12 – Trajetória de equilíbrio para Pórtico de Lee com imperfeições

geométricas iniciais: (a) P x v; (b) P x u. .................................................................................. 78

Figura 4.13 – Pórtico de Lee estaiado com imperfeições iniciais. ........................................... 79

Figura 4.14 – Trajetória de equilíbrio para Pórtico de Lee estaiado com

imperfeições geométricas iniciais: (a) P x v; (b) P x u. ............................................................ 79

xi

Figura 4.15 – Deformada do Pórtico de Lee estaiado com imperfeições

geométricas iniciais: (a) Sem cabo; (b) Com cabo. .................................................................. 80

Figura 4.16 – Coluna estaiada imperfeita com dois cabos a 60°. ............................................. 80

Figura 4.17 – Trajetória de equilíbrio para coluna estaiada com geometria inicial

deformada e tensionamento diferenciado: (a) P x v; (b) P x u. ................................................ 81

Figura 4.18 – Deformada coluna estaiada com tensionamento diferenciado. .......................... 82

Figura 4.19 – Coluna estaiada com tensionamentos diferentes. ............................................... 82

Figura 4.20 – Trajetória de equilíbrio para coluna estaiada: (a) P x v; (b) P x u. ..................... 82

Figura 4.21 – Curva de sensibilidade a imperfeições geométricas inicias para T1 =

10 kN e T2 variável. ................................................................................................................. 83

Figura 4.22 – Coluna estaiada imperfeita 0 ≤ δ ≤ 1,0 m. ......................................................... 83

Figura 4.23 – Trajetória de equilíbrio para coluna estaiada: (a) P x v; (b) P x u. ..................... 84

Figura 4.24 – Curva de sensibilidade a imperfeições geométricas inicias para 0 ≤ δ

≤ 0,9 m. ..................................................................................................................................... 84

Figura 4.25 – Coluna engastada de seção constante sem força axial. ...................................... 85

Figura 4.26 – Modos de Vibração Normalizados: (a) 1° Modo; (b) 2° Modo; (c) 3°

Modo. ........................................................................................................................................ 86

Figura 4.27 – Coluna Estaiada: (a) Dois Cabos; (b) Quatro Cabos. ......................................... 87

Figura 4.28 – Variação da freqüência natural em função do ângulo α. .................................... 87

Figura 4.29 – Modos de vibração normalizados para T1 = 10 kN e: (a) T2 = 10

kN; (b) T2 = 9,8 kN; (c) T2 = 9,6 kN; (d) T2 = 8,0 kN. ........................................................... 88

Figura 4.30 – Modos de vibração normalizados: (a) δ = 0,0 m; (b) δ = 0,2 m; (c) δ

= 1,0 m. ..................................................................................................................................... 89

Figura 4.31 – Coluna engastada com carregamento súbito P................................................... 89

Figura 4.32 – Resposta no tempo da coluna engastada para um carregamento

súbito P. .................................................................................................................................... 90

Figura 4.33 – Curva de Budianski para coluna engastada sob carga axial súbita P. ................ 91

Figura 4.34 – Coluna estaiada com dois cabos a 60°. .............................................................. 91

Figura 4.35 – Resposta no tempo da coluna estaiada para um carregamento súbito

P. ............................................................................................................................................... 92

Figura 4.36 – Curva de Budianski para coluna estaiada sob carga axial súbita P.................... 92

Figura 4.37 – Coluna engastada com carregamento harmônico )cos( tPP a Ω= . .................. 93

xii


harmônico )cos( tPP a Ω= . ..................................................................................................... 93

Figura 4.39 – Curva de Budianski para coluna engastada sob carga axial

harmônica )cos( tPP a Ω= . ..................................................................................................... 94

Figura 4.40 – Coluna engastada com carregamento harmônico )2cos( tPP a Ω= . ................ 94


harmônico )2cos( tPP a Ω= . .................................................................................................. 95

Figura 4.42 – Curva de Budianski para coluna engastada sob carga axial

harmônica )2cos( tPP a Ω= . ................................................................................................... 96

Figura 4.43 – Coluna estaiada com dois cabos a 60° e carga axial harmônica

)cos( tPP a Ω= . ....................................................................................................................... 96

Figura 4.44 – Resposta no tempo da coluna estaiada para um carregamento

harmônico )cos( tPP a Ω= . ..................................................................................................... 97

Figura 4.45 – Curva de Budianski para coluna estaiada sob carga axial harmônica

)cos( tPP a Ω= . ....................................................................................................................... 97

Figura 4.46 – Resposta no tempo da coluna estaiada para um carregamento

harmônico )2cos( tPP a Ω= . .................................................................................................. 98


)2cos( tPP a Ω= . .................................................................................................................... 99

Figura 4.48 – Coluna estaiada com duas molas lineares. ......................................................... 99

Figura 4.49 – Caminhos pós-críticos. ..................................................................................... 100

Figura 4.50 – Coluna estaiada com dois cabos a 32° e carga axial periódica

)cos( tPP a Ω= . ..................................................................................................................... 101

Figura 4.51 – Fronteira de Instabilidade Paramétrica............................................................. 102

Figura 4.52 – Mecanismo de escape para 75,0=Ω (seção 1): (a) Diagrama de

bifurcação; (b) Detalhe do diagrama de bifurcação. ............................................................... 103


bifurcação; (b) Resposta no tempo A; (c) Resposta no tempo B; (d) Plano fase e

mapeamento de Poincaré B; (e) Resposta no tempo C; (f) Plano fase e

mapeamento de Poincaré C. ................................................................................................... 104

xiii


bifurcação; (b) Resposta no tempo A; (c) Resposta no tempo B; (d) Plano fase e

mapeamento de Poincaré B. ................................................................................................... 105

Figura 4.55 – Mecanismo de escape para 80,1=Ω (seção 4): diagrama de

bifurcação. .............................................................................................................................. 105


bifurcação; (b) Resposta no tempo do ponto A; (c) Resposta no tempo do ponto B;

(d) Plano fase e mapeamento de Poincaré do ponto B. .......................................................... 106







Figura 4.59 – Diagramas de bifurcação para estudo de erosão da bacia de atração. .............. 109

Figura 4.60 – Evolução da bacia de atração para o diagrama de bifurcação da

Fig.4.59a. ................................................................................................................................ 110

Figura 4.61 – Evolução da bacia de atração para o diagrama de bifurcação da

Fig.4.59b. ................................................................................................................................ 111

Figura 4.62 – Resposta no tempo, plano fase e mapeamento de Poincaré do ponto

A (a); Resposta no tempo, plano fase e mapeamento de Poincaré do ponto B (b)................. 111

xiv

LISTA DE TABELAS

Tabela 4.1 – Propriedades físicas e geométricas do Pórtico de Lee e do cabo. ....................... 71

Tabela 4.2 – Propriedades físicas e geométricas da coluna circular. ....................................... 73

Tabela 4.3 – Coordenadas que definem a geometria inicial do pórtico de Lee da

Fig. 4.11. ................................................................................................................................... 77

Tabela 4.4 – Propriedades físicas e geométricas da coluna de seção constante e

descarregada. ............................................................................................................................ 85

Tabela 4.5 – Comparação das freqüências naturais. ................................................................. 86

Tabela 4.6 – Freqüência natural. .............................................................................................. 88

Tabela 4.7 – Parâmetros físicos e geométricos da coluna e das molas. ................................. 101

xv

LISTA DE SÍMBOLOS E ABREVIATURAS

SÍMBOLOS LATINOS

A Seção transversal do elemento.

a0 - a1 Coeficientes da equação que determina os deslocamentos axiais.

a1 - a2 - a3 - a4 - a5 Fatores de forma ou constantes de interação.

b - h Base e a altura do elemento.

b0 - b1 - b2 - b3 Coeficientes da equação que determina os deslocamentos transversais.

dm Elemento infinitesimal de massa.

ds Comprimento infinitesimal de um elemento da coluna.

E Módulo de elasticidade do material.

0F Força que age na mola.

Fi Força interna.

Fr Carga de referência.

g Força residual.

H Ponto de aplicação das cargas axiais.

H1 ... H6 Funções de interpolação.

h Ponto de fixação das molas.

I - xI Momento de inércia da seção.

Id Número de iterações desejadas para o processo iterativo. tI Número de iterações necessárias para convergência no passo de carga.

K Número de incrementos.

xvi

1K Constante da mola

ijK1 - ijK2 Componentes das parcelas das matrizes de rigidez não-lineares dependentes do deslocamento.

ijLK Componentes da parcela da matriz de rigidez linear.

ijKτ Componentes da parcela da matriz de rigidez geométrica.

k Iteração corrente.

L Comprimento do elemento. tL - t+∆tL Comprimento do elemento na configuração t e tt ∆+ .

Mt Momento fletor na seção transversal na configuração t.

M1 - M2 Momento fletor no primeiro e segundo nó do elemento.

m Massa da barra conectada às molas.

Cm Massa por unidade de comprimento L, do elemento de cabo/treliça.

Pm Massa por unidade de comprimento L, do elemento de viga/coluna.

Nci Número de células que dividem cada coordenada.

Nss Número de iterações do estado permanente.

Nt Número de iterações do estado transiente. np Número de cargas concentradas aplicadas na coluna.

P - p Cargas axiais.

Pt Força axial seção transversal na configuração t.

aP Parcela da carga crítica.

Pb Mapeamento do sistema realizado através do método de Runge-Kutta de quarta ordem.

Pcr Carga crítica.

eP Carga crítica de Euler.

Qt Força cisalhante seção transversal na configuração t.

1xyS Parcela linear das tensões cisalhantes incrementais.

nxyS Parcela não-linear das tensões cisalhantes incrementais.

s Distância do ponto à base da torre.

T Período de vibração natural da estrutura.

t Última configuração de equilíbrio processada.

xvii

tt ∆+ Configuração de equilíbrio processada no passo corrente.

U Energia interna de deformação.

0U Parcela de U∆ que corresponde às forças acumuladas até a configuração de equilíbrio t.

1U - 2U Parcelas de U∆ que resultarão nas matrizes de rigidezes não-lineares.

LU Parcela de U∆ que dá origem à matriz de rigidez linear.

τU Parcela de U∆ que corresponde à influência das deformações iniciais e dá origem à matriz de rigidez geométrica.

u Deslocamentos horizontais do ponto.

iu Deslocamentos generalizados.

u1 - v1 - u2 - v2 Componentes de deslocamento nodais.

v Deslocamentos verticais do ponto.

W Trabalho das forças externas.

x Distância ao primeiro nó do elemento.

YGL - XGL Eixos cartesianos das coordenadas globais. 0Y - 0X Eixos cartesianos das coordenadas na configuração inicial 0=t . tY - tX Eixos cartesianos das coordenadas na configuração deformada tt ∆+ .

y Distância de um ponto em relação à linha neutra da seção.

y - x Sistema local de coordenadas para o elemento.

y Centro de gravidade.

SÍMBOLOS GREGOS

α Inclinação do cabo em relação à base.

α Parâmetro de controle no método da força bruta.

α Ângulo entre o sistema de coordenadas global XGL e YGL e local y e x.

α1 Inclinação da mola em relação à base.

α - β Constantes de amortecimento proporcional de rigidez e de massa.

minα - maxα Limites para variação do parâmetro de controle.

δ Imperfeição geométrica inicial devido aos deslocamentos horizontais.

δ - φ - φ2 Deslocamentos naturais.

xviii

xxe∆ Parcela linear do incremento de deformação axial.

xye∆ Parcela linear do incremento de deformação cisalhante.

XXε∆ Incremento da deformação axial.

XYε∆ Incremento da deformação cisalhante.

t∆l - ∆l Incrementos do comprimento de arco no passo de carga anterior e no passo de carga corrente

l∆ Alongamento da mola.

0l∆ Alongamento inicial da mola.

θl∆ Alongamento causado durante a rotação θ da barra.

∆M Incremento do momento fletor.

xxη∆ Parcela não-linear do incremento de deformação axial.

xyη∆ Parcela não-linear do incremento de deformação cisalhante.

∆Π - π∆ Variação da energia potencial total.

t∆ Incremento do tempo.

T∆ Variação da energia cinética.

u∆ Deslocamento axial de um ponto na linha neutra.

u∆ Deslocamento axial de um ponto distante y da linha neutra da seção.

1u∆ - 2u∆ - 3u∆ Incremento dos deslocamentos nodais.

∆u1 - ∆v1 - ∆θ1 Deslocamento axial, transversal e de rotação do nó 1.

∆u2 - ∆v2 - ∆θ2 Deslocamento axial, transversal e de rotação do nó 2.

U∆ Variação da energia interna de deformação.

v∆ Deslocamento vertical de um ponto.

W∆ Variação do trabalho das forças externas.

ε - ζ Fatores de tolerância.

ξ Razão de amortecimento do sistema.

λ Parâmetro de carga.

λ0 Valor inicial do paramento de carga. 0λ∆ Incremento inicial do paramento de carga.

1λ∆ - 2λ∆ - 3λ∆ Incrementos do parâmetro de carga.

xix

∆λk Incremento do paramento de carga na iteração corrente.

δλ Correção do parâmetro de carga avaliado ao longo do ciclo iterativo.

λk - λk-1 Parâmetro de carga avaliado na iteração k e na iteração anterior.

1Pλ - pξ - 0ξ - 1ξ Coeficientes da equação de caminho pós-crítico adimensional.

Π Energia potencial total.

xxtτ Tensão axial na configuração t.

xytτ Tensão cisalhante na configuração t.

xxtτ∆ Tensões axiais incrementais.

xytτ∆ Tensões cisalhantes incrementais.

Ψ Rotação de corpo rígido.

L Lagrangiano.

ρ Densidade do material.

ω Freqüência de vibração da estrutura.

0ω Freqüência natural de vibração da estrutura.

1ω - 2ω Freqüência natural da estrutura e a máxima freqüência da estrutura.

Ω Parâmetro de freqüência.

VETORES E MATRIZES

C Matriz de amortecimento consistente.

iFt Vetor de forças internas no instante t.

GLFi Vetor de forças internas no sistema global.

)(tIF Forças de inércia dependentes no tempo.

)(tDF Forças de amortecimento dependentes no tempo.

)(tSF Forças elásticas dependentes no tempo.

)(tFr Forças externas dependentes no tempo.

^F Vetor de carregamento efetivo.

Fi Vetor de forças internas.

xx

Fr Vetor de forças de referência.

g Vetor de força residual.

H Matriz das funções de interpolação ou de forma do elemento.

I Matriz identidade.

LK Matriz de rigidez linear.

τK Parcelas não-lineares da matriz de rigidez da estrutura correspondentes à parcela relativa às tensões iniciais.

TK Matriz de rigidez tangente, TK = LK + τK .

^K Matriz de rigidez efetiva.

1K - 2K Parcelas não-lineares da matriz de rigidez da estrutura dependentes dos deslocamentos.

M Matriz de massa consistente.

PM Matriz de massa do elemento de pórtico plano.

CM Matriz de massa do elemento de cabo.

R Matriz de rotação do elemento.

u Vetor de deslocamento total.

u0 Vetor inicial dos deslocamentos nodais.

u - •

u - ••

u Vetores de deslocamento, velocidade e aceleração do sistema.

••

u0 - •

u0 - u0 Vetores iniciais de aceleração, velocidade e deslocamento das coordenadas nodais do sistema.

••∆+ utt -

•∆+ utt - utt ∆+

Vetores de aceleração, velocidade e deslocamento das coordenadas nodais do sistema na configuração tt ∆+ .

Xo Vetor de condições iniciais.

∆d Vetor de deslocamentos incrementais.

i∆F Vetor de incremento das forças internas.

∆R Vetor de forças residuais.

∆u Vetor de deslocamentos nodais incrementais.

n∆u Deslocamentos naturais incrementais do elemento.

GL∆u Vetor de deslocamentos no sistema global de coordenadas.

k∆u Vetor de deslocamentos nodais incrementais avaliados na iteração k. 0∆u Vetor de incremento inicial dos deslocamentos nodais.

xxi

∆y Deslocamentos incrementais.

uδ Vetor de deslocamentos nodais residuais.

δy Deslocamentos incrementais corrigidos.

δug Correção obtida da aplicação do método de Newton-Raphson com a estratégia convencional de incremento do parâmetro de carga constante.

1δur Vetor de deslocamentos nodais, devido o carregamento de referência no primeiro incremento e primeira iteração.

t-1δur Vetor de deslocamentos nodais, devido o carregamento de referência na primeira iteração do incremento anterior.

tδur Vetor de deslocamentos nodais, devido o carregamento de referência na primeira iteração do incremento atual.

δur Vetor de deslocamentos nodais tangentes.

••

θ - •

θ - θ Vetores de aceleração, velocidade e deslocamento das coordenadas livres da estrutura.

Φ Vetor de amplitudes dos deslocamentos nodais.

Ω Matriz espectral, diagonal.

SIGLAS

MEF Método dos Elementos Finitos.

PEP Do inglês Pointwise Equilibrium Polynomial.

RLA Referencial Lagrangiano Atualizado.

RLT Referencial Lagrangiano Total.

1

1 INTRODUÇÃO

1.1 JUSTIFICATIVA

Devido ao avanço dos métodos numéricos aplicados à engenharia, ao emprego de

novos materiais, à aplicação de novas tecnologias e ao uso de técnicas modernas de

construção, as estruturas têm se tornado cada vez mais leves, esbeltas e flexíveis.

Neste contexto, a utilização das estruturas estaiadas tornou-se uma realidade

devido à sua capacidade de suportar grandes cargas mesmo sofrendo grandes deslocamentos,

sendo a solução estrutural freqüentemente empregada na concepção de torres elevadas

(Fig. 1.1)1, antenas de transmissão (Fig. 1.2)2, suportes para coletores de energia solar, nas

estruturas off-shore, na sustentação de tabuleiros de pontes (Fig. 1.3)3, em estruturas especiais

como guindastes e nas coberturas de grandes espaços, tais como estádios, feiras itinerantes,

galpões industriais, entre outras.

Nestes sistemas, a estrutura aporticada ou reticulada, é estaiada por cabos, os

quais, em razão de sua alta eficiência quando submetidos à tração, proporcionam elevada

rigidez ao conjunto, atendendo plenamente concepções arquitetônicas e estruturais exigidas.

Em contrapartida, essa concepção estrutural, impõe fortes não-linearidades

geométricas em seu comportamento. Logo, para ter projetos confiáveis, é necessário conhecer

a fundo o comportamento destas estruturas na presença de cargas estáticas e, principalmente,

de cargas dinâmicas tais como a ação do vento e terremotos, possíveis rupturas de cabos,

ocupação humana, atuação do gelo, entre outros.

1 Fonte: <http://en.wikipedia.org/wiki/KNMI-mast_Cabauw>. Acesso em: 22 ago. 2007. 2 Fonte: <http://www.ronneke.nl/zendmast_lopik/nederlandse_zendmasten.htm>. Acesso em: 22 ago. 2007. 3 Fonte: <http://www.nelsonmarins.com.br/fotos/g2/galeria2.htm>. Acesso em: 06 mai. 2008.

2

Na literatura, existe uma quantidade reduzida de trabalhos que descrevem o

comportamento das estruturas estaiadas submetidas a carregamentos externos dependentes do

tempo. Assim, esse estudo justifica-se pela necessidade de se observar e compreender melhor

os aspectos relacionados à estabilidade dinâmica não-linear dessas estruturas, principalmente

das torres estaiadas.

Figura 1.1 – Torre estaiada Mast Cabauw.

Figura 1.2 – Torre estaiada Zendmast Lopik.

Figura 1.3 – Ponte estaiada Juscelino Kubitschek.

3

1.2 REVISÃO BIBLIOGRÁFICA

Existem na literatura, trabalhos mostrando a interação entre elementos de pórticos

e cabos. Grande parte destes contempla casos específicos e as conclusões são baseadas

exclusivamente nos dados obtidos da solução numérica destes modelos. A seguir, são

apresentados alguns trabalhos relacionados a este tema.

Neves (1990) apresentou uma ferramenta para análise estática não-linear

geométrica e dinâmica linear de pontes estaiadas. Utilizando um modelo tridimensional de

ponte estaiada com três vãos, o autor demonstrou que a modelagem requer considerações na

mudança da geometria, que os cabos influenciam de maneira significativa a resposta estrutural

devido à sua não-linearidade e que não é necessário utilizar elementos mais refinados que os

de cabo/treliça para discretizar os cabos. Em termos de características dinâmicas, ficou

demonstrado que a influência da tensão inicial nos cabos é tanto maior quanto maior for a

flexibilidade do sistema.

Kahla (1997), para compreender a resposta dinâmica não-linear de torres

estaiadas, modelou uma torre tridimensional com 152,4 m de altura e simulou,

simultaneamente, a ação do vento com a ruptura de um dos cabos. Nesta análise, foram

utilizados elementos de treliça espacial na modelagem do mastro e de catenária elástica

tridimensional para os cabos. Os resultados evidenciaram a falha dos cabos na compressão e

que a estrutura colapsou quando um conjunto de cabos perdeu sua capacidade de suportar as

solicitações de tração.

Xu et al. (1997) estabeleceram um elemento finito tridimensional para analise

dinâmica do sistema torre-cabo da ponte suspensa Tsing Ma (Fig. 1.4)1. As duas torres da

ponte foram modeladas com elementos de viga tridimensional de Timonshenko e os cabos

com elementos de cabo tridimensional de três nós. Os resultados foram verificados com as

propriedades dinâmicas do sistema “in loco” e mostraram que, para freqüências naturais

baixas, os modos de vibração do sistema podem ser razoavelmente separados em in-plane e

out-of-plane. Mostraram ainda que as interações dinâmicas entre as torres e os cabos são

significativas para freqüências naturais globais e que existe um grande número de freqüências

naturais locais em que apenas os cabos vibram e as torres permanecem imóveis ou com

movimentos relativamente pequenos.

1 Fonte: <http://depedraecal.blogspot.com/2005/04/ponte-suspensa-tsing-ma-um-feito.html>. Acesso em: 31

jul. 2007.

4

No que se refere à modelagem, Wahba et al. (1998a) utilizaram seis exemplos de

torres estaiadas e compararam dois modelos computacionais distintos: no primeiro,

modelaram o mastro com elementos de treliça espacial e cabos com elementos de cabo não-

linear, enquanto que no segundo, utilizaram elementos de viga/coluna no mastro.

Compararam também os resultados obtidos com um modelo no qual a torre foi representada

como uma viga com suportes elásticos não-lineares. Os autores concluíram que o modelo em

treliça espacial não apresentou vantagem alguma sobre o de viga/coluna, sendo que no

segundo modelo, houve considerável redução no número de elementos e graus de liberdade.

Para o modelo de viga em suportes elásticos não-lineares, verificaram que este apresentou

deslocamentos menores, quando comparado com os modelos de elementos finitos.

Figura 1.4 – Ponte suspensa Tsing Ma.

Dando continuidade ao estudo, Wahba et al. (1998b) avaliaram a resposta

dinâmica dos seis exemplos analisados em Wahba et al. (1998a). As freqüências naturais

obtidas nos modelos concordaram entre si. Entretanto, a modelagem com elementos de

viga/coluna resultou em considerável ganho computacional. Com o propósito de verificar a

modelagem dos elementos finitos, dois módulos de aço triangulares na escala do modelo da

torre estaiada foram fabricadas e testadas em uma mesa vibratória. Os resultados

experimentais tiveram boa aproximação com os resultados das análises via elementos finitos.

Foi mostrado também que o acréscimo na tensão inicial dos cabos resultou em um acréscimo

significativo nas freqüências naturais.

5

Karoumi (1999) aplicou o elemento de catenária de dois nós para modelagem de

pontes suspensas ou estaiadas. Em particular, estudou-se o comportamento da ponte suspensa

de Great Belt (Fig. 1.5)1 construída entre Kastrup (na costa dinamarquesa) e Lernacken (no

litoral sueco), a qual possui vão suspenso de 1624 m e mastros com altura de 254 m. O autor

demonstrou que a modelagem dos cabos com vários elementos de catenária foi precisa para a

análise dinâmica, pois os elementos de catenária incluem o efeito da pré-tensão, além de

modelarem corretamente a mudança de geometria e de considerarem o peso próprio do cabo.

As análises da freqüência natural mostraram boa concordância dos resultados quando

comparados com dados publicados previamente.

Figura 1.5 – Ponte suspensa Great Belt.

Kahla (2000) investigou os efeitos da ruptura de cabos em torres estaiadas, sem

contudo, levar em consideração a ação do carregamento devido ao vento. Os resultados

revelaram que os cabos não excedem sua capacidade de suportar as solicitações de tração,

mas que um conjunto de membros do mastro falha na compressão, principalmente os

inferiores, o que pode levar ao colapso da estrutura.

Millar e Barghian (2000) apresentaram um estudo dos comportamentos estático e

dinâmico de estruturas flexíveis que exibem saltos dinâmicos. Um exemplo simples com

cabo, analisado usando os programas de elementos finitos UMIST e FINELE, foi suficiente

para mostrar que problemas não-lineares estáticos que exibem saltos dinâmicos podem ser 1 Fonte: <http://www.copenhagenpictures.dk/grt_blt.html>. Acesso em: 31 jul. 2007.

6

eficazmente analisados com uma abordagem dinâmica, pois as forças de amortecimento

podem ser anuladas na equação de movimento, ficando apenas os termos da inércia.

Preocupado com a limitada atenção dada ao comportamento sísmico de torres

estaiadas, Amari (2002) buscou propor alguns indicadores para ajudar os projetistas a

decidirem se efeitos sísmicos são importantes e se uma análise dinâmica detalhada da

estrutura é requerida. O estudo baseou-se nas análises sísmicas não-lineares de oito torres de

telecomunicação estaiadas com alturas variando entre 150 e 607 m. Resumindo os indicadores

de sensibilidade sísmica propostos, tem-se que torres com altura de 150 a 350 m podem

desenvolver significativo esforço cortante na base, na ordem de 40 a 80% do peso da torre.

Para torres acima de 350 m, não foi possível sugerir uma simples estimativa do esforço

cortante e, portanto, uma análise dinâmica detalhada é recomendada. Ainda para torres acima

de 350 m de altura sujeitas a acelerações laterais, uma análise dinâmica detalhada também é

recomendada. Por fim, independente da altura da torre, quando esta possui cabos sem tensão

inicial significativa uma análise dinâmica detalhada faz-se necessária.

Oliveira et al. (2002) realizaram uma análise comparativa entre o elemento finito

isoparamétrico de cabo com dois nós, desenvolvido a partir de uma formulação variacional, e

o elemento de catenária desenvolvido, a partir de expressões exatas oriundas da equação da

catenária elástica. Foi realizada uma análise não-linear verificando o comportamento estático

desses elementos quando submetidos a algumas condições de carregamento. Nos resultados

obtidos, foram observadas pequenas diferenças entre o elemento de catenária e os elementos

das formulações Lagrangeana Total e Atualizada em favor do elemento isoparamétrico.

Entretanto, devido à maior simplicidade no tratamento matemático de sua formulação, o

elemento de catenária requer menor esforço computacional. Além disso, o elemento de

catenária apresentou resultados mais consistentes para um carregamento concentrado.

Cheng et al. (2002), baseado no conceito de instabilidade do ponto limite, propôs

um elemento finito não-linear avançado para analisar a estabilidade aerostática de pontes

estaiadas. Foi utilizado um procedimento incremental para determinar com precisão a

velocidade crítica do vento, bem como para reduzir consideravelmente o tempo

computacional no cálculo. Na modelagem de cada cabo, foi usado um único elemento de

catenária com dois nós. O algoritmo para formar a matriz de rigidez tangente do elemento foi

descrito por Karoumi (1999). O estudo de caso da segunda ponte estaiada de Santou Bay na

China, validou a eficiência do método. Os autores concluíram que a resposta do deslocamento

durante aplicação das cargas de vento demonstrou forte não-linearidade.

7

Chan et al. (2002) apresentaram uma análise de segunda ordem da estabilidade de

colunas estaiadas protendidas com imperfeições iniciais. Um elemento do tipo pointwise

equilibrium polynomial (PEP) foi usado em associação com o elemento de cabo proposto por

Wei et al. (1999). Exemplos numéricos demonstraram a versatilidade do procedimento

numérico proposto e foi apresentado um estudo paramétrico para avaliar a influência do

carregamento, da imperfeição inicial, do tipo de apoio e do diâmetro do cabo no modo de

flambagem. Dos resultados, observou-se que a capacidade de flambagem de uma coluna pode

ser significativamente acrescida pelo uso de apoios e cabos pré-tensionados, o elemento PEP

modela bem as imperfeições iniciais e o procedimento incremental interativo pôde determinar

corretamente o comportamento não-linear de uma imperfeição na coluna.

Yan-Li et al. (2003) aplicaram o método discreto de vibrações na análise da

resposta de uma torre estaiada sob ventos induzidos. Quando considerados cabos não-lineares,

foi associada à teoria de truncamento de Gauss. O método discreto apresentou boa

convergência e boa estabilidade, sendo as respostas próximas às obtidas com um modelo

ensaiado em túnel de vento, atestando assim a confiabilidade do método para torres estaiadas.

Harikrishna et al. (2003) descreveram com maior profundidade as características do vento,

medindo em escala real a resposta estrutural de uma torre estaiada com 50 m de altura em

condições de vento ambiente.

Para estudar a estabilidade estática e dinâmica, Pasquetti (2003), com base no

princípio da energia potencial mínima, apresentou uma formulação discreta para um modelo

plano de torre estaiada, onde os cabos foram modelados com elementos de mola linear,

elementos de mola não-linear e como cabos inextensíveis. Na análise estática foi determinada

a carga crítica e o caminho pós-crítico, associados ao estudo paramétrico dos diversos

parâmetros físicos e geométricos determinantes para a estabilidade da torre. Na análise

dinâmica foi realizado o estudo paramétrico da freqüência natural. O estudo mostrou que a

torre pode apresentar diversos comportamentos típicos de sistemas não-lineares, tais como

saltos, bifurcações de período e caos.

Campos Filho (2004) apresentou uma formulação via método dos elementos

finitos (MEF) não-lineares, com referencial Lagrangiano atualizado, para modelar o

comportamento estático de estruturas estaiadas. Na modelagem dos elementos de pórtico, foi

utilizado o elemento de viga e na dos cabos o elemento de treliça. Foi aplicado o método de

Newton-Raphson, associado ao método do comprimento de arco ou método dos

deslocamentos generalizados, para determinar as configurações de equilíbrio da estrutura. Os

8

resultados mostraram que não houve diferenças significativas nos resultados das trajetórias

obtidas com as duas técnicas numéricas. Foi realizada, ainda, uma análise paramétrica para o

estudo da influência do posicionamento, do tensionamento e da rigidez dos cabos. Estes

fatores mostraram ser determinantes no comportamento estrutural, havendo um forte

acréscimo no valor da carga crítica e perda de rigidez do sistema após a flambagem.

No tocante aos mecanismos de controle de vibrações, Yau e Yang (2004)

estudaram a redução da vibração de pontes estaiadas sujeitas à passagem de trens de alta

velocidade com o emprego de amortecedores de massa sintonizados. O trem foi modelado

com uma série de massas-molas, o tabuleiro e as torres da ponte por elementos não-lineares

de viga/coluna e os cabos por elementos de treliça. Os exemplos numéricos utilizados no

estudo demonstraram que o sistema proposto supre com eficiência os picos de ressonância.

Law et al. (2005) desenvolveram um método para indicar a variação no tempo da

carga de vento na resposta estrutural. Os autores utilizaram uma torre estaiada de 50 m de

altura para validar a eficiência do método. As simulações numéricas mostraram que o método

proposto representa com precisão a resposta estrutural para carga de vento. Contudo, em

algumas direções, a carga de vento medida apresentou ruído no espectro da freqüência.

Freire et al. (2006) avaliaram os efeitos não-lineares na análise estática de uma

ponte de aço estaiada. Foram usadas uma análise linear e uma análise pseudo-linear baseada

na modificação dos módulos de elasticidades dos elementos de treliça, bem como uma análise

não-linear com vários modelos numéricos de elementos finitos. Os resultados mostraram que

a curvatura dos cabos origina os efeitos não-lineares mais importantes, podendo ser decisivos

no comportamento global dessas estruturas, especialmente quando grandes deslocamentos

levam ao surgimento de tensões nos cabos. Os resultados mostraram também que os

deslocamentos obtidos pela análise não-linear, usando elementos de catenária elásticos,

podem ser o dobro daqueles obtidos usando cabos modelados com elementos de treliça.

Segundo os autores, com o aumento da deflexão do tabuleiro, o efeito da curvatura do cabo na

resposta não-linear diminuiu e os elementos de treliça puderam ser usados na modelagem dos

cabos. Devido à grande não-linearidade de pontes estaiadas com grandes vãos, a análise linear

não forneceu resultados satisfatórios. Quanto à modificação dos módulos de elasticidade, os

autores recomendaram cautela na utilização, pois o método provou ser limitado ou impróprio.

Andreu et al. (2006) apresentaram uma formulação de um elemento de catenária

que permite a simulação de redes compostas por múltiplos cabos. A formulação proposta,

originada de uma modificação das equações convencionais para cabos inextensíveis,

9

assegurou o equilíbrio depois da deformação do cabo. O uso de expressões analíticas para

descrição geométrica do elemento e a correspondente matriz de rigidez tangente forneceram

alta estabilidade numérica. Em comparação com resultados numéricos e experimentais

encontrados por diferentes pesquisadores, o elemento foi eficiente e preciso.

Orlando (2006) estudou o comportamento não-linear de um sistema torre-pêndulo

submetido a um carregamento harmônico, onde abordou aspectos gerais ligados à estabilidade

dinâmica e ao controle de vibrações. Foi apresentada uma formulação não-linear para obter o

funcional de energia, bem como as freqüências naturais e os modos de vibração do sistema.

Para descrever o comportamento na vizinhança da freqüência fundamental e obter as

equações de movimento, foi adotado, com base na análise modal do sistema, um modelo com

dois graus de liberdade. Uma análise paramétrica das oscilações não-lineares demonstrou que

o absorsor pendular passivo pode reduzir ou amplificar a resposta da coluna. No estudo da

influência da não-linearidade geométrica do pêndulo ficou evidente que a não-linearidade não

pode ser desprezada nessa classe de problema. Por fim, com base nos resultados, foi proposto

um absorsor pendular híbrido, por ser um controle mais eficiente que o passivo e por não

requerer grande gasto de energia.

1.3 OBJETIVOS

Essa dissertação está inserida na linha de pesquisa em Instabilidade e Dinâmica de

Estruturas da Escola de Engenharia Civil da Universidade Federal de Goiás. Com esta

pesquisa, pretende-se aprimorar as ferramentas para análise estática e dinâmica de estruturas

estaiadas planas existentes na escola.

O primeiro objetivo deste trabalho é aprofundar o estudo do comportamento

estático de estruturas estaiadas planas, avaliando os efeitos das imperfeições iniciais no

regime não-linear e estudar a instabilidade dinâmica de torres estaiadas carregadas axialmente

tanto por cargas súbitas quanto harmônicas.

Em segundo lugar, usando um modelo discreto desenvolvido para um sistema de

torre estaiada por molas, tem-se como objetivo examinar o comportamento não-linear estático

e dinâmico, bem como a instabilidade dinâmica destas estruturas, analisando os diversos tipos

de bifurcação e bacias de atração associados às fronteiras de estabilidade.

10

1.4 ESCOPO

No capítulo dois descreve-se, com base no funcional de energia, a formulação do

elemento finito de pórtico plano e de cabo, utilizados na análise não-linear geométrica das

estruturas, deduzindo as matrizes de rigidez e de massa desses elementos. Apresenta-se

também, uma formulação não-linear estática e dinâmica do modelo discreto torre-cabo, onde

os cabos são modelados por molas lineares. Destacam-se nesse capítulo uma descrição

resumida dos sistemas de referência utilizados, das relações deformação-deslocamento, das

relações constitutivas, da determinação das freqüências naturais e da consideração do

amortecimento.

No capítulo três apresentam-se, de forma sucinta, alguns conceitos fundamentais

da teoria da estabilidade estrutural e os métodos numéricos usados para calcular as respostas

estática e dinâmica das estruturas estaiadas planas. Para a análise foram empregadas diversas

ferramentas numéricas, a saber: iterações do tipo Newton-Raphson associadas à técnica do

comprimento de arco para caracterizar a solução não-linear, integração das equações

diferenciais usando o método de Runge-Kutta, método de Newmark para obter a resposta no

tempo na análise dinâmica, algoritmo para o cálculo das fronteiras de instabilidade e cálculo

dos diagramas de bifurcação, bem como das bacias de atração usando o método da força

bruta. Neste trabalho, as implementações computacionais para o MEF foram realizadas

através da linguagem de programação C++ e desenvolvimento do modelo discreto foi

realizado usando-se o programa de álgebra simbólica MAPLE.

No capítulo quatro apresentam-se os exemplos numéricos utilizados para validar

as formulações matemáticas e para avaliar os caminhos pós-críticos, as freqüências naturais,

os modos de vibração de colunas perfeitas e imperfeitas e a instabilidade dinâmica sob carga

súbita e harmônica. Descrevem-se as oscilações não-lineares, os mecanismos de escape e a

evolução da estabilidade global através das bacias de atração de alguns exemplos.

Finalmente, no capítulo cinco, são apresentadas as conclusões obtidas no estudo,

bem como algumas sugestões para trabalhos futuros.

11

2 FORMULAÇÃO MATEMÁTICA

2.1 INTRODUÇÃO

O comportamento de uma estrutura submetida a um conjunto de ações externas

pode ser definido como a relação existente entre a intensidade, a combinação e o histórico de

aplicação dessas ações e os efeitos provocados por elas na estrutura em termos de tensões,

deformações, deslocamentos, etc. A determinação desse comportamento requer a

consideração conjunta das equações descritas a seguir:

a) Equações de Equilíbrio: consideram as forças aplicadas, os esforços e as

tensões;

b) Relações Constitutivas: descrevem o comportamento dos materiais que

constituem a estrutura, relacionando esforços ou tensões e deslocamentos ou

deformações;

c) Relações Cinemáticas: consideram as relações entre deformações e

deslocamentos;

d) Equações de Compatibilidade: envolvem deslocamentos ou deformações e são

destinadas a garantir que a estrutura respeite as ligações dos vários elementos

entre si e os vínculos externos.

A análise da maioria das estruturas está baseada na premissa de que as hipóteses

de comportamento físico e geométrico linear atendem de forma satisfatória a determinação do

comportamento das estruturas. No entanto, em alguns casos, esta análise fornece uma

avaliação incorreta com relação ao comportamento estrutural e uma análise não-linear faz-se

necessária por permitir avaliar melhor a estrutura e os materiais que a compõem.

12

No presente trabalho, uma análise geométrica não-linear será realizada e a teoria a

ser apresentada será baseada nos estudos realizados por Silveira (1995), Galvão (2000),

Oliveira (2002) e Campos Filho (2004), que implementaram formulações geometricamente

não-lineares para análise de sistemas estruturais planos formados por barras ou por cabos.

Entre as formulações implementadas, será utilizada a formulação proposta por Yang e Kuo

(1994), a qual se baseia nas relações de deformação-deslocamento geradas a partir do tensor

de Green-Lagrange na sua forma completa. Serão utilizados elementos de viga/coluna para

representação dos elementos componentes do pórtico e elementos de cabo/treliça para a

representação dos cabos.

2.2 SISTEMAS DE REFERÊNCIAS

Ao se considerar a não-linearidade geométrica, as equações de equilíbrio são

definidas na configuração deformada e a hipótese dos pequenos deslocamentos deixa de ser

válida, tendo em vista a ocorrência de grandes deslocamentos e/ou de grandes deformações.

Desta forma, deve ser adotada uma formulação incremental iterativa para que se conheça o

comportamento da estrutura à medida que se varia o carregamento.

Como as variáveis descrevem as alterações no comportamento da estrutura ao

longo do processo de deformação, na análise não-linear não se considera um sistema de

coordenadas estacionário. Assim, grande parte das formulações pelo MEF com não-

linearidade geométrica encontradas na literatura baseia-se em referenciais Lagrangianos.

Nesses referenciais três tipos de configurações podem ser concebidos para medir os

deslocamentos decorrentes de um dado carregamento: a configuração inicial, a última

configuração deformada t e a configuração deformada corrente tt ∆+ . Por hipótese, assume-

se que todas as variáveis de estado, tais como tensões, deformações e deslocamentos,

juntamente com a história de carregamento, são conhecidas na configuração t. A partir daí, a

questão principal passa a ser a formulação de um processo incremental para determinar todas

essas variáveis de estado para o corpo na configuração tt ∆+ , considerando que o

carregamento externo atuando na configuração t tenha sofrido um pequeno acréscimo de

valor. O passo que caracteriza o processo de deformação do corpo de t para tt ∆+ é

tipicamente referido como um passo incremental.

Dependendo da configuração anterior selecionada como referência para a

obtenção do estado de equilíbrio do corpo na configuração deformada corrente, tt ∆+ , dois

tipos de referenciais Lagrangianos podem ser identificados: o referencial Lagrangiano

13

atualizado (RLA), onde a última configuração t de equilíbrio é selecionada como o estado de

referência, e o referencial Lagrangiano total (RLT) que utiliza a configuração inicial

indeformada para o mesmo propósito.

Para o RLT os deslocamentos são medidos em relação à configuração inicial

indeformada da estrutura, como mostra o esquema da Fig. 2.1. No RLA os deslocamentos são

medidos em relação à última configuração de equilíbrio obtida no processo incremental, ou

seja, em relação a um referencial que é atualizado a cada incremento de carga, conforme a

Fig. 2.2, onde YGL e XGL são os eixos cartesianos das coordenadas globais, 0Y e 0X são os

eixos cartesianos das coordenadas na configuração inicial 0=t , tY e tX são os eixos

cartesianos das coordenadas na configuração deformada corrente t e, finalmente, u1, v1, u2 e

v2 são as componentes de deslocamento nodais.

Figura 2.1 – Referencial Lagrangiano Total.

Figura 2.2 – Referencial Lagrangiano Atualizado.

14

Por terem características distintas, diferentes tensores de tensão e deformação

devem ser utilizados para cada formulação. Contudo, se forem utilizadas relações tensão-

deformação apropriadas, os dois procedimentos levarão a resultados equivalentes. Para o

presente trabalho, foi adotado o RLA.

Uma descrição mais aprofundada sobre as formas de referencial Lagrangiano

pode ser encontrada em Bathe (1996) e Crisfield (1991).

2.3 RELAÇÃO DEFORMAÇÃO-DESLOCAMENTO

A formulação proposta por Yang e Kuo (1994) baseia-se nas relações de

deformação-deslocamento, geradas a partir do tensor de Green-Lagrange na sua forma

completa. Consideram-se, além das relações correspondentes às deformações axiais, as

correspondentes às deformações de cisalhamento dadas por:

,21

21

,21 22

⎥⎦

⎤⎢⎣

⎡ ∆∆+

∆∆+⎥

⎦

⎤⎢⎣

⎡ ∆+

∆=∆

⎥⎥⎦

⎤

⎢⎢⎣

⎡⎟⎠⎞

⎜⎝⎛ ∆

+⎟⎟⎠

⎞⎜⎜⎝

⎛ ∆+

∆=∆

dxvd

dyvd

dxud

dyud

dxvd

dyud

dxvd

dxud

dxud

XY

XX

ε

ε

(2.1)

onde u∆ é o deslocamento axial de um ponto distante y da linha neutra da seção, v∆ é o

deslocamento vertical, XXε∆ é o incremento da deformação axial e XYε∆ é o incremento da

deformação cisalhante. Separando-se as parcelas lineares das não-lineares da Eq. 2.1 tem-se:

,,

xyxyXY

xxxxXX

ee

ηεηε

∆+∆=∆∆+∆=∆

(2.2)

onde as componentes de cada parcela são definidas como:

.21,

21

,21,

22

⎥⎦

⎤⎢⎣

⎡ ∆∆+

∆∆=∆⎥

⎦

⎤⎢⎣

⎡ ∆+

∆=∆

⎥⎥⎦

⎤

⎢⎢⎣

⎡⎟⎠⎞

⎜⎝⎛ ∆

+⎟⎟⎠

⎞⎜⎜⎝

⎛ ∆=∆

∆=∆

dxvd

dxvd

dxud

dyud

dxvd

dyude

dxvd

dxud

dxude

xyxy

xxxx

η

η

(2.3)

Admitindo-se a teoria de flexão de Euler-Bernoulli, a qual considera que as seções

transversais perpendiculares ao eixo da barra antes da flexão continuem planas, indeformadas

e perpendiculares ao eixo, após a flexão, tem-se:

15

dxvdyuu ∆

−∆=∆ , (2.4)

onde a primeira parcela é uma conseqüência dos esforços normais atuantes, sendo constante

ao longo da seção; a segunda parcela é devida aos esforços de flexão que, por sua vez, varia

linearmente com a distância à linha neutra (LN), como pode ser visto na Fig. 2.3.

A C

B D∆xx

L.N.

(a)

M

O

L.N.

A'

B'

C'

D'

M

(b)

L.N.

∆u

∆u

y

dxd∆v

xdv∆d

A' C'

BD

(c)

Figura 2.3 – Deformação na flexão: (a) Elemento indeformado; (b) Elemento deformado; (c)

Elemento infinitesimal.

Substituindo-se a Eq. 2.4 na Eq. 2.3 e desprezando-se os termos não-lineares

relativos à deformação por cisalhamento, XYε∆ , obtêm-se:

.21

,021

,221

,

2

2

22

2

22

2

22

2

2

⎥⎦

⎤⎢⎣

⎡⎟⎟⎠

⎞⎜⎜⎝

⎛ ∆∆+

∆∆−=∆

=⎥⎦⎤

⎢⎣⎡ ∆

+∆

−=∆

⎥⎥⎦

⎤

⎢⎢⎣

⎡⎟⎠⎞

⎜⎝⎛ ∆

+⎟⎟⎠

⎞⎜⎜⎝

⎛ ∆+

∆∆−⎟

⎠⎞

⎜⎝⎛ ∆

=∆

∆−

∆=∆

dxvd

dxvdy

dxvd

dxud

dxvd

dxvde

dxvd

dxvdy

dxvd

dxudy

dxud

dxvdy

dxude

xy

xy

xx

xx

η

η

(2.5)

16

2.4 RELAÇÕES CONSTITUTIVAS

Considerando-se o estado de tensões ou de deformações na última configuração

de equilíbrio obtida no processo de solução incremental, ou seja, na configuração t, os

esforços resultantes podem ser obtidos por:

,

,

,

∫

∫

∫

=

=

=

Axx

tt

Axy

tt

Axx

tt

dAyM

dAQ

dAP

τ

τ

τ

(2.6)

onde Pt é a força axial, Qt é a força cisalhante e Mt é o momento fletor na seção

transversal.

A partir das forças de extremidades representadas na Fig. 2.4, pode-se obter os

esforços solicitantes ao longo do elemento pelas relações:

Figura 2.4 – Esforços no elemento.

( ) ( ),; 2121

1 LMMQx

LMMMM tt +

−=+

+−= (2.7)

sendo Mt e Qt o momento fletor e o esforço cortante, respectivamente, do elemento de barra

em uma seção localizada a uma distância x do nó inicial.

Segundo Yang e Kuo (1994), para um elemento no qual a hipótese de Euler-

Bernoulli tenha sido adotada, apenas as tensões axiais incrementais, xxtτ∆ , podem ser obtidas

diretamente da lei constitutiva dada:

xxxxt E ετ ∆=∆ . (2.8)

17

As tensões cisalhantes incrementais xytτ∆ devem ser determinadas a partir da

condição de equilíbrio:

0=∂

∂+

∂∂ ∆∆

y

xyt

x

xxt ττ

. (2.9)

Como yIM

xxt ⎟

⎠⎞

⎜⎝⎛ ∆

=∆ τ , a Eq. 2.9 pode ser reescrita na forma:

∫∆

−=∆ dyIy

dxMd

xytτ , (2.10)

onde y é a distância em relação à linha neutra e I o momento de inércia da seção.

Para um elemento com seção retangular de dimensões b e h representando,

respectivamente, a base e a altura do elemento, a aplicação da condição de contorno 0=∆xy

tτ

para 2hy ±= conduz a:

⎟⎟⎠

⎞⎜⎜⎝

⎛−⎟

⎠⎞

⎜⎝⎛∆

−=∆

421 2

2 hyIdx

Mdxy

tτ . (2.11)

Como dAyMA

xx∫= τ obtém-se:

3

322

3

322

4242 dxvd

dxudhyE

dxvdhyE

xyt ∆∆

⎟⎟⎠

⎞⎜⎜⎝

⎛−+

∆⎟⎟⎠

⎞⎜⎜⎝

⎛−=∆ τ , (2.12)

ou seja,

nxyxyxy

t SS +=∆ 1τ , (2.13)

sendo:

3

322

3

3221

42

42

dxvd

dxudhyES

dxvdhyES

nxy

xy

∆∆⎟⎟⎠

⎞⎜⎜⎝

⎛−=

∆⎟⎟⎠

⎞⎜⎜⎝

⎛−=

, (2.14)

onde 1xyS é a parcela linear e n

xyS é a parcela não-linear das tensões cisalhantes incrementais, E

é o módulo de elasticidade longitudinal do material e y é a distância à linha neutra.

18

2.5 FUNCIONAL DE ENERGIA

Em um corpo tridimensional que sofre grandes deslocamentos quando submetido

à ação de forças externas, a avaliação da posição final do equilíbrio é feita mediante a análise

das sucessivas deformações sofridas no processo de carregamento, até que se atinja a

configuração de equilíbrio entre o trabalho das forças internas e externas. A Fig. 2.5

representa os deslocamentos deste corpo e suas respectivas configurações de equilíbrio.

A energia potencial total na configuração tt ∆+ pode ser calculada como:

( ) ( )WUWUWWUUWU ttttttttt ∆+∆++=∆++∆+=+=Π ∆+∆+∆+ , (2.15)

onde Π é a energia potencial total, U é a energia interna de deformação e sua respectiva

variação U∆ e W é o trabalho das forças externas e sua respectiva variação W∆ .

Figura 2.5 – Movimento de um corpo tridimensional.

Da primeira variação da energia potencial total dada na Eq. 2.15 obtém-se:

( ) ( )WUWUttt ∆+∆++=Π∆+ δδδ . (2.16)

Como o corpo encontra-se em equilíbrio na configuração t, tem-se que

( ) 0=+ WUtδ . Assim:

( ) 0=∆+∆=Π∆+ WUtt δδ , (2.17)

com isso, o acréscimo no funcional de energia na configuração tt ∆+ é definido como:

WU ∆+∆=∆Π , (2.18)

assim, a energia interna de deformação, U∆ , pode ser escrita como:

19

dVdUV

ijij

ijijt

ijt

∫ ∫ ⎥⎥⎦

⎤

⎢⎢⎣

⎡=∆

∆+ εε

ε

ετ . (2.19)

Expandindo a Eq. 2.19 em função das tensões e das deformações axiais e

cisalhantes definidas nas Eqs. 2.5, 2.8 e 2.13, chega-se a:

( ) ∫∫ ⎟⎠⎞

⎜⎝⎛ ∆+∆+∆+∆=∆ ∆

vxyxy

txx

vxxxx

txxxx

t dVEdVU ετεετετ 22

2 2 . (2.20)

O trabalho realizado pelas forças externas nodais, W∆ , é definido como:

⎥⎦

⎤⎢⎣

⎡+−=−=∆ ∫∫∫ ∆+

AA

t

A

tt dAdAdAW ∆u∆Fi∆uFi∆uFi , (2.21)

onde iFt é o vetor de forças internas no instante t, ∆u é o vetor de deslocamentos nodais

incrementais e i∆F é o vetor de incremento das forças internas.

ALVES (1995) propõe que o funcional de energia seja escrito como:

⎥⎦

⎤⎢⎣

⎡+−++++=∆Π ∫∫

AA

tL dAdAUUUUU ∆u∆Fi∆uFi210 τ , (2.22)

onde 0U representa a parcela de U∆ que corresponde às forças acumuladas até a

configuração de equilíbrio t, τU representa a parcela de U∆ que corresponde à influência das

deformações iniciais e que dá origem à matriz de rigidez geométrica, LU representa a parcela

de U∆ que dá origem à matriz de rigidez linear LK e, finalmente, 1U e 2U representam as

parcelas de U∆ que resultarão nas matrizes de rigidezes não-lineares, sendo:

( )

.22

,2

,2

,21

,2

22

11

2

0

dVSEU

dVSeEU

dVU

dVeU

dVeeU

Vxy

nxyxx

Vxyxyxxxx

xyxyt

Vxxxx

t

V

xxL

Vxyxy

txxxx

t

∫

∫

∫

∫

∫

⎟⎠⎞

⎜⎝⎛ ∆+∆=

∆+∆∆=

∆+∆=

∆=

∆+∆=

ηη

ηη

ητητ

ττ

τ (2.23)

20

Com o objetivo de obter uma equação na forma incremental, ALVES (1995)

efetuou uma simplificação baseada no princípio dos trabalhos virtuais da forma:

( )dVeedAV

xyxyt

xxxxt

A

t ∫∫ ∆+∆= ττ 2∆uFi , (2.24)

ou então:

∫=A

t dAU ∆uFi0 , (2.25)

onde o lado direito da Eq. 2.25 corresponde ao trabalho virtual das forças externas acumulada

até a configuração t. Assim, o funcional de energia dado na Eq. 2.22 pode ser reescrito na

forma:

∫−+++=∆ΠA

L dAUUUU ∆u∆Fi21τ . (2.26)

Substituindo-se convenientemente as Eqs. 2.5, 2.6, 2.12 na Eq. 2.23 e, após

simplificações, as parcelas da energia interna de deformação podem ser definidas como:

.2

48

,

23

2

,21

21

,

21

3

3

2

22

2

2

2

2

2222

0

22

2

3

3

2

2

2

0

22

1

0

2

2

2

0

2

002

2

0

2

2

222

dxdx

vddx

uddx

udEIdx

vddx

udEI

dxvd

dxvd

dxudEI

dxvd

dxudEAU

dxdx

vddx

vddx

udEI

dxvd

dxudEI

dxud

dxvd

dxudEAU

dxdx

vdEIdxdx

udEAU

dxdx

vddx

udQdxdx

vddx

udM

dxdx

vdAI

dxvd

dxudPU

L

L

LL

L

LL

L

∆∆∆+⎟⎟

⎠

⎞⎜⎜⎝

⎛ ∆∆+

+⎟⎟⎠

⎞⎜⎜⎝

⎛ ∆⎟⎟⎠

⎞⎜⎜⎝

⎛ ∆+

∆+⎟⎟

⎠

⎞⎜⎜⎝

⎛ ∆+

∆=

∆∆∆+

+⎟⎟⎠

⎞⎜⎜⎝

⎛ ∆∆+

∆⎟⎟⎠

⎞⎜⎜⎝

⎛ ∆+

∆=

⎟⎟⎠

⎞⎜⎜⎝

⎛ ∆+

∆=

∆∆−

∆∆−

−⎥⎥⎦

⎤

⎢⎢⎣

⎡⎟⎟⎠

⎞⎜⎜⎝

⎛ ∆+

∆+

∆=

∫

∫

∫∫

∫∫

∫τ

(2.27)

21

2.6 FORMULAÇÃO DO ELEMENTO FINITO

2.6.1 FORMULAÇÃO DO ELEMENTO DE PÓRTICO PLANO

Para a análise dos componentes do pórtico, optou-se por um elemento finito plano

do tipo viga-coluna, conforme mostra a Fig. 2.6. Este elemento é esquematizado como um

segmento reto, limitado pelos nós 1 e 2, que se deforma no plano de definição da estrutura.

Para cada elemento define-se um sistema local de coordenadas x e y, que faz um ângulo α em

relação ao sistema de coordenadas global XGL e YGL. Os deslocamentos nodais no sistema de

coordenadas local ∆u1, ∆v1, e ∆θ1 são, respectivamente, o deslocamento axial, o deslocamento

transversal e a rotação do nó 1 e ∆u2, ∆v2, e ∆θ2 são, respectivamente, o deslocamento axial, o

deslocamento transversal e a rotação do nó 2.

Figura 2.6 – Elemento de Pórtico.

Pela hipótese básica do MEF, para que haja convergência dos resultados, é

necessário que as funções aproximadoras escolhidas para representar os deslocamentos nodais

sejam contínuas no domínio do elemento finito. Ao observar os termos da Eq. 2.27, nota-se a

necessidade de se realizar integrações no domínio do elemento por meio de operações:

∫L

n

n

dxdx

vd

0

, (2.28)

onde v representa genericamente os deslocamentos.

Pelo cálculo diferencial, a Eq. 2.28 é definida somente se v for contínua até a

ordem (n-1), ou seja, v e as derivadas de v até a ordem (n-1) forem contínuas. Assim, para que

haja continuidade de deslocamentos, representados pelas translações e pela rotação nas bordas

22

dos elementos adjacentes, é suficiente aproximar o deslocamento axial, u∆ , por uma função

linear, enquanto que, para o deslocamento transversal, v∆ , admitindo-se dxv∆=∆θ , deve-

se adotar uma função do terceiro grau. Dessa forma:

,

,3

32

210

10

xbxbxbbv

xaau

+++=∆

+=∆ (2.29)

onde a0, a1, b0, b1, b2 e b3 são coeficientes determinados por meio da aplicação das condições

de contorno do elemento em x=0 (∆u=∆u1, ∆v=∆v1, ∆θ1=∆v/dx) e em x=L (∆u=∆u2, ∆v=∆v2,

∆θ2=∆v/dx). Substituindo-se os coeficientes a0, a1, b0, b1, b2 e b3 da Eq. 2.29 e rearranjando-se

convenientemente os termos das equações, chega-se a:

,1,

2625413

2211

θθ ∆+∆+∆+∆=∆∆+∆=∆

HvHHvHvuHuHu

(2.30)

onde H1, H2, H3, H4, H5 e H6 são funções de interpolação dadas por:

.,23,2

,231,,1

2

33

63

3

2

2

52

3

4

3

3

2

2

321

Lx

LxH

Lx

LxH

Lx

LxxH

Lx

LxH

LxH

LxH

+−=−=+−=

+−==−= (2.31)

Assim, as translações u∆ e v∆ e a rotação θ∆ de um ponto do elemento situado

a uma distância x do nó 1 podem ser escritas na forma matricial apresentada:

∆uH∆d = , (2.32)

onde [ ]Tvu ∆∆=∆d , [ ]Tvuvu 222111 θθ ∆∆∆∆∆∆=∆u e H tem a forma:

⎥⎦

⎤⎢⎣

⎡=⎥

⎦

⎤⎢⎣

⎡=

6543

21

000000

HHHHHH

NN

v

uH . (2.33)

A transformação de ∆u para o sistema de coordenadas global XGL e YGL, que é o

referencial comum a todos os elementos, é dada por:

∆uR∆u TGL = , (2.34)

sendo R a matriz de rotação do elemento:

23

⎥⎥⎥⎥⎥⎥⎥⎥

⎦

⎤

⎢⎢⎢⎢⎢⎢⎢⎢

⎣

⎡

−

−

=

1000000cossen0000sencos0000001000000cossen0000sencos

αααα

αααα

R . (2.35)

2.6.2 FORMULAÇÃO DO ELEMENTO DE CABO

Para a análise dos componentes do cabo, optou-se pelo elemento de treliça,

conforme mostra a Fig. 2.7. Este elemento é esquematizado como um segmento reto, limitado

pelos nós 1 e 2. Da mesma forma que o elemento de pórtico, para cada elemento de treliça

define-se um sistema local de coordenadas x e y, que faz um ângulo α em relação ao sistema

de coordenadas global XGL e YGL. Os deslocamentos nodais no sistema de coordenadas local

1u∆ e 1v∆ são, respectivamente, o deslocamento axial e o deslocamento transversal do nó 1, e

2u∆ e 2v∆ são, respectivamente, o deslocamento axial e o deslocamento transversal do nó 2.

Como o elemento de treliça só apresenta deformações axiais, ele é comumente

representado com deslocamentos u∆ apenas. No entanto, os deslocamentos v∆ são

necessários para representar deslocamentos de corpo rígido do elemento. Assim, adota-se uma

função linear para representar o deslocamento u∆ :

xaau 10 +=∆ , (2.36)

onde 0a e 1a são coeficientes a serem determinados por meio da aplicação das condições de

contorno do elemento em ( )10 uux ∆=∆= e em ( )2uuLx ∆=∆= .

Figura 2.7 – Elemento de treliça.

24

Substituindo-se os coeficientes 0a e 1a encontrados na Eq. 2.36 e rearranjando-se

convenientemente os termos das equações, pode-se escrever:

2211 uHuHu ∆+∆=∆ , (2.37)

onde H1 e H2 são funções de interpolação dadas por:

.

,1

2

1

LxHLxH

=

−= (2.38)

As translações v∆ que permitem o deslocamento de corpo rígido serão

interpoladas pelas mesmas funções de interpolação utilizadas para o deslocamento u∆ , pois a

barra deve permanecer retilínea após a deformação. Assim:

2211 vHvHv ∆+∆=∆ . (2.39)

Dessa forma, os deslocamentos em uma seção situada a uma distância x do nó 1

podem ser obtidos como:

∆uH∆d = , (2.40)

onde [ ]Tvu ∆∆=∆d , [ ]Tvuvu 2211 ∆∆∆∆=∆u e H tem a forma:

⎥⎦

⎤⎢⎣

⎡=⎥

⎦

⎤⎢⎣

⎡=

43

21

0000

HHHH

NN

v

uH . (2.41)

A transformação de ∆u para o sistema de coordenadas global XGL e YGL, que é o

referencial comum a todos os elementos, é dada por:

∆uR∆u TGL = , (2.42)

sendo R a matriz de rotação do elemento:

⎥⎥⎥⎥

⎦

⎤

⎢⎢⎢⎢

⎣

⎡

−

−=

αααα

αααα

cossen00sencos00

00cossen00sencos

R . (2.43)

25

2.7 MATRIZ DE RIGIDEZ DOS ELEMENTOS

Seja uma expansão em série de Taylor do funcional de energia potencial total,

∆Π , para um grau de liberdade em função dos deslocamentos nodais u∆ :

44

43

3

32

2

2

0 241

61

21 u

uu

uu

uu

un ∆∆∂∆Π∂

+∆∆∂∆Π∂

+∆∆∂∆Π∂

+∆∆∂∆Π∂

+∆Π=∆Π δδδδ , (2.44)

considerando-se que o elemento encontra-se em equilíbrio para uma determinada

configuração t em u∆ , então 00 =∆Π . Assim, a Eq. 2.44 pode ser escrita como:

uu

uuu

uuu

un ∆∆∂∆Π∂

+∆⎥⎦

⎤⎢⎣

⎡∂

∆∂∆Π∂

+∂∆∂∆Π∂

+∆∂∆Π∂

∆=∆Π δδδ 24

4

3

3

2

2

241

61

21 . (2.45)

A parcela da Eq. 2.45 relativa à derivada de 1ª ordem do funcional de energia

potencial total, ∆Π , representa o termo linear referente ao trabalho da força residual,

provenientes da diferença entre a força interna e a força externa:

FrFigu

λ−==∆∂∆Π∂ , (2.46)

onde g é a força residual, λ é o parâmetro de carga, Fr é a carga de referência e Fi é a força

interna.

Para os casos com mais de um grau de liberdade o funcional de energia potencial

total, ∆Π , pode ser escrito na forma da Eq. 2.47, apresentada por ALVES (1995):

( ) ( )∆uFrFi∆uKKKK∆u λτ −+⎥⎦⎤

⎢⎣⎡ +++=∆Π 21 24

161

21

LT

n , (2.47)

sendo LK , τK , 1K e 2K as parcelas da matriz de rigidez da estrutura, obtidas diretamente

das expressões da energia interna de deformação definidas pelos termos dados na Eq. 2.27. A

matriz de rigidez LK corresponde à parcela linear da matriz de rigidez, a matriz τK

corresponde à parcela relativa às tensões iniciais e as matrizes 1K e 2Κ correspondem às

parcelas não-lineares dependentes dos deslocamentos. Os elementos que compõem as

referidas matrizes devem ser calculados como:

26

.

,

,

,

11

14

2

13

1

2

2

uuuuuu

UK

uuuu

UK

uuU

K

uuU

K

kkji

ij

kkji

ij

jiij

ji

LijL

∆∂∆∆∂∆∂∆∂∆∂

∂=

∆∆∂∆∂∆∂

∂=

∆∂∆∂∂

=

∆∂∆∂∂

=

ττ

(2.48)

Do princípio da energia potencial total estacionária, tem-se que a condição de

equilíbrio do sistema em uma configuração t qualquer é:

FrFi∆uKKKΚ λτ =+⎥⎦⎤

⎢⎣⎡ +++ 21 6

121

L , (2.49)

ou representada genericamente na forma:

FrFi λ=)(u , (2.50)

onde Fi∆FiFi +=)(u e o incremento das forças internas é definido como:

∆uKKKK∆Fi ⎥⎦⎤

⎢⎣⎡ +++= 21 6

121

τL . (2.51)

Segundo Yang e Kuo (1994), não é necessário utilizar os termos de ordem elevada

1K e 2K no cálculo das forças internas iniciais, mas esses termos devem ser incluídos no

cálculo dos incrementos das forças internas. Assim, a matriz de rigidez tangente TK pode ser

obtida da seguinte forma:

[ ]τKKK += LT . (2.52)

2.8 VETOR DE FORÇAS INTERNAS

Uma das condições necessárias para a convergência do MEF baseia-se na correta

representação dos deslocamentos de corpo rígido pelas funções aproximadoras dos

deslocamentos.

Ao observar as funções aproximadoras adotadas para o elemento de pórtico,

verifica-se que as mesmas não descrevem adequadamente o movimento de rotação de corpo

27

rígido, pois a relação dxv∆=∆θ só é válida para o caso de pequenas rotações. Com isso a

condição de equilíbrio deixa de ser atendida.

Para reduzir os efeitos dessa incompatibilidade foram adotados os procedimentos

descritos a seguir:

a) Adoção de uma formulação com RLA para diminuir os problemas decorrentes

dos movimentos de corpo rígido. Assim, a deformada a ser calculada não se

distancia daquela tomada como referência e as parcelas de rotações de corpo

rígido são divididas em partes menores, sendo melhor aproximadas pelas

funções de interpolação;

b) Adoção da geometria deformada do elemento para o cálculo das forças

internas da estrutura.

Na formulação tratada neste capítulo, o incremento das forças internas foi obtido

através dos deslocamentos naturais incrementais, isto é, através dos deslocamentos que

realmente causam deformações, medidos em relação à última configuração de equilíbrio

obtida no processo de solução (YANG e KUO, 1994). Nesta abordagem, determina-se o vetor

de forças, Fitt ∆+ , de maneira incremental, obtendo-se a cada passo o acréscimo nas forças

internas. Para isso, considera-se válida a expressão:

FiFiFi tttt ∆∆+ += , (2.53)

onde Fit∆ é obtido diretamente por:

nt ∆uKFi =∆ , (2.54)

sendo K a matriz de rigidez dada pela Eq. 2.55:

⎥⎦⎤

⎢⎣⎡ +++= 21 6

121 KKKKK τL . (2.55)

O vetor de deslocamentos naturais incrementais do elemento, n∆u , no sistema

local, é definido como:

[ ]Tn 21 000 φδφ=∆u , (2.56)

onde, a partir da Fig. 2.8, chega-se à:

28

⎥⎦⎤

⎢⎣⎡ ∆−∆

=Ψ

∆−∆=∆−∆=

−

Lvv

vvvuuu

t121

12

12

tan

, (2.57)

sendo Ψ a rotação de corpo rígido que o elemento sofre e ∆u1, ∆v1, ∆u2 e ∆v2 as componente

de deslocamento nodais.

(a)

(b)

Figura 2.8 – Cálculo das forças internas: (a) deslocamentos incrementais; (b) geometria

deformada.

Assim, a partir dos valores intermediários u , v e Ψ , chega-se às expressões dos

deslocamentos naturais:

.,, 2211 Ψ−∆=Ψ−∆=−= ∆+ θφθφδ LL ttt (2.58)

Utilizando-se a Eq. 2.53 juntamente com a Eq. 2.54, chega-se ao vetor de forças

nodais absorvidas pela estrutura na configuração de equilíbrio tt ∆+ , determinado no sistema

local para o elemento genérico.

O vetor de forças internas deve ser transformado para o sistema global por meio

da matriz de rotação, R.

O vetor de forças internas do sistema é, então, obtido na forma global, somando-

se os esforços internos absorvidos por cada elemento, devidamente transformados para o

sistema global de coordenadas da forma:

FiRFi ttTGL

tt ∆+∆+ ∑= . (2.59)

29

2.9 ACOPLAMENTO DOS ELEMENTOS DE CABO E DE PÓRTICO

O problema de acoplamento entre os elementos de cabo (dois deslocamentos

livres por nó) e os elementos de pórtico plano (dois deslocamentos e uma rotação livres por

nó) é resolvido com a compatibilização dos graus de liberdade dos nós dos elementos de

cabos com os nós dos elementos de pórtico. Os elementos de cabo submetidos à tensão de

compressão não são considerados na resolução do sistema estrutural.

2.10 MATRIZ DE MASSA E DE AMORTECIMENTO

A determinação da matriz de massa e de amortecimento para um elemento está

baseada no método das massas consistentes apresentado por Paz (1992). Neste método, os

coeficientes das matrizes são calculados, aplicando o princípio dos trabalhos virtuais, a fim de

se determinar os efeitos axiais, de flexão, das forças internas, de velocidade e de aceleração

nos nós do elemento. Deste procedimento encontra-se que:

∫=L

0

T HHM ρ , (2.60)

∫=L

0

T HHC ξ , (2.61)

onde M é a matriz de massa consistente, C a matriz de amortecimento consistente, H a

matriz com as funções de forma, ρ a densidade do material e ξ a razão de amortecimento

viscoso distribuído ao longo do comprimento L do elemento.

Expandindo a Eq. 2.60 constrói-se a matriz de massa do elemento de pórtico plano

PM e a matriz de massa do elemento de cabo CM , escritas como:

⎥⎥⎥⎥⎥⎥⎥⎥

⎦

⎤

⎢⎢⎢⎢⎢⎢⎢⎢

⎣

⎡

−−−−

−−

=

²4220²313022156013540001400070

²3130²422013540221560007000140

420

LLLLLL

LLLLLL

LmPPM , (2.62)

30

⎥⎥⎥⎥

⎦

⎤

⎢⎢⎢⎢

⎣

⎡

=

2010020110200102

6LmC

CM , (2.63)

onde Pm é a massa por unidade de comprimento L, do elemento de viga/coluna e Cm é a

massa por unidade de comprimento L, do elemento de cabo/treliça.

2.11 AMORTECIMENTO

O amortecimento está presente em todos os sistemas oscilatórios, mas é difícil

prever seu real valor. Para isso, admitem-se modelos ideais de amortecimento que resultam

em respostas satisfatórias. Dentre esses modelos, a força de amortecimento viscoso conduz a

um tratamento matemático simples, razão pela qual é adotada no presente trabalho.

O amortecimento proporcional ou de Rayleigh é usualmente utilizado para

calcular a matriz de amortecimento C. Este esquema é uma combinação linear das matrizes de

rigidez K e massa M, isto é:

MKC βα += , (2.64)

onde α e β são, respectivamente, as constantes de amortecimento proporcional de rigidez e

de massa. A relação entre α , β e a freqüência ω é dada pela a razão de amortecimento ξ:

⎟⎠⎞

⎜⎝⎛ +=

ωβωαξ

21 . (2.65)

As constantes de amortecimento α e β são determinadas para duas diferentes

freqüências ( 1ω e 2ω ). Escolhendo a razão de amortecimento ξ e resolvendo simultaneamente

as equações:

( ) ( )( ) ( )2

1221221

21

2212

2

2

ωωωξωξωωβ

ωωωξωξα

−−=

−−= , (2.66)

sendo 1ω e 2ω , a freqüência natural da estrutura e a máxima freqüência de vibração obtidas,

respectivamente.

Considerando um sistema em vibração livre, o valor de ξ determina o caráter

oscilatório do sistema, passando-se a ter um movimento harmônico amortecido ou até sem

31

caráter oscilatório. Para o parâmetro ξ<1,0 tem-se um movimento oscilatório subamortecido,

quando ξ>1,0 o movimento é superamortecido e para ξ =1,0 tem-se o caso crítico. Neste

trabalho, adotou-se a razão de amortecimento, ξ, de 1%.

2.12 EQUAÇÃO DINÂMICA NÃO-LINEAR DE EQUILÍBRIO

Nas seções anteriores deste capítulo, as equações foram escritas em função das

coordenadas nodais. A equação de movimento, como função destas coordenadas, pode ser

escrita impondo condições de equilíbrio dinâmico entre as forças de inércia, as forças de

amortecimento, as forças elásticas e as forças externas, isto é:

)()()()( tttt FrFFF SDI λ=++ . (2.67)

As três forças no primeiro membro da Eq. 2.67 são expressas em função das

matrizes de massa M , de amortecimento C e de rigidez K do sistema. Assim, combinando

as Eqs. 2.60, 2.61 e 2.50, obtém-se a equação dinâmica não-linear de equilíbrio, na forma:

)()(61

21

21 tuL FrFi∆uKKKΚuCuM λτ =+⎥⎦⎤

⎢⎣⎡ +++++

•••

, (2.68)

onde, M é a matrizes de massa consistente, C é a matriz de amortecimento viscoso, LK

corresponde à parcela linear da matriz de rigidez, τK corresponde à parcela relativa às

tensões iniciais, 1K e 2Κ correspondem às parcelas não-lineares dependentes dos

deslocamentos, )(uFi é o vetor de forças internas dependentes dos deslocamentos, )(tFr é o

vetor de carga de referência dependentes do tempo, ••

u , •

u e ∆u são respectivamente, os

vetores de aceleração, de velocidade e de deslocamento das coordenadas nodais do sistema, e

λ é o parâmetro de carga.

Os métodos numéricos empregados neste trabalho para resolução da Eq. 2.68 são

descritos no capítulo três.

2.13 FREQÜÊNCIAS NATURAIS

Para avaliar a possibilidade de ocorrência de uma condição indesejável de quase

ressonância, faz-se necessário conhecer as freqüências dominantes de excitação do

carregamento externo a fim de se fazer uma comparação com as freqüências naturais da

estrutura. Estas freqüências naturais são obtidas a partir das vibrações livres e não

32

amortecidas, as quais são expressas pela equação diferencial de movimento, dada na Eq. 2.68,

na ausência do amortecimento e das forças externas, ou seja, de forma genérica:

0=+••

uKuM , (2.69)

cuja solução pode ser:

)cos( tωΦu = , (2.70)

onde ω é uma freqüência natural de vibração e Φ é o vetor de amplitudes dos deslocamentos

nodais.

Da substituição da Eq. 2.70 em 2.69 chega-se à seguinte equação característica do

problema de autovalor:

02 =− MK ω . (2.71)

A importância, dentro da análise estrutural, dos resultados de uma análise de

vibração livre, que são os modos de vibração e as freqüências naturais de vibração, advêm dos

seguintes fatos (NEVES, 1990):

a) As freqüências naturais e seus respectivos modos de vibração estão ligados

intrinsecamente às propriedades dinâmicas da estrutura, já que esses resultados

são obtidos das equações diferenciais que representam o equilíbrio entre as

forças elásticas e as forças inerciais;

b) As freqüências naturais possibilitam verificar se a estrutura corre risco de

entrar em ressonância, pelo confronto entre as freqüências de um determinado

carregamento variável no tempo e as freqüências naturais do sistema;

c) O conhecimento das freqüências naturais é importante para verificação dos

limites de utilização de algumas estruturas como, por exemplo, de edifícios

altos ou pontes estaiadas sob a ação de forças ambientais.

2.14 MODELO DISCRETO

A análise do modelo discreto estudado por Pasquetti (2003) é realizada a partir da

expressão da energia potencial total. Uma apresentação clara da aplicação do método da

energia para a análise de estabilidade estrutural é apresentada, por exemplo, em Croll e

Walker (1975) e Thompsom e Hunt (1984).

33

2.14.1 FORMULAÇÃO ESTÁTICA

A Fig. 2.9 apresenta uma barra rígida de comprimento L presa a duas molas

lineares e sujeita a dois tipos de carregamentos: carga axial P, disposta a uma altura H do

apoio da barra, e outra carga axial p, devida ao peso próprio, situada no centro de gravidade,

y , da barra. Cada mola tem inclinação α e está fixada à barra a uma distância h do apoio da

mesma. Mostra-se, também, a geometria da torre e o carregamento. Na Fig. 2.9a apresenta-se

a torre na configuração fundamental de equilíbrio cuja estabilidade se deseja analisar e na

Fig. 2.9b a estrutura sujeita a uma perturbação cinematicamente admissível, θ.

Figura 2.9 – Modelo de molas: (a) configuração inicial; (b) configuração perturbada.

A variação da energia potencial total, π∆ , é dada por:

WU ∆−∆=∆π , (2.72)

onde ∆U é a variação da energia interna de deformação e ∆W é a variação do trabalho das

cargas externas.

A energia interna de deformação da mola devido ao alongamento l∆ é:

l∆= 121 KU , (2.73)

onde 1K é a constante da mola.

O alongamento l∆ da mola é a soma de um alongamento inicial 0l∆ , com a mola

pré-tensionada em θ=0, e o alongamento causado durante a rotação θ da barra, θl∆ . Assim:

0lll ∆+∆=∆ θ , (2.74)

sendo:

34

1

00 K

F=∆l , (2.75)

onde 0F é a força que age na mola em θ=0, ou seja, de pré-tensionamento.

A variação do comprimento de uma mola i devido à rotação θ é:

( ) ( ) ( ) ( ) ⎟⎟⎟

⎠

⎞

⎜⎜⎜

⎝

⎛+−+⎟⎟

⎠

⎞⎜⎜⎝

⎛±=∆ 1

tan1cossen

tan1

22

2

iiii L

αθθ

αγθl , (2.76)

onde ( )θsen é positivo para mola que alonga e negativo para que encurta, θ é positivo no

sentido apresentado na Fig.2.9b e Lh=γ na Fig. 2.9a.

A variação do trabalho das forças externas W∆ , produzido pelas cargas externas

concentradas iP e pelo peso próprio p, durante uma rotação θ da barra é dada por:

( )( ) ( )( )∑=

−Γ+−=∆np

iii LPLypW

1cos1cos1 θθ , (2.77)

onde np é o número de cargas concentradas aplicadas na coluna e, ambos na Fig. 2.9a,

LHy = e Ly=Γ .

Substituindo as Eqs. 2.77 e 2.73 na Eq. 2.72, chega-se à seguinte expressão da

energia potencial total:

( )( )

( )( )∑

∑

=

=

−Γ+

+−−⎥⎦

⎤⎢⎣

⎡⎟⎟⎠

⎞⎜⎜⎝

⎛+∆+∆=∆

np

iii

nm

i i

iiiii

LP

LypKF

FK

1

1

00

2

cos1

cos1221

θ

θπ θθ ll

, (2.78)

onde nm é o número de molas.

Derivando a expressão da energia potencial total (Eq. 2.78) em relação à θ,

obtém-se a equação não-linear de equilíbrio do caminho pós-crítico. Isolando a carga

concentrada P1, tem-se:

( )( )

( ) ⎥⎥⎥

⎦

⎤

⎢⎢⎢

⎣

⎡

−

Γ−⎟⎟⎠

⎞⎜⎜⎝

⎛+∆+∆

Γ= ∑∑

==

θ

θθ

θθθ

sen

sen121

21

dd

sen1

11

00

2

11

yp

PLK

FFK

P

np

iii

nm

i i

iiiii ll

, (2.79)

35

A equação de equilíbrio pós-crítico de equilíbrio (Eq. 2.79) pode ser escrita em

uma forma adimensional. Para isso, considera-se o caso mostrado na Fig. 2.9, em que as

únicas cargas são P e p, e no qual há apenas duas molas. Substituindo-se o alongamento

correspondente a cada mola (Eq. 2.76) e considerando as molas idênticas, o caminho pós-

crítico, escrito em forma adimensional, é dado por:

( )( ) ( ) ( ) ⎟

⎟⎠

⎞⎜⎜⎝

⎛⎟⎟⎠

⎞⎜⎜⎝

⎛ −−

−+⎟

⎠

⎞⎜⎝

⎛−=

BD

AC

BAP11

101 sen1

sen111

sencos

ααξξ

θθλ , (2.80)

para:

( ) ( ) ( )( ) ( ) ( )( ) ( ) ( )( ) ( ) ( ) ,sensencos21

,sensencos21,tantansen21

,tantansen21

11

11

12

1

12

1

θααθαα

ααθ

ααθ

−=+=

+−=

++=

DCB

A

(2.81)

sendo:

.

,

,

,

211

1

1010

1

1

111

γξ

γξ

ξ

λ

=

=

=

Γ=

LKF

LKypLK

P

p

P

(2.82)

2.14.2 FORMULAÇÃO DINÂMICA

A equação de movimento pode ser obtida a partir da função de Lagrange ou

Lagrangiano L , que é dado pela diferença entre a variação da energia cinética T∆ , e a

variação da energia potencial total π∆ .

π∆−∆= TL . (2.83)

Para a Fig. 2.10, a massa de um elemento infinitesimal da coluna de comprimento

ds é dada por:

dsAdVdm ρρ == , (2.84)

onde ρ é a densidade do material e A é a área da seção transversal da coluna.

36

Figura 2.10 – Posição do elemento de barra com comprimento infinitesimal, ds.

A energia cinética é dada por:

∫=∆L

AdsvT0

2

2ρ , (2.85)

onde L é o comprimento da coluna e v é a velocidade.

É necessário determinar a posição do elemento infinitesimal de massa dm e

comprimento ds segundo as direções x e y indicadas na Fig. 2.10, a saber:

( )θsensx = , (2.86)

( )θcossy = , (2.87)

onde s é a distância do ponto à base da torre e θ, neste caso, é uma função do tempo.

Derivando as equações 2.86 e 2.87 em relação ao tempo, obtêm-se as

componentes do vetor velocidade nas direções x e y:

( ) ( )→•→•→

−= jsisv θθθθ sencos . (2.88)

Substituindo-se 2.88 em 2.84 e integrando-se ao longo de s, chega-se à expressão

que descreve a energia cinética para uma coluna de seção transversal constante:

23

61 •

=∆ θρ LAT . (2.89)

Substituindo-se a energia potencial total (Eq. 2.72) e a energia cinética (Eq. 2.89),

no Lagrangiano (Eq. 2.83), e aplicando-se o princípio de Hamilton, chega-se à equação de

37

movimento que descreve o equilíbrio dinâmico de um sistema conservativo, ou seja, sem

amortecimento, com vários deslocamentos nodais livres:

0=∂

∂−

∂

∂

∂

∂•

θθ

LL

t , (2.90)

sendo L o Lagrangiano.

Para um sistema conservativo, tem-se, pelo princípio da conservação da energia,

que a soma da energia cinética, T∆ , com a energia potencial total, π∆ , do sistema é uma

constante:

CT =∆+∆ π , (2.91)

onde C depende apenas das condições iniciais do sistema, ou seja, do deslocamento inicial

θ (0) e da velocidade inicial •

θ (0), da torre.

No caso de um sistema não conservativo, ou seja, com amortecimento, apenas

uma parcela da energia cinética é transformada em energia potencial e vice-versa, o restante é

dissipado pelo amortecimento. A equação de movimento é obtida igualando-se a Eq. 2.90 à

força de amortecimento, c:

•

•=

∂

∂−

∂

∂

∂

∂ θθθ

ct

LL . (2.92)

Para sistemas submetidos a oscilações de pequena amplitude, os termos da

equação de Lagrange podem ser reescritos, a fim de se chegar à equação geral de movimento

dinâmico dada por:

0=∆

++•••

θπθθ

ddcm , (2.93)

Finalmente, as equações dimensional e adimensional que descrevem o movimento

da barra quando a mesma esta conectada a molas são dadas, respectivamente, por:

( )

( ) ( ) 0sensen121

21

dd

sen32

3

1

0

10

2

1

133

=⎥⎦⎤−Γ−⎟⎟

⎠

⎞+

⎢⎣

⎡⎜⎝⎛ +∆+∆

Γ+⎟⎟

⎠

⎞⎜⎜⎝

⎛+

∑

∑

=

=

•••

θθ

θθθξωρθρ

θθ

ypPLK

F

FKPLALA

np

iii

i

i

nm

iiiii ll

, (2.94)

38

( )( )

( ) ( )0

sen1

sen1

11sencos

32

31

111

1

011

=−⎟⎟⎠

⎞⎟⎟⎠

⎞⎜⎜⎝

⎛ −−

−+

⎜⎜⎝

⎛+⎟⎟

⎠

⎞⎜⎜⎝

⎛−+⎟⎟

⎠

⎞⎜⎜⎝

⎛+

•••

PBD

AC

BAKLA

KLA

λαα

ξ

ξθθθξωρθρ

, (2.95)

onde ••

θ e •

θ são, respectivamente, a aceleração e a velocidade das coordenadas livres da

estrutura.

A Eq. 2.93 se resolve, neste trabalho, empregando os métodos numéricos

descritos no capítulo que se segue.

39

3 MÉTODOS NUMÉRICOS

3.1 INTRODUÇÃO

Na análise estrutural, a noção de estabilidade aparece diretamente associada ao

conceito de equilíbrio estático ou dinâmico e suas configurações. A estabilidade dessas

configurações pode ser avaliada através do estudo do comportamento da estrutura quando ela

sofre uma perturbação.

Uma forma de se representar as configurações de equilíbrio estático de um

determinado sistema estrutural consiste no traçado de um diagrama, fundamentalmente não-

linear, denominado trajetória de equilíbrio. Para representar as configurações de equilíbrio

dinâmico pode-se traçar, a partir das soluções de equilíbrio dos sistemas contínuos, a resposta

no tempo, as seções de Poincaré e as bacias de atração.

Neste capítulo serão mostrados alguns conceitos necessários à compreensão da

teoria da estabilidade estrutural e dos resultados obtidos neste trabalho, além dos métodos

numéricos usados na determinação da solução de sistemas de equações não-lineares, que

permitem obter as várias configurações de equilíbrio estático e dinâmico das estruturas

estudadas, apresentando diversos tipos de bifurcações e a estabilidade das soluções.

Descrevem-se também os métodos numéricos empregados na análise de sistemas estáticos e

dinâmicos, sendo que, os algoritmos utilizados foram obtidos a partir de Galvão (2000) e Del

Prado (2001).

40

3.2 PRINCÍPIOS DA TEORIA DE ESTABILIDADE ESTRUTURAL

Na análise não-linear geométrica de estabilidade, as equações de equilíbrio são

estabelecidas na configuração deformada, no entanto, realiza-se a linearização dessas

equações com relação aos deslocamentos envolvidos e o equilíbrio é considerado na

configuração indeformada. Esse tipo de análise permite encontrar a chamada solução

fundamental e fornece o valor das cargas críticas e de seus respectivos modos de flambagem.

Outro aspecto importante da análise da estabilidade estrutural envolve a definição

do comportamento do sistema após a passagem para o caminho pós-crítico. Nesta análise,

denominada análise não-linear de estabilidade, são considerados os termos não-lineares nas

equações de equilíbrio, que permitem estudar o comportamento da estrutura quando se têm

grandes deslocamentos ou grandes deformações.

O estudo da estabilidade baseado no critério da energia relaciona a energia interna

de deformação e o trabalho realizado pelas forças externas da estrutura (Item 2.5). Para

atender a condição de equilíbrio segundo o princípio da energia potencial total estacionária, o

funcional da energia potencial total Π , deve encontrar-se numa posição de mínimo, assim:

,0,,0

,

=∂

Π∂∀=

∂Π∂

=Π

+=Π

∑i

iii u

uuu

WU

δδδ (3.1)

onde Π representa a energia potencial do sistema, U a energia interna de deformação, W o

trabalho das forças externas e iu representa os deslocamentos generalizados (translações e

rotações).

O tipo de equilíbrio de uma dada configuração é definido pela segunda variação

da energia potencial. Caso se tenha 02 >Πδ , o equilíbrio é dito estável, se 02 <Πδ o

equilíbrio é instável e para 02 =Πδ o equilíbrio é indiferente, isto é, o sistema pode ser tanto

estável quanto instável.

A grande dificuldade de se obter a solução exata das equações, que regem os

fenômenos de instabilidade na maioria dos sistemas estruturais contínuos, fez com que fossem

desenvolvidos métodos numéricos aproximados de análise de estruturas. Esses métodos

baseiam-se na substituição da análise da estrutural real contínua por uma análise estrutural em

pontos discretos, fornecendo uma solução aproximada para o sistema contínuo original.

Dentre os métodos mais utilizados têm-se:

41

a) Método das diferenças finitas;

b) Método de Engesser-Newmark;

c) Método de Galerkin;

d) Método de Rayleigh-Ritz;

e) Método dos elementos finitos (MEF).

Enquanto os três primeiros são aplicados às equações diferenciais de equilíbrio, os

dois últimos envolvem a consideração da energia potencial. O MEF pode igualmente ser

formulado a partir das equações de equilíbrio. Uma descrição mais aprofundada sobre a teoria

da estabilidade das estruturas pode ser encontrada em Bazant e Cedolin (1991).

No método MEF, ao se discretizar uma estrutura contínua, a determinação do seu

comportamento não-linear é transferida para a resolução de um sistema de equações de

equilíbrio com a forma da Eq. 2.50.

A solução do sistema de equações não-linear da Eq. 2.50 é obtida por meio de

pequenos incrementos de deslocamentos que conduzem à condição de equilíbrio do sistema.

Assim, para cada incremento de deslocamento o sistema é resolvido em sua forma linearizada.

F∆uK = , (3.2)

onde K representa a matriz de rigidez tangente, ∆u é o vetor dos deslocamentos nodais

generalizados e F é o vetor das forças nodais generalizadas calculadas a partir das cargas

aplicadas à estrutura.

3.3 DETERMINAÇÃO DA TRAJETÓRIA DE EQUILÍBRIO

Por ser uma representação gráfica da resposta de uma estrutura submetida à ação

de forças externas, a obtenção da trajetória de equilíbrio utilizando-se o MEF depende da

solução do sistema de equações definidos na Eq. 3.2. Dependendo do tipo de análise estrutural

a ser realizada (linear ou não-linear geométrica), a matriz de rigidez K pode apresentar uma

dependência linear ou não-linear com relação aos deslocamentos nodais ∆u. Na análise não-

linear geométrica da estrutura, a dependência da matriz de rigidez em relação aos

deslocamentos é altamente não-linear e a solução do sistema representado na Eq. 2.50 é

obtida usualmente através de técnicas com procedimentos incrementais e iterativos.

Estas técnicas de solução devem ser apropriadas e capazes de superar os

problemas numéricos associados ao comportamento não-linear dos sistemas estruturais. No

42

contexto da implementação computacional, devem ter a capacidade de detectar pontos

críticos, tais como pontos limites e pontos de bifurcação, e seguir a trajetória de equilíbrio

além desses pontos.

Fenômenos associados ao comportamento não-linear de uma estrutura podem

ocorrer quando a estrutura encontra-se na configuração de equilíbrio situada em um ponto

limite e é submetida a um ligeiro aumento de carga. Nesta situação, alguns sistemas

estruturais desenvolvem uma redução progressiva da rigidez da estrutura até que esta redução

se anule em outro ponto limite. Este comportamento denominado de salto ocorre de duas

formas:

a) Salto dinâmico sob controle de carga (snap-through), conforme mostrado na

Fig. 3.1, caracteriza-se pela redução da rigidez da estrutura acompanhada da

redução da carga, tendo, no entanto, um acréscimo dos deslocamentos;

b) Salto dinâmico sob controle do deslocamento (snap-back), conforme mostrado

na Fig. 3.2, caracteriza-se pela redução da rigidez da estrutura acompanhada

da redução da carga. No entanto, ocorre a inversão do crescimento dos

deslocamentos e a trajetória de equilíbrio retrocede com uma curvatura

positiva.

Deslocamento

Carga

A B

Figura 3.1 – Trajetória de equilíbrio com salto dinâmico do tipo snap-through.

Deslocamento

A

B

Carga

Figura 3.2 – Trajetória de equilíbrio com salto dinâmico do tipo snap-back.

43

Para evitar problemas que possam ocorrer devido ao comportamento não-linear

das estruturas e determinar toda a trajetória de equilíbrio, os métodos numéricos devem:

a) Ser auto-adaptativos na mudança do sentido de crescimento de carga nos

pontos limites;

b) Manter a estabilidade numérica das iterações tanto na região estável quanto na

região instável da trajetória de equilíbrio;

c) Ajustar automaticamente os tamanhos dos passos para retratar a variação na

rigidez da estrutura.

Diversos procedimentos de solução têm sido propostos para traçar as trajetórias de

equilíbrio até e além dos pontos críticos. Muitos são baseados em variantes do método de

Newton-Raphson e incorporam diferentes técnicas, tais como: controle de deslocamento,

controle de energia, comprimento do arco, controle de deslocamento generalizado, entre

outras.

No presente trabalho, para a determinação dos pontos críticos e obtenção da

trajetória de equilíbrio da estrutura, foi utilizado o método de Newton-Raphson, tendo como

técnicas de controle o método do comprimento de arco.

3.4 MÉTODO DE SOLUÇÃO PARA ANÁLISE NÃO-LINEAR

Na análise não-linear das estruturas, a solução da Eq. 2.50 é obtida de forma

incremental, ou seja, para uma seqüência de incrementos do parâmetro de carga 1λ∆ , 2λ∆ ,

3λ∆ , é calculada uma seqüência de incremento dos deslocamentos nodais 1u∆ , 2u∆ , 3u∆ ...

Como Fi é uma função não-linear dos deslocamentos, a solução do problema

( λ∆ , ∆u ) não satisfaz inicialmente a Eq. 2.50. Tem-se então uma força residual g definida

como:

)(uFiFrg −= λ .

(3.3)

Assim, define-se como solução da Eq. 2.50, o vetor de deslocamentos nodais u

que minimiza o resíduo g estabelecido na Eq. 3.3. Para isso, os algoritmos existentes

apresentam como passo fundamental a avaliação das forças residuais através da relação:

guK =δT , (3.4)

44

onde TK é a matriz de rigidez tangente do sistema estrutural e uδ é o vetor de deslocamentos

nodais residuais.

O cálculo do vetor de deslocamentos nodais u é realizado por meio de um

processo preditor-corretor, conforme apresentado a seguir:

a) A primeira etapa é definida a partir da configuração de equilíbrio da estrutura.

Seleciona-se um incremento inicial do parâmetro de carga, 0λ∆ , com base na

estratégia de incremento de carga utilizada e, a seguir, determina-se o

incremento inicial dos deslocamentos nodais 0u∆ . As aproximações 0λ∆ e 0u∆ caracterizam a chamada solução incremental predita.

b) Na segunda etapa procura-se, através de estratégias de iterações do tipo

Newton-Raphson, corrigir a solução incremental inicialmente proposta na

etapa anterior, com o objetivo de restaurar o equilíbrio da estrutura o mais

rápido possível.

3.4.1 MÉTODO DE NEWTON-RAPHSON

O método de Newton-Raphson provavelmente é o algoritmo numérico iterativo

mais utilizado na solução de equações não-lineares. Muitos métodos utilizados para a solução

de equações não-lineares podem ser considerados como variações deste método. A técnica do

método consiste em manter o parâmetro de carga λ sempre positivo e constante ao longo de

um ciclo iterativo de um incremento de carga. Para que seja satisfeito o equilíbrio dos

elementos finitos deve-se encontrar a solução da equação:

0=− FiR , (3.5)

onde R é o vetor de cargas externas aplicadas nos nós e Fi é o vetor de forças internas

também aplicadas nos nós, correspondentes às tensões atuantes no elemento.

Como o vetor de forças internas Fi, é dependente dos deslocamentos nodais de

uma forma não-linear, torna-se necessário aplicar um método iterativo para se obter a solução

da Eq. 3.5. Assim, na iteração k, para um determinado deslocamento nodal u é obtido um

vetor resíduo, ∆R , na Eq. 3.5 dado por:

1−−= kk FiR∆R , (3.6)

onde o vetor resíduo k∆R é obtido por:

45

kkk ∆R∆uK =−1 , (3.7)

onde 1−kK é a matriz de rigidez tangente definida em função dos deslocamentos nodais da

iteração anterior k-1 e k∆u o incremento dos deslocamentos nodais na iteração k.

Assim, os deslocamentos nodais totais na iteração k são calculados pela equação:

1−+= kkk u∆uu . (3.8)

O processo iterativo apresentado nas Eqs. 3.6 e 3.7, continua até que seja satisfeito

um critério de convergência. É importante destacar uma característica importante do método

de Newton-Raphson. A cada iteração é calculada uma nova matriz de rigidez, implicando em

formação e triangularização da matriz de rigidez, exigindo considerável tempo computacional

no processo de solução. O método de Newton-Raphson modificado reduz este problema,

empregando-se uma mesma matriz de rigidez tangente até a convergência.

De acordo com Yang e Kuo (1994), bem como, Crisfield (1991), falhas do

método de Newton-Raphson são detectadas na resolução de problemas envolvendo pontos

limites. Isto pode ser atribuído ao fato de que as iterações são realizadas com carga constante.

Como pode ser visto na Fig. 3.3a e Fig. 3.3b, a direção das iterações é controlada pela linha

horizontal que representa o nível atual das cargas aplicadas. Sempre que as cargas aplicadas

excedem as cargas correspondentes ao ponto limite em um determinado incremento, a linha

horizontal que guia a direção das iterações não retorna a trajetória de equilíbrio. Como

resultado, a convergência não pode ser alcançada pelo método nos pontos próximos aos

pontos limites. Devido a essa razão o método de Newton-Raphson não é recomendado para

determinar a resposta da estrutura além dos pontos limites, sendo usados nestes casos métodos

associados ao Newton-Raphson.

3.4.2 SOLUÇÃO INCREMENTAL PREDITA

A obtenção da solução incremental inicial ou predita tem como passo fundamental

a definição do parâmetro de carga inicial ∆λ0. A seleção automática do tamanho do

incremento desse parâmetro é importante e deve refletir o grau de não-linearidade corrente do

sistema estrutural em estudo. Em outras palavras, uma estratégia eficiente de incremento

automático de carga deve satisfazer basicamente aos seguintes requisitos:

a) Produzir grandes incrementos quando a resposta da estrutura for

aproximadamente linear;

46

b) Gerar pequenos incrementos quando a resposta da estrutura for fortemente

não-linear;

c) Ser capaz de escolher o sinal correto para o incremento, introduzindo medidas

capazes de detectar quando os pontos limites são ultrapassados.

Rt

∆R

u1∆

1

∆R

2

2∆u

t u ut+∆t

inclinação: K1t+∆t

t+∆tinclinação: 2K

Deslocamento

Carga

t+∆tR

(a)

tR1

ut

∆u

Deslocamento

2u∆

ut+∆t

Rt+∆t

∆R1

Carga

t+∆t

2∆R

inclinação: K1

(b)

Figura 3.3 – Método de Newton-Raphson: (a) Padrão; (b) Modificado.

Para a obtenção da solução incremental predita ∆λ0 e ∆u0, primeiramente, monta-

se a matriz de rigidez tangente da estrutura K, usando informações da última configuração de

equilíbrio. Após a definição de K, resolve-se o sistema de equações lineares:

FruK =rδ . (3.9)

Obtêm-se então os deslocamentos nodais tangentes δur. Pelo fato de Fr ser um

vetor de magnitude arbitrária, os deslocamentos nodais incrementais tangenciais ∆u0, serão

obtidos escalonando-se δur, como indicado:

rδu∆u 00 λ∆= , (3.10)

onde ∆λ0 selecionado automaticamente pela estratégia de incremento de carga adotada.

Nessa etapa, o parâmetro de carga e os deslocamentos nodais totais são

atualizados como mostrado nas equações:

,,

01

01

uuu ∆+=

∆+=−

−

kk

kk λλλ (3.11)

47

onde λk-1 e uk-1 caracterizam o ponto de equilíbrio obtido no último passo de carga, como

indicado na Fig. 3.4. Devido a não-linearidade do problema, as soluções dadas nas Eq. 3.11

raramente satisfazem à condição de equilíbrio total do sistema. Assim, iterações de correção

tornam-se necessárias para restaurar o equilíbrio.

u

λ

k

k

∆u

∆

∆

2u

1u

0

u1δ

restrição

uδ 2

Deslocamento

∆λ 1∆λ

0δλ

∆λ2

2

1δλ

Carga

nova solução

solução preditaTrajetória de equilíbrio

solução anterior

Figura 3.4 – Solução incremental-iterativa em sistema com um grau de liberdade.

3.4.3 CICLO DE ITERAÇÕES

No esquema tradicional do método de Newton-Raphson, o parâmetro de carga λ é

mantido constante durante o ciclo iterativo. Porém, caso se pretenda acompanhar todo o

traçado da trajetória de equilíbrio, com possíveis passagens pelos pontos limites, é necessário

permitir a variação de λ a cada iteração.

Seguindo essa proposta, onde é permitida a variação do parâmetro de carga, a

alteração nos deslocamentos nodais é representada pela equação de equilíbrio:

guK −=− kk δ1 . (3.12)

Observando-se a Eq. 3.12 nota-se que g é função dos deslocamentos nodais totais

uk, calculados na última iteração, e do valor corrente do parâmetro de carga total λk dado por:

kkk δλλλ += −1 , (3.13)

onde δλk é a correção do parâmetro de carga.

48

Substituindo-se a Eq. 2.50 e a Eq. 3.13 na Eq. 3.12, chega-se a:

( )( )FiFruK −+−= −− kkkk δ δλλ 11 , (3.14)

onde o produto Fr1−kλ caracteriza o vetor das forças externas total atuante na última

iteração.

Reescrevendo-se a Eq. 3.14 obtém-se:

FrguK kkkk δ δλ+−= −− 11 , (3.15)

que é a equação necessária para se trabalhar durante o ciclo iterativo. Da Eq. 3.15, tem-se que

os deslocamentos nodais iterativos podem ser decompostos em duas parcelas:

kr

kkg

k δδδ uuu δλ+= , (3.16)

onde:

.

,1

11

FrKu

gKu−

−−

=

−=k

r

kkg

δ

δ (3.17)

Deve-se notar que δugk é a correção que teria sido obtida da aplicação do método

de Newton-Raphson com a estratégia convencional de incremento do parâmetro de carga

constante e que δurk é o vetor de deslocamentos iterativos resultantes da aplicação de Fr.

Caso seja adotado o método de Newton-Raphson modificado, o vetor de

deslocamentos δurk será igual ao vetor de deslocamentos tangentes δur, calculado no

Item 3.4.2, não se modificando durante as iterações, pois K não se altera.

O valor da correção do parâmetro de carga δλk, na Eq. 3.12, é determinado

seguindo a estratégia de iteração adotada, onde se introduz uma equação de restrição que deve

ser respeitada a cada iteração. Com a determinação de δλk retorna-se à Eq. 3.16 para o cálculo

da correção dos deslocamentos.

Com a obtenção da solução iterativa (δλk, δuk), faz-se a atualização das variáveis

incrementais do problema.

.,

1

1

kkk

kkk

δuuu +∆=∆

+∆=∆−

− δλλλ (3.18)

Para o parâmetro de carga e os deslocamentos nodais totais tem-se que:

49

.,

1

1

kkk

kkk

uuu ∆+=

∆+=−

− λλλ (3.19)

Os procedimentos descritos da Eq. 3.15 à Eq. 3.19 são repetidos até que o critério

de convergência seja satisfeito.

3.4.4 MÉTODO DE NEWTON-RAPHSON ASSOCIADO À TÉCNICA DO COMPRIMENTO

DE ARCO

A idéia básica do método consiste em tratar o parâmetro de carga como uma

variável adicional, controlando não o incremento do parâmetro de carga, nem o incremento de

uma determinada componente do deslocamento, mas sim o comprimento do vetor que une o

ponto conhecido da trajetória ao ponto incógnito desejado.

Assim, para equilibrar o número de equações e o número de incógnitas, uma

equação de restrição é somada as equações de equilíbrio originais, tendo a forma:

2222 l∆=∆+ Fr∆u λβα , (3.20)

onde ∆u é a norma do vetor de deslocamento incremental, λ∆ é incremento do parâmetro

de carga, Fr é a norma do vetor de carga de referência, l∆ é comprimento do arco e,

finalmente, α e β são fatores de escala ajustáveis que podem ser empregados para

homogeneizar as dimensões e a magnitude numérica das parcelas da equação. Uma descrição

gráfica do método do comprimento de arco aplicado a um sistema com um grau de liberdade

pode ser vista na Fig. 3.5.

Há diferentes versões da técnica do controle do arco, correspondendo a diferentes

valores atribuídos aos fatores de escala α e β . No caso da técnica do comprimento de arco

cilíndrico, tem-se 1=α e 0=β . Desta forma, a Eq. 3.20 pode ser reescrita como:

22 l∆=∆u , (3.21)

sendo:

kkk u∆u∆u δ+= −1 . (3.22)

50

Trajetória de equilíbrio

Deslocamento

solução anterior

solução predita

Carga

nova solução

restrição

∆l

Figura 3.5 – Técnica do comprimento de arco aplicada a um problema com um grau de

liberdade.

Realizando-se a substituição da Eq. 3.21 na Eq. 3.22, chega-se a uma equação

quadrática em função de δλ , ou seja:

02 =++ CBA δλδλ , (3.23)

onde os coeficientes A, B e C têm a seguinte forma:

( )( )( )( ) ( )( ) .

,2

,

211

1

lC

B

A

kg

kTkg

k

kg

kTr

rTr

∆−++=

+=

=

−−

−

u∆uu∆u

u∆uu

uu

δδ

δδ

δδ

(3.24)

A estratégia de incremento de carga baseado no incremento do comprimento de

arco como proposto por Crisfield (1991) pode ser definida como:

2/1

⎟⎠⎞

⎜⎝⎛∆=∆

II

ll tdt , (3.25)

onde t∆l e ∆l representam, respectivamente, os incrementos do comprimento de arco no passo

de carga anterior (valor conhecido) e no passo de carga corrente (incógnito), Id é o número de

iterações desejadas para o processo iterativo corrente, especificado pelo usuário, e tI é o

número de iterações necessárias para convergência no passo de carga anterior.

Através da Eq. 3.25 e da condição de restrição escrita para a solução incremental

inicial, tem-se:

51

( ) 200 lT∆=∆u∆u . (3.26)

Substituindo-se a Eq. 3.10 na Eq. 3.26, chega-se à expressão do incremento inicial

do parâmetro de carga:

rTr δδl

uu∆

±=∆ 0λ . (3.27)

O critério utilizado para escolher o sinal correto na Eq. 3.27 é o mesmo sugerido

por Crisfield (1991). No programa computacional utilizado neste trabalho, o usuário deve

especificar ∆λ0 como dado de entrada, sendo este valor usado em seguida para calcular ∆u0

através da Eq. 3.10. Substituindo-se, então, ∆u0 na Eq. 3.26, chega-se a ∆l1. Para os

incrementos de carga seguintes, calcula-se automaticamente ∆l através da Eq. 3.25.

A determinação da correção do parâmetro de carga, δλ, tem a função de otimizar

a convergência do processo iterativo, tornando a estratégia eficiente computacionalmente.

Isso significa que, para um dado passo de carga, a configuração de equilíbrio do sistema

estrutural em estudo deve ser obtida da forma mais rápida possível.

Para a técnica do comprimento de arco, a solução corretiva é determinada

calculando-se a solução da equação quadrática em δλ dada na Eq. 3.23. Com a resolução,

chega-se aos dois valores de δλ1 e δλ2, de forma que se deve escolher entre as soluções:

( )

( ) .

,

21

2

11

1

rkg

k

rkg

k

uuuu

uuuu

δδλδ

δδλδ

++∆=∆

++∆=∆−

−

(3.28)

Entre as duas soluções apresentadas na Eq. 3.28, escolhe-se aquela que

corresponda à menor mudança na direção do vetor de deslocamento da iteração k, ∆uk, em

relação ao vetor de deslocamento anterior ∆uk-1. Essa escolha deve prevenir um possível

retorno, o que faria a solução regredir ao longo do caminho já calculado. Um procedimento

utilizado para esse fim consiste em se obter o menor ângulo entre ∆uk e ∆uk-1. Isto equivale a

pesquisar o valor máximo do co-seno do ângulo θ.

( )( )

( )( ) ( )( )( ) ( )( ).cos

,cos

2

1

2,12

11

2,1

2

1

2,1

ll

l

rTkk

gTkTk

kTk

∆∆

+∆

+∆∆=

∆∆∆

=

−−−

−

uuuuu

uu

δδλ

δθ

θ (3.29)

52

Como a Eq. 3.23 é quadrática, ela poderá ter raízes imaginárias se ACB 42 − for

negativo. Essa situação existirá quando o incremento inicial do parâmetro de carga for muito

grande, ou se a estrutura possuir múltiplos caminhos de equilíbrio em torno de um ponto.

3.4.5 CRITÉRIOS DE CONVERGÊNCIA

O processo iterativo descrito no Item 3.4.4 deste capítulo termina indicando uma

nova posição de equilíbrio para a estrutura em análise quando pelo menos um dos três

critérios de convergência for atendido.

a) O primeiro critério é baseado em relações de forças e é calculado no início da

iteração corrente utilizando parâmetros da iteração anterior.

( )

( ) ζλ

ζ ≤∆

=−

−

Fr

g1

1

1 k

k

, (3.30)

onde ( )1−kg é a norma euclidiana do vetor das forças desequilibradas, que é calculada

usando-se o parâmetro de carga e os deslocamentos nodais totais da iteração anterior, ( )Fr1−∆ kλ é a norma euclidiana do vetor de incremento de carregamento externo e ζ é o fator

de tolerância.

b) O segundo critério de convergência obedece à relação de deslocamentos e é

sempre verificado no final da iteração corrente.

ζδ

ζ ≤∆

=ku

u2 , (3.31)

onde uδ é a norma euclidiana dos deslocamentos iterativos residuais, ku∆ é a norma

euclidiana dos deslocamentos incrementais, que são obtidos após a correção do processo

iterativo e ζ segue a mesma definição do primeiro critério.

c) O terceiro critério de convergência consiste em obedecer a ambas as relações

de força e de deslocamento dadas na Eq. 3.33 e na Eq. 3.34. Assim, este

critério é verificado se:

,, 21 ζζζζ ≤≤ (3.32)

onde ζ, ζ1, ζ2 são definidos nas alíneas (a) e (b).

53

3.5 MÉTODOS DE SOLUÇÃO PARA ANÁLISE DINÂMICA

3.5.1 CONCEITOS BÁSICOS

Um sistema dinâmico pode ser definido como um conjunto de equações

agrupadas por alguma interdependência, de modo que existam relações de causa e efeito entre

fenômenos que ocorrem e os elementos desse conjunto, quando algumas grandezas que

caracterizam seus objetos constituintes variam no tempo.

Assim, um sistema pode ser de tempo contínuo ou de tempo discreto. Um sistema

é de tempo discreto se o tempo t é um número inteiro, ou seja, se Zt ∈ . A evolução de um

sistema de tempo discreto é governada por uma ou mais equações de diferenças finitas na qual

o tempo t só pode assumir determinados valores. Já em um sistema onde a variável temporal é

contínua, o tempo t é um numero real, ou seja, Rt ∈ . A evolução de um sistema contínuo é

comandada por uma ou mais equações diferenciais.

Um sistema pode ser considerado como autônomo se as variáveis dependentes do

sistema não forem explicitamente dependentes da variável t, por exemplo, o sistema de

equações diferenciais da forma ),( µxFx =&& não apresenta o tempo de forma explícita e, nesse

caso, o sistema de equações é dito autônomo. Caso o sistema de equações diferenciais seja

explicitamente dependente do tempo na forma ),,( tx µFx =&& , o sistema é chamado de não-

autônomo.

3.5.2 EVOLUÇÃO DA RESPOSTA NO TEMPO E PLANO FASE

Observar a evolução no tempo da resposta de um sistema dinâmico submetido a

uma determinada condição inicial é fundamental para o completo conhecimento de seu

comportamento.

As respostas no tempo mostram como variam os deslocamentos e as velocidades

de um determinado grau de liberdade. Neste contexto, o algoritmo do tipo Newmark é o

método mais utilizado para integração implícita, na variável tempo, da equação de equilíbrio

dinâmico, escrita na Eq. 2.68 e reescrita de forma genérica no instante tt ∆+ como:

FruKuCuM λtttttttt ∆+∆+•

∆+••

∆+ =++ . (3.33)

Segundo o método de Newmark, as expressões para o deslocamento e a

54

velocidade no tempo tt ∆+ são expressas por:

( ) tt tttttt ∆+∆−+=••

∆+••••

∆+ uuuu γγ1 , (3.34)

22

21 ttt ttttttt ∆+∆⎟

⎠⎞

⎜⎝⎛ −+∆+=

••∆+

•••∆+ uuuuu ββ , (3.35)

adotando neste modelo 21=γ e 41=β o que representa uma aceleração constante no

intervalo de tempo t∆ igual a média entre a aceleração no tempo final e inicial. O intervalo de

tempo utilizado neste algoritmo é pequeno ( 5Tt ≤∆ ), a fim de garantir a estabilidade

numérica do método, sendo T o período de vibração natural da estrutura.

Das Eqs. 3.33, 3.34 e 3.35, chega-se ao seguinte sistema:

( )

⎟⎠⎞

⎜⎝⎛ +++

+⎟⎠⎞

⎜⎝⎛ +++=++

•••

•••∆+∆+

uuuC

uuuMFruCMK

ttt

ttttttt

aaa

aaaaa

541

32010 λ , (3.36)

ou de forma mais resumida:

^^F∆uK = , (3.37)

sendo ^

K a matriz de rigidez efetiva, ^F o vetor de carregamento efetivo e a1, a2, a3, a4 e a5

as constantes de interação dadas por:

.212

5,14,2

13

,12,1,10

2

2

⎟⎟⎠

⎞⎜⎜⎝

⎛−

∆=−=

∆=

∆=

∆=

∆=

ββγ

β

ββγ

β

taat

a

ta

ta

ta

(3.38)

Um resumo dos procedimentos do método de Newmark para problemas

estruturais é apresentado a seguir:

a) De posse da matriz de massa M, de rigidez tangente K e de amortecimento

proporcional C, inicializa-se os vetores ••

u0 , •

u0 e u0 ;

b) Para cada passo de tempo t∆ calculam-se as constantes de interação

5e4,3,2,1,0 aaaaaa ;

55

c) Monta-se a matriz de rigidez efetiva ^

K , calcula-se o vetor de carregamento

efetivo ^F e avalia-se os deslocamentos nodais incrementais

^^F∆uK = ;

d) No processo interativo para obtenção do equilíbrio dinâmico (k=1..) calcula-se

as aproximações para acelerações, velocidades e deslocamentos:

kktt

ttkktt

ttkktt

a

aaa

aaa

∆uuu

uu∆uu

uu∆uu

+=

−−=

−−=

∆+

••••∆+

•••••∆+

0

541

320

, (3.39)

calcula-se o vetor de forças internas:

11 ++∆+ += ktktt ∆uKFF ii , (3.40)

calcula-se o vetor de forças residuais:

⎟⎠⎞

⎜⎝⎛ ++−= ∆+

•∆+

••∆+∆++∆+ ktttkttkttttktt

iFuCuMFrR λ1 , (3.41)

calcula-se a correlação dos deslocamentos incrementais:

11^

+∆++ = kttk RδuK , (3.42)

calcula-se os deslocamentos incrementais:

11 ++ += kkk δu∆u∆u , (3.43)

verifica-se a convergência:

≤+ +

+

1

1

kt

k

∆uu

∆u fator de tolerância, (3.44)

e reavalia-se os deslocamentos, velocidades e acelerações;

e) Para o próximo instante avalia-se o vetor de forças internas iii FFF tttt ∆∆+ +=

e seleciona-se um novo intervalo de tempo t∆ .

O plano fase descreve a variação, para um determinado grau de liberdade, da

velocidade em função de seu próprio deslocamento. Seja, pois um sistema dinâmico discreto

com k graus de liberdade cujo movimento é descrito por k equações diferenciais de segunda

56

ordem, da forma:

,.,..,2,1,,...,,,,,...,,, 321321 kiuuuuuuuu kki =⎟⎠⎞

⎜⎝⎛=

••••••

Uu (3.45)

onde iU são funções não-lineares dos deslocamentos generalizados ( iu ) e das velocidades

generalizadas ( i

•

u ) que definem o comportamento do sistema em um dado instante.

Considerando as velocidades i

•

u como variáveis auxiliares, pode-se transformar a

Eq. 3.45 em um sistema de n = 2k equações diferenciais de movimento de primeira ordem, a

saber:

( )

( )

( ) ,,...,,

,,...,,

,,...,,

21

2122

2111

nnn

n

n

xxxFdt

dx

xxxFdt

dx

xxxFdtdx

=

=

=

M

(3.46)

onde x1, x2, ..., xn são as componentes do vetor x de dimensão n (x ∈ Rn) tendo como

elementos os deslocamentos e as velocidades generalizadas, que é chamado de vetor de

estado, pois descreve o estado do sistema em um instante t; e F1, F2, ..., Fn são as

componentes do vetor F de funções não-lineares do sistema.

O sistema de equações diferenciais de primeira ordem varia de acordo com o

tempo, portanto, um ponto P (x1, x2, ..., xn) qualquer move-se no espaço solução no qual o

tempo aparece explicitamente gerando uma trajetória denominada resposta no tempo. É

comum avaliar a evolução ao longo do tempo em um espaço de menor dimensão onde o

tempo aparece implicitamente. Este espaço é denominado plano fase e para sua determinação

toma-se como coordenadas as componentes do vetor de estado x.

O número de graus de liberdade envolvido na análise influencia na dimensão do

espaço solução. Por exemplo, sendo k o número de graus de liberdade do sistema estrutural, o

espaço solução tem dimensão 2k + 1, enquanto que, o plano fase tem dimensão 2k. Estes

espaços e suas respectivas projeções permitem ao projetista visualizar a evolução do sistema

dinâmico ao longo do tempo e compreender seu comportamento. A Fig. 3.6 ilustra a resposta

no tempo e o correspondente plano fase para um sistema com um grau de liberdade.

57

Tempo Deslocamento

Vel

ocid

ade

Tempo Deslocamento

Vel

ocid

ade

(a)

Deslocamento

Vel

ocid

ade

Deslocamento

Vel

ocid

ade

Deslocamento

Vel

ocid

ade

(b)

Figura 3.6 – Resposta no tempo (a) e plano fase (b).

A mudança de um estado do sistema em to = 0 a um outro estado em um tempo t é

determinada pelas equações diferenciais de movimento e denominada fluxo fase. As

configurações do sistema em cada intervalo de tempo dependem das condições iniciais

impostas. Essa influência das condições iniciais na resposta é uma das características dos

sistemas dinâmicos regidos por equações diferenciais não-lineares. Quando o tempo ∞→t e

o comportamento do sistema é assintótico, diz-se que o estado do sistema é permanente ou

steady state, enquanto que o comportamento do sistema antes de alcançar o estado

permanente é chamado de estado transiente.

Assim, um ponto escolhido arbitrariamente no instante inicial, descreve uma

trajetória tangente, em cada ponto, ao campo vetorial da velocidade, conforme a Fig. 3.7. A

representação do movimento no espaço fase é chamada trajetória ou órbita e o conjunto de

todos os movimentos possíveis é denominado fluxo fase. Quando se estudam sistemas com

mais de duas dimensões, a representação da resposta do sistema é feita através de projeções.

x1

x3

x2

t = t0

F(x)

x1

x3

x2

t = t0

F(x)

Figura 3.7 – Trajetória e velocidade no espaço fase 3D.

58

3.5.3 SOLUÇÕES DE EQUILÍBRIO

Seja o sistema autônomo:

),( µxFx =& · (3.47)

Neste sistema o vetor F não depende explicitamente da variável tempo t. Assim, a

Eq. 3.47 é chamada de invariante com o tempo ou estacionária. Isto significa que x(t) é uma

solução do sistema e então x(t+τ ) é também uma solução do sistema para qualquer valor

arbitrário de τ .

Soluções do sistema de equações diferenciais são denominadas pontos fixos,

sendo obtidas a partir do sistema de equações:

0),( =µxF . (3.48)

Nos problemas de análise estrutural, o vetor x é constituído por componentes de

deslocamentos e de velocidade, implicando em que, em um ponto fixo, as componentes de

velocidade e de aceleração são nulas, o que corresponde fisicamente a uma posição de

equilíbrio do sistema. Sendo assim, as soluções do sistema são conhecidas como soluções

estacionárias, soluções constantes ou soluções permanentes. Informações mais detalhadas

podem ser encontradas em Del Prado (2001) e Nayfeh e Balanchandran (1995).

3.5.4 SOLUÇÕES PERIÓDICAS

Uma solução periódica é uma solução dinâmica que possui uma freqüência básica

característica f. Uma solução ( )tXx= de um sistema é dita periódica com um período mínimo

ou básico T se ( ) ( )tTt XX =+ e ( ) ( )tt XX ≠+τ para T<<τ0 .

3.5.5 SEÇÕES DE POINCARÉ

A identificação das soluções periódicas, bem como o estudo de sua estabilidade,

são tarefas na maioria das vezes complexas, pois o sistema é contínuo e descrito por equações

diferenciais ordinárias.

Com o intuito de simplificar essa questão, transforma-se o sistema contínuo em

um modelo discreto representado por um mapeamento. Em geral, no estudo sobre estabilidade

de soluções periódicas, utiliza-se o mapeamento de Poincaré.

59

A seção de Poincaré é uma hiper-superfície imersa no espaço de estado que é

transversal ao fluxo definido pelo sistema de equações diferenciais. Para um sistema

autônomo, diz-se que a seção é transversal ao fluxo se:

( ) ( ) 0≠xFxnT , (3.49)

onde ( )xnT é o vetor normal à seção no ponto x.

No espaço do sistema não-autônomo, a condição de transversalidade é dada por

( )( ) ( ) 0, ≠ttT xFxn . (3.50)

No espaço n-dimensional, uma hiper-superfície é uma superfície de dimensão

menor que n e cada ponto na seção será descrito por n – 1 coordenadas. A Fig. 3.8 ilustra um

mapeamento de Poincaré típico para um sistema de duas equações diferenciais de primeira

ordem onde o plano inicial Σ é denominado plano de Poincaré ou plano fase. Os planos

paralelos ao plano fase Σ , em intervalos t∆ , são as seções de Poincaré.

x

y

t

P0 P3

P1

P2

seçã

o 1

seçã

o 2

seçã

o 3

Σ

∆t ∆t ∆t

x

y

t

P0 P3

P1

P2

seçã

o 1

seçã

o 2

seçã

o 3

Σ

∆t ∆t ∆t

Figura 3.8 – Planos de Poincaré.

O procedimento matemático que toma o conjunto de pontos seqüenciais, obtido

nas intersecções das seções de Poincaré com o fluxo fase do sistema, levando-o para o espaço

de condições iniciais, é dito mapeamento de Poincaré. Para problemas forçados com uma

excitação periódica, o intervalo regular de tempo t∆ é igual ao período da força excitadora

(MACHADO, 1993).

60

Quando um ponto é marcado exatamente sobre outro já existente durante o

processo de mapeamento, a órbita é dita fechada e este ponto é chamado de ponto fixo.

Necessariamente, uma solução periódica apresenta um número finito de pontos fixos, cada um

representando uma trajetória periódica com período igual ao número de intervalos de tempo

t∆ entre duas marcações coincidentes.

Se os novos pontos jamais coincidem com algum marcado anteriormente, a órbita

é dita aberta e o movimento é não-periódico. Nesta situação, o mapa apresenta infinitos

pontos geometricamente dispostos dentro de uma região limitada do espaço fase, o que indica

a presença de um movimento quase periódico e um atrator estranho ou caótico.

3.5.6 BIFURCAÇÃO DAS SOLUÇÕES DE EQUILÍBRIO

As equações que regem os sistemas dinâmicos são compostas por parâmetros de

controle. Em determinadas ocasiões, a variação gradual de certos parâmetros de controle pode

provocar, em um determinado ponto de instabilidade, denominado ponto crítico, uma

mudança brusca no comportamento global do sistema, denominada bifurcação. As

bifurcações, em geral, implicam na mudança no número e no tipo de soluções, ou seja, a partir

de um ponto de bifurcação podem surgir novas soluções diferentes ou similares, bem como

desaparecer uma família de soluções presentes antes do ponto de bifurcação. Ainda, uma

bifurcação pode mudar a estabilidade da solução de estável para instável ou vice-versa.

Quando se varia, de maneira quase incremental, um ou mais parâmetros de

controle ao longo de uma região que contém um ponto crítico, podem surgir bifurcações do

tipo contínua ou suave e do tipo descontínua ou catastrófica. As bifurcações são

caracterizadas pelo aparecimento de um novo atrator local que cresce suavemente após o

ponto de bifurcação. As bifurcações descontínuas são caracterizadas pelo salto repentino para

um novo grupo de soluções de equilíbrio, que podem estar até mesmo no infinito.

Segundo Nayfeh e Balachandran (1995), ao se variar um dos parâmetros de

controle de um ponto fixo estável de um sistema dinâmico contínuo, o ponto fixo pode perder

a estabilidade por uma das bifurcações seguintes: bifurcação nó-sela ou dobra cíclica,

bifurcação pitchfork ou simétrica, bifurcação transcrítica ou assimétrica e bifurcação de Hopf.

Essas bifurcações podem ser classificadas, ainda, em estáticas ou dinâmicas.

Uma bifurcação é dita estática quando nos pontos de bifurcação convergem

pontos fixos ou soluções de equilíbrio. A este grupo pertencem as bifurcações nó-sela,

61

simétrica e transcrítica. Na Fig. 3.9 estão ilustrados os cenários possíveis para as bifurcações

estáticas, onde as linhas tracejadas indicam soluções instáveis, enquanto as linhas contínuas

representam soluções estáveis.

A Fig. 3.9a ilustra uma bifurcação do tipo nó-sela, enquanto a Fig. 3.9b representa

uma bifurcação transcrítica. Já a Fig. 3.9c ilustra uma bifurcação simétrica do tipo super-

crítica, isto é, após o ponto crítico emergem dois ramos de soluções estáveis e outro ramo de

soluções instáveis. Finalmente a Fig. 3.9d mostra uma bifurcação simétrica sub-crítica, ou

seja, uma bifurcação simétrica instável, a qual apresenta localmente, de um lado do ponto de

bifurcação, dois ramos de pontos fixos instáveis mais um ramo de pontos fixos estáveis e, do

outro lado, um ramo de pontos fixos instáveis.

Des

loca

men

to

µ

Des

loca

men

to

µ (a)

Des

loca

men

to

µ

Des

loca

men

to

µ (b)

Des

loca

men

to

µ

Des

loca

men

to

µ (c)

Des

loca

men

to

µ

Des

loca

men

to

µ (d)

Figura 3.9 - Bifurcações estáticas das soluções de equilíbrio (DEL PRADO, 2001).

A bifurcação de um sistema autônomo é dita dinâmica quando convergem, para os

pontos de bifurcação, pontos fixos e soluções periódicas. A este grupo pertence à bifurcação

de Hopf. Na Fig. 3.10, são mostrados bifurcações de Hopf super-crítica e sub-crítica.

A Fig. 3.11 ilustra, para um sistema com um grau de liberdade, conforme

Thompson e Stewart (1993), como variam os autovalores à medida que se varia a rigidez e o

amortecimento efetivo da estrutura. Nesta figura, à medida que a rigidez estrutural diminui,

um foco estável se transforma em um nó estável, o qual se transforma em um ponto de sela

62

quando o maior autovalor cruza o eixo imaginário no ponto de bifurcação. Para sistemas

estruturais não-amortecidos, a transformação ocorre de um centro para um ponto de sela.

Desta forma, a perda de estabilidade estática do sistema acontece quando a rigidez efetiva

passa de positiva a negativa, representada pela seta horizontal na Fig. 3.11.

µµ

µ

Desloca

mento

Vel

ocid

ade

µ

Desloca

mento

Vel

ocid

ade

(a)

Des

loca

men

toD

eslo

cam

ento

µ

Desloca

mento

Vel

ocid

ade

µ

Desloca

mento

Vel

ocid

ade

(b)

Figura 3.10 - Bifurcação de Hopf: (a) super-crítica; (b) sub-crítica (DEL PRADO, 2001).

A Fig. 3.11 apresenta também a variação dos autovalores no fenômeno de

instabilidade dinâmica. Esse fenômeno ocorre em função do amortecimento efetivo. À medida

que diminui o amortecimento, passando-o de positivo para negativo, os autovalores

transformam-se de focos estáveis para focos instáveis, dando origem a um movimento

oscilatório de amplitude crescente. Este processo de perda de estabilidade é ilustrado pela seta

vertical na Fig. 3.11.

Um diagrama de bifurcação de um sistema dinâmico ilustra como os pontos fixos

do mapeamento de Poincaré variam de acordo com a mudança incremental de um

determinado parâmetro de controle. Os diagramas de bifurcação e os respectivos autovalores

de cada ponto fixo possibilitam detectar e caracterizar os tipos de bifurcações presentes em

um sistema dinâmico.

Neste trabalho foi utilizado o método numérico da força bruta obtido a partir de

Del Prado (2001). O método da força bruta é, sem dúvida, um dos métodos mais simples para

o traçado dos diagramas de bifurcação devido à sua simplicidade de programação e à sua

63

abrangência para encontrar vários tipos de soluções, como pontos fixos, soluções periódicas,

soluções quase periódicas e movimentos caóticos, que atestam a simplicidade do método.

Porém, o método apresenta uma convergência extremamente lenta para sistemas levemente

amortecidos, o que pode dificultar a obtenção das soluções permanentes do sistema, já que o

período de tempo pertinente à parte transiente é desconhecido.

I

R

Foco

I

R

Sela

Rigidez

Amortecimento

D > 0 D = 0

I

R

I

R

NóI

R

Nó

I

R

Foco

D < 0

I

R

Sela

I

R

I

R

Centro

I

R

Centro

I

R

Sela

I

R

I

R

Nó

I

R

Nó

Inst

abili

dade

Din

âmic

aInstabilidade Estática

I

R

Foco

I

R

Foco

I

R

Sela

I

R

Sela

Rigidez

Amortecimento

D > 0 D = 0

I

R

I

R

I

R

Nó

I

R

NóI

R

Nó

I

R

Nó

I

R

Foco

I

R

Foco

D < 0

I

R

Sela

I

R

Sela

I

R

I

R

I

R

Centro

I

R

Centro

I

R

Centro

I

R

Centro

I

R

Sela

I

R

Sela

I

R

I

R

I

R

Nó

I

R

Nó

I

R

Nó

I

R

Nó

Inst

abili

dade

Din

âmic

aIn

stab

ilida

de D

inâm

icaInstabilidade EstáticaInstabilidade Estática

Figura 3.11 - Estrutura das raízes para a equação 0=++ cxxbx &&& , onde b é o amortecimento

efetivo e c a rigidez efetiva. cbD 422 −= .

Basicamente, no método da força bruta, para se identificar o atrator associado a

um dado parâmetro de controle, escolhe-se um conjunto de condições iniciais e integra-se o

sistema por um tempo suficientemente longo para se chegar à resposta permanente, embora

não se tenha a garantia de que a integração, após o período de tempo especificado, vai

convergir para o atrator desejado. A utilização do método segue os seguintes passos:

Seja α o parâmetro de controle em questão, minα e maxα os limites para variação

do parâmetro de controle e K o número de incrementos. A partir destes dados, obtém-se o

incremento do parâmetro de controle para o traçado do diagrama de bifurcação, a saber:

64

( ) KkK

Kk ...,,1,1

1 minmaxmin =

−−

−±=αα

αα 3.51

Para a k-ésima iteração, o sistema de equações diferenciais é integrado

numericamente, sendo as primeiras Nt iterações assumidas como referentes à fase transiente e,

portanto, não consideradas. Nas Nss iterações seguintes assume-se que o sistema já se encontra

no estado permanente e armazenam-se as coordenadas do mapa de Poincaré associadas a cada

valor de α , o que permite traçar as diversas projeções do diagrama de bifurcação.

Cabe ressaltar que, em um diagrama assim obtido, não se pode diferenciar uma

resposta com um período muito longo de uma resposta quase periódica ou de uma resposta

caótica. Para isto, é necessário empregar outras técnicas de análise.

3.5.7 FRONTEIRAS DE INSTABILIDADE E ESCAPE DE SISTEMAS FORÇADOS

Para estruturas submetidas a carregamento harmônico, é importante realizar um

estudo detalhado da resposta do sistema no espaço dos parâmetros de controle para detectar

possíveis perdas de estabilidade. Assim, os algoritmos descritos a seguir visam identificar as

fronteiras de estabilidade associadas à força externa, o tempo de escape da bacia de atração

pré-flambagem e a periodicidade da resposta dinâmica.

Na Fig. 3.12 é apresentado um diagrama de bifurcação típico. Ao se variar um

determinado parâmetro de controle, um comportamento possível é representado pela reta de

estabilidade. Esse trecho ilustra que, para valores do parâmetro menores que o valor crítico,

denominado instabilidade paramétrica, a amplitude da resposta dinâmica do sistema decresce,

convergindo para solução fundamental.

O ponto de instabilidade paramétrica representa o valor do parâmetro de controle

a partir do qual o sistema começa a apresentar vibrações que podem ir aumentando até atingir

o escape pré-flambagem. Quando se fala em escape, entende-se que as condições do

parâmetro de controle provocam vibrações de grande amplitude. Nota-se na Fig. 3.12 que as

soluções do sistema sofrem mudanças no período até o escape para um atrator remoto.

Para a determinação das fronteiras de instabilidade paramétrica e de escape foi

utilizado o algoritmo apresentado por Del Prado (2001). Neste algoritmo, parte-se das

condições iniciais e de um valor para o parâmetro de freqüência Ω . Integra-se o sistema por

um longo período de tempo, avaliando se a vibração final da resposta é superior a um valor

pré-fixado.

65

Des

loca

men

to

Parâmetro de controle

Estabilidade

InstabilidadeParamétrica

Soluções periódicas e caóticasno vale pré-flambagem

Escape do vale potencial

Outras soluçõeslimitadas

Des

loca

men

to

Parâmetro de controle

Estabilidade

InstabilidadeParamétrica

Soluções periódicas e caóticasno vale pré-flambagem

Escape do vale potencial

Outras soluçõeslimitadas

Figura 3.12 – Diagrama de bifurcação típico.

Caso não haja o escape, incrementa-se o valor da amplitude de excitação até que

ocorra o mesmo. Quando isto acontecer, anula-se a amplitude da excitação e incrementa-se a

freqüência da excitação.

O valor que é pré-fixado para realizar a comparação com a amplitude de vibração

da resposta é determinado em função do tipo de fronteira a ser traçada. Caso seja a fronteira

de instabilidade paramétrica, o valor a ser considerado é nulo e, caso seja a fronteira de

escape, o valor a ser considerado depende do pré-carregamento estático que, por sua vez, gera

valores distintos em função do caminho pós-crítico da estrutura.

3.5.8 ANÁLISE GLOBAL DE SISTEMAS DINÂMICOS – DETERMINAÇÃO DAS BACIAS

DE ATRAÇÃO

A análise global de um sistema dinâmico tem como objetivo identificar, para uma

determinada região do espaço de estado, todas as soluções periódicas estáveis e instáveis do

sistema, obtendo-se, assim, a bacia de atração das soluções assintoticamente estáveis.

Para a determinação da bacia de atração das soluções periódicas estáveis, o tempo

considerado nos cálculos deve tender ao infinito. Isso é feito para eliminar a resposta

transiente de forma que a resposta do sistema se estabilize. Em termos numéricos essa medida

é impossível. Adota-se então um tempo múltiplo do período da força excitadora no qual a

66

bacia de atração é calculada. Esse procedimento implica em uma aproximação do que seria

real e é tão boa quanto maior for o tempo adotado.

No presente trabalho foi utilizado o método da força bruta para o mapeamento das

células, que é o método mais simples para se traçar as bacias de atração. Este método consiste

em discretizar o espaço onde se deseja calcular a bacia de atração em um conjunto de células

para os quais se tomam as propriedades do centro da célula como as propriedades da célula

como um todo. Para calcular a bacia de atração através do método da força bruta foi utilizado

o algoritmo apresentado por Del Prado (2001).

3.6 IMPLEMENTAÇÃO COMPUTACIONAL

Nos procedimentos computacionais adotados na implementação do esquema de

solução não-linear foi utilizado como base para este trabalho o programa de análise estrutural

desenvolvido em linguagem Fortran, por Silveira (1995) e Galvão (2000), e posteriormente

utilizado por Campos Filho (2004) em linguagem C++.

Basicamente, o programa é composto por cinco módulos principais:

a) No módulo um (entrada de dados) é feita a leitura dos dados necessários à

análise do problema. Inicialmente são lidos os parâmetros que definem o

modelo de elementos finitos adotado para o sistema estrutural e, em seguida,

são lidos os parâmetros que controlam a estratégia de solução não-linear

estática ou dinâmica;

b) No módulo dois (análise do problema não-linear estático) é realizada a análise

do sistema estrutural em estudo. Para cada incremento de carga, é resolvido o

sistema definido na Eq. 2.50;

c) No módulo três (freqüências naturais e modos de vibração), a solução do

problema de autovalor apresentado na Eq. 2.71, é obtida;

d) No módulo quatro (análise do problema dinâmico linear ou não-linear), é

encontrada a resposta no tempo do sistema estrutural, resolvendo a equação de

equilíbrio dinâmico, dada na Eq. 2.68;

e) No módulo cinco (saída de resultados) é feita a impressão dos resultados em

arquivos de saída utilizados para a impressão de um relatório da análise

efetuada, para a impressão da trajetória de equilíbrio e das deformadas da

estrutura ou das freqüências naturais, modos de vibração e resposta no tempo.

67

O procedimento inicial a ser realizado pelo programa principal é a abertura dos

arquivos de entrada e saída de dados, bem como, a leitura dos dados de entrada. Esses dados

definem a geometria do modelo estrutural, o número de pontos nodais, o número de

elementos, condições de contorno, carregamentos externos, as propriedades físicas e

geométricas das seções dos elementos que compõem a estrutura.

Caso a opção seja pela análise não-linear estática, chama-se a sub-rotina destinada

à leitura dos dados necessários a mesma, como: a estratégia de solução (Newton-Raphson

padrão ou modificado), o valor inicial do parâmetro de carga, o número de incrementos, o

critério de convergência (força, deslocamento ou ambos), o número máximo de iterações por

incremento, entre outros.

Em seguida, têm-se a sub-rotina destinada à montagem do vetor de cargas.

Quando o sistema proposto na análise possua cabos pré-tensionados, este processo será

efetuado em dois estágios. No primeiro estágio, as forças induzidas no pórtico pelo pré-

tensionamento dos cabos são determinadas e a condição de equilíbrio é verificada. No

segundo estágio, o carregamento externo é incrementado às forças encontradas anteriormente.

Montado o vetor de cargas, inicia-se o processo incremental-iterativo de solução, que é

resumido a seguir:

a) Chamada da rotina responsável pela montagem da matriz de rigidez e

obtenção dos deslocamentos nodais tangenciais;

b) Chamada da rotina responsável pela determinação da solução predita ∆λ0, e

∆u0 de acordo com a estratégia de incremento de carga escolhida;

c) Correção da solução predita pelo processo iterativo através de rotina

específica;

d) Correção dos parâmetros necessários ao próximo incremento de carga;

e) Se o número de passos de carga é menor do que o desejado recomeça-se o

processo.

Os resultados da análise não-linear estática são apresentados em arquivos de

saída, contendo: dados necessários à impressão das trajetórias de equilíbrio, dados para

visualização das configurações deformadas do sistema estrutural e o arquivo de relatório

contendo todas as informações do processo de solução. A Fig. 3.13 mostra o fluxograma do

programa principal.

68

Figura 3.13 – Fluxograma do programa principal para a análise não-linear.

Análise não-linear dinâmica

Análise linear dinâmica

Análise dinâmica Análise não-linear estática

Não

Sim

INÍCIO

Leitura dos dados da estrutura

Leitura dos dados para análise não-linear

geométrica

Montagem do vetor de cargas

Configuração inicial: tu e tλ

Matriz de Rigidez: K

Solução predita: ∆u0 e ∆λ0

Vetor de forças internas: Fi Ciclo iterativo

itera ≤ n° máximo de iterações

Vetor de forças residuais: g = λ Fr – Fi(u)

ζ ≤ ζadm

Cálculo de: δuk e δλk

Atualizam-se as variáveis incrementais totais

Novo incremento Ciclo incremental-iterativo

inc ≤ n° máximo de incrementos

Montagem do vetor de cargas

Matriz de Rigidez: K Matriz de massa: M

Método de Choleski Decomposição da matriz M

M = ST S

Obtenção da matriz Q = S-1

Obtenção da matriz A = QT K Q

Método de Jacobi Autovalores na diagonal da matriz A Autovetores nas colunas da matriz B

Obtenção dos autovetores

Arquivo de Saída

FIM DO PROCESSAMENTO

Leitura dos dados para integração direta

Leitura dos dados para integração não-linear Newmark linear

Newmark não-linear Arquivo de Saída

FIM DO PROCESSAMENTO Arquivo de Saída FIM DO

PROCESSAMENTO

69

Se a opção for pela análise dinâmica, chama-se a sub-rotina destinada à montagem

do vetor de cargas e, em seguida, inicia-se o cálculo da freqüência natural e dos modos de

vibração da estrutura, a saber:

a) Chama-se a rotina para a montagem da atriz de rigidez linear geométrica da

estrutura e a rotina responsável pela montagem da matriz de massa

consistente;

b) Chama-se a rotina destinada ao cálculo dos autovalores e autovetores pelo

método de Jacobi;

c) Chama-se a rotina que grava no arquivo de saída as freqüências de vibração e

os modos de vibração normalizados.

Para a análise dinâmica linear, chama-se em seguida, a sub-rotina destinada à

leitura dos dados necessários à integração direta, tais como: tipo de carregamento (harmônico

ou não), função de excitação (seno ou co-seno), a freqüência de excitação, entre outros. Na

análise dinâmica não-linear, após o cálculo das freqüências naturais e modos de vibração,

chama-se a sub-rotina destinada à leitura dos dados necessários à integração não-linear, como

por exemplo, as constantes de amortecimento proporcional de rigidez e de massa.

Ao final desta etapa, os resultados da análise linear ou não-linear dinâmica são

apresentados em arquivos de saída, contendo: as freqüências naturais, os modos de vibração

normalizados, os pontos para se traçar a resposta no tempo, dados para visualização dos

modos de vibração e o arquivo de relatório contendo todas as informações do processo de

solução.

70

4 RESULTADOS NUMÉRICOS

4.1 INTRODUÇÃO

Neste capítulo são apresentados exemplos numéricos para validar as formulações

matemáticas, para avaliar os caminhos pós-críticos, as freqüências naturais e os modos de

vibração de colunas perfeitas e imperfeitas. Avalia-se, também, a instabilidade dinâmica sob

carga súbita e harmônica, os mecanismos de escape e as bacias de atração de alguns

exemplos.

4.2 ANÁLISE ESTÁTICA NÃO-LINEAR VIA MÉTODO DOS ELEMENTOS FINITOS

4.2.1 VALIDAÇÃO DO MÉTODO

O pórtico de Lee é um pórtico em forma de L bi-rotulado, constituído por duas

barras e submetido a uma força vertical P situada a 51 do nó de ligação das barras, conforme

mostrado na Fig. 4.1. Por tratar-se de um sistema estrutural cuja trajetória de equilíbrio

apresenta comportamento não-linear complexo, esta estrutura é amplamente utilizada na

literatura para validar formulações matemáticas e estratégias de solução não-linear. Na Tabela

4.1 apresentam-se as características físicas e geométricas da estrutura analisada.

Neste trabalho, analisou-se, primeiramente, a forma clássica do pórtico de Lee,

isto é, o pórtico sem associação com cabos, modelando-o com 6 e 20 elementos. Avaliou-se o

comportamento da estrutura utilizando o método de Newton-Raphson associado à estratégia

do comprimento de arco com o primeiro incremento 0,10 =∆λ e tolerância 0401x1 −=ε .

71

Tabela 4.1 – Propriedades físicas e geométricas do Pórtico de Lee e do cabo.

Propriedades Pórtico de Lee Cabo

Área (A) 6,00 cm² 6,00 x 10 - 05 cm²

Momento de Inércia ( xI ) 2,00 cm4 -

Módulo de Elasticidade (E) 7,20 kN/cm² 7,20 kN/cm²

Ο

P

u

v

120

cm

24 cm 96 cm

(a)

120

cm

24 cm

v

P

Ο u

30 c

m

96 cm

(b)

Figura 4.1 – Pórtico de Lee: (a) Sem cabo; (b) Com cabo.

A Fig. 4.2 apresenta as trajetórias de equilíbrio para o pórtico de Lee,

relacionando a carga vertical P aos deslocamentos verticais (v) e horizontais (u) do ponto de

aplicação da carga vertical.

0 20 40 60 80 100

v [cm]

-1.5

-1

-0.5

0

0.5

1

1.5

2

2.5

P [kgf]

06 Elementos20 ElementosSchweizerhof e Wriggers (1986)

(a)

0 20 40 60 80 100

u [cm]

-1.5

-1

-0.5

0

0.5

1

1.5

2

2.5

P [kgf]

06 Elementos20 ElementosSchweizerhof e Wriggers (1986)

(b)

Figura 4.2 – Trajetória de equilíbrio para Pórtico de Lee: (a) P x v; (b) P x u.

72

Observa-se que os resultados são significativamente influenciados pela

modelagem da estrutura. Percebe-se que uma modelagem com poucos elementos não é

suficiente para determinar corretamente o comportamento da estrutura. No modelo estrutural

com 20 elementos os resultados são coincidentes com os obtidos por Schweizerhof e

Wriggers (1986) que também utilizaram elementos finitos.

Em seguida foi acrescentado um elemento de cabo ao pórtico de Lee, mantendo a

modelagem do pórtico com 20 elementos idênticos (10 para cada barra).

A Fig. 4.3 apresenta as trajetórias de equilíbrio para o pórtico de Lee com e sem

cabo, relacionando a carga vertical P aos deslocamentos verticais (v) e horizontais (u) do

ponto de aplicação da carga vertical.

0 20 40 60 80 100

v [cm]

-1

-0.5

0

0.5

1

1.5

2

2.5

P [kgf]

Pórtico Lee EstaiadoPórtico Lee

Pontos da Deformada

AB

C

D

E

F

A'B'

C'

D' E'

F'

(a)

0 20 40 60 80 100

u [cm]

-1

-0.5

0

0.5

1

1.5

2

2.5

P [kgf]

Pórtico Lee EstaiadoPórtico Lee

Pontos da Deformada

A

B

C

D E

F

A'B'

C' D'

E'

F'

(b)

Figura 4.3 – Trajetória de equilíbrio para o Pórtico de Lee estaiado: (a) P x v; (b) P x u.

É possível observar que, no início, a trajetória de equilíbrio de ambos os modelos

seguem o mesmo caminho, entretanto, no sistema estrutural com cabo, o pórtico atinge um

ponto limite de carga superior ao do modelo do pórtico sem o cabo. Percebe-se que, no

modelo do pórtico com cabo, a perda de rigidez, além de ser menor, ocorre de forma menos

brusca, obtendo também, ganho de rigidez mais rápidamente. Verifica-se, ainda, o fenômeno

snap-back nas trajetórias de equilíbrio de ambas as estruturas.

Pode-se afirmar que, nos estágios iniciais e avançados da trajetória de equilíbrio, a

influência do cabo no comportamento não-linear é pequena visto que, as curvas de ambos os

modelos, são bastante semelhantes. A Fig. 4.4 apresenta as configurações deformadas da

estrutura para o modelo sem cabo (a) e com cabo (b).

73

AB C D

E

F

(a)

A' B' C'

D'

E'

F'

(b)

Figura 4.4 – Deformada do Pórtico de Lee: (a) Sem cabo; (b) Com cabo.

A Fig. 4.5 mostra uma coluna com seção circular oca engastada em sua base e

livre na extremidade superior. Na parte superior, age uma força vertical conservativa de

compressão P. Como o movimento do topo da coluna não é prescrito, é necessário provocar

uma pequena perturbação capaz de iniciar seu deslocamento lateral; essa perturbação, na

prática, decorre de imperfeições geométricas da estrutura ou da dificuldade em se obter um

carregamento axial coincidente com o centro de gravidade da seção transversal. Dessa forma,

foi admitida uma excentricidade equivalente a 0,1% do comprimento da coluna, que produz

um momento fletor M, transversal ao plano de flambagem da peça, aplicado no topo da

coluna. A estrutura foi modelada com uma malha uniforme de 10 elementos de pórtico. As

propriedades físicas e geométricas da coluna circular podem ser vistas na Tabela 4.2.

Tabela 4.2 – Propriedades físicas e geométricas da coluna circular.

Propriedades Coluna Circular Oca

Diâmetro Interno ( iφ ) 0,475 m

Diâmetro Externo ( eφ ) 0,500 m

Área (A) 1,914 x 10 -02 m²

Momento de Inércia ( xI ) 5,691 x 10 -04 m4

Módulo de Elasticidade (E) 1,18 x 10 +08 kN/m²

74

Figura 4.5 – Coluna engastada.

Na determinação das trajetórias de equilíbrio da estruturas, utilizou-se o método

de Newton-Raphson associado à estratégia de iteração baseada no método do controle dos

deslocamentos generalizados com tolerância 01,0=ε .

A Fig. 4.6 apresenta a trajetória de equilíbrio para a coluna engastada,

relacionando o parâmetro de carga vertical EIPL2 aos deslocamentos verticais ( Lv ) e

horizontais ( Lu ) do ponto de aplicação da carga vertical. Na Fig. 4.7 são mostradas suas

configurações deformadas.

0 0.4 0.8 1.2 1.6

v/L

0

2

4

6

8

PL²/EI

10 Elementos

Pontos da Deformada

AB

C

D

(a)

-1 -0.8 -0.6 -0.4 -0.2 0

u/L

0

2

4

6

8

PL²/EI

10 Elementos

Pontos da Deformada

AB

C

D

(b)

Figura 4.6 – Trajetória de equilíbrio na extremidade livre (topo) da coluna: (a) P x v; (b) P x u.

P

vθ u

M = 1 N/m

L = 100 m

75

Observa-se que as trajetórias de equilíbrio para a coluna engastada-livre

apresentam um contínuo ganho de rigidez, sendo mais acentuado na ocorrência de grandes

deslocamentos. Além disso, apresenta também ponto de bifurcação estável, com parâmetro

43,22 =EIPL . A carga crítica de Euler para coluna engastada-livre é de 22 4LEIPe π=

(BAZANT e CEDOLIN, 1991), isto é, 47,2422 == πEILPe , ou seja, 1,62% maior que a

obtida neste trabalho.

A

B

C

D

Figura 4.7 – Deformada da coluna engastada-livre.

Em seguida, como mostra a Fig. 4.8, foram acrescentados à coluna dois cabos

com inclinação =α 60°, tensionamento 10=T kN, diâmetro 018,0=φ m e módulo de

elasticidade longitudinal E = 1 x 10 +08 kN/m². Cada cabo foi discretizado com um elemento,

mantendo-se a modelagem anterior da coluna.

Figura 4.8 – Coluna estaiada com dois cabos a 60°.

L=100 m

α α

M = 1 N/m

θv

u

P

76

A Fig. 4.9 mostra as trajetórias de equilíbrio da coluna estaiada, para os

deslocamentos vertical ( Lv ) e horizontal ( Lu ) da extremidade superior da coluna. Para fins

de comparação, apresentam-se, as trajetórias de equilíbrio da coluna sem os cabos.

0 0.4 0.8 1.2

v/L

0

5

10

15

20

PL²/EI

Coluna Com CabosColuna Sem Cabos

Pontos da Deformada

A

B

C D E

(a)

-1 -0.8 -0.6 -0.4 -0.2 0

u/L

0

5

10

15

20

PL²/EI

Coluna Com CabosColuna Sem Cabos

Pontos da Deformada

A

B

C

D E

(b)

Figura 4.9 – Trajetória de equilíbrio para a coluna engastada estaiada: (a) P x v; (b) P x u.

Nota-se que o valor limite da carga para a coluna com os cabos foi,

aproximadamente, 18 vezes maior que a carga crítica da coluna sem cabo. Contudo, na coluna

sem cabo, o caminho pós-crítico apresenta ganho de rigidez, ao passo que, na coluna com

cabo, após a flambagem, há perda de rigidez para posterior ganho de rigidez.

Na Fig. 4.10 mostra-se a evolução da deformada da coluna estaiada com base

engastada à medida que se incrementa a carga. Nesta deformada, os elementos de cabo que

estão submetidos à compressão foram eliminados, pois os cabos não apresentam rigidez à

compressão.

AB

C

D

E

Figura 4.10 – Deformada da coluna estaiada.

77

4.2.2 EFEITO DAS IMPERFEIÇÕES GEOMÉTRICAS INICIAIS

Para avaliar os efeitos de uma imperfeição geométrica inicial no pórtico de Lee,

considera-se que suas barras, como mostrado na Fig. 4.11, possuem, antes de ser solicitada,

uma ligeira curvatura, descrita por 21 pares de coordenadas apresentados na Tabela 4.3.

Tabela 4.3 – Coordenadas que definem a geometria inicial do pórtico de Lee da Fig. 4.11.

Coordenadas Coordenadas Coordenadas

Nó x y Nó x y Nó x y

1 0,000 0,000 8 6,476 83,003 15 48,059 122,359

2 2,910 11,638 9 4,471 94,830 16 60,058 122,445

3 5,488 23,354 10 2,282 106,625 17 72,057 122,272

4 7,441 35,190 11 0,253 118,448 18 84,050 121,888

5 8,554 47,134 12 12,145 120,052 19 96,038 121,347

6 8,724 59,129 13 24,089 121,209 20 108,021 120,700

7 7,979 71,101 14 36,065 121,963 21 120,000 120,000

Logo, na modelagem, foram mantidos os 20 elementos de pórtico e foram

conservadas as propriedades físicas e geométricas dadas pela Tabela 4.1. Manteve-se,

também, a estratégia de iteração do tipo comprimento de arco com o primeiro incremento

1,00 =∆λ e tolerância de 0410x1 −=ε .

Figura 4.11 – Pórtico de Lee com imperfeições iniciais.

P

uv

24 cm 96 cm

θ

120 cm

59,129 cm

60,058 cm

2,445cm

8,724 cm

1

6

11 12

16

21

78

A Fig. 4.12 apresenta as trajetórias de equilíbrio para o pórtico de Lee com

imperfeições iniciais, relacionando a carga vertical P aos deslocamentos verticais (v) e

horizontais (u) do ponto de aplicação da carga vertical.

0 20 40 60 80 100

v [cm]

-1

-0.5

0

0.5

1

1.5

2

2.5

P [kgf]

Pórtico de Lee ImperfeitoPórtico de Lee

Pontos da Deformada

A

B

C

D

E

(a)

0 20 40 60 80 100

u [cm]

-1

-0.5

0

0.5

1

1.5

2

2.5

P [kgf]

Pórtico de Lee ImperfeitoPórtico de Lee

Pontos da Deformada

A

B

C D

E

(b)

Figura 4.12 – Trajetória de equilíbrio para Pórtico de Lee com imperfeições geométricas

iniciais: (a) P x v; (b) P x u.

Observa-se neste caso que a trajetória de equilíbrio do pórtico foi fortemente

influenciada pela imperfeição geométrica inicial, diferenciando-se completamente da

trajetória de equilíbrio para o pórtico perfeito. No estágio inicial do caminho pós-crítico, o

pórtico imperfeito apresenta maior ganho de rigidez, mas atinge um ponto limite com menor

carga. Após o ponto limite, há perda de rigidez seguida do fenômeno snap-back mais

acentuado que no pórtico perfeito, sendo que o salto foi menor nos deslocamentos verticais do

que nos horizontais.

Assim como no tratamento dado ao pórtico de Lee perfeito, acrescenta-se no

pórtico de Lee imperfeito um elemento de cabo, a fim de se analisar também o

comportamento não-linear do pórtico de Lee estaiado com imperfeições geométricas iniciais,

representado na Fig. 4.13. A Fig. 4.14 apresenta as trajetórias de equilíbrio para este caso,

relacionando a carga vertical P aos deslocamentos verticais (v) e horizontais (u) do ponto de

aplicação da carga vertical.

Observa-se que o ponto limite de carga no modelo do pórtico com cabo é

praticamente o mesmo do pórtico sem cabo. Entretanto, a carga mínima atingida no modelo

do pórtico com cabo é, em módulo, menor do que no modelo sem cabo, indicando assim uma

perda de rigidez também menor. Posteriormente, no modelo do pórtico com cabo ocorre um

79

ganho de rigidez mais rápido em comparação com o modelo sem cabo. Verifica-se ainda que,

nos estágios iniciais e avançados da trajetória de equilíbrio, a influência do cabo no

comportamento não-linear é pequena, devido à similaridade das curvas. Na Fig. 4.15 são

mostradas as configurações deformadas da estrutura.

Figura 4.13 – Pórtico de Lee estaiado com imperfeições iniciais.

0 20 40 60 80 100

v [cm]

-1

-0.5

0

0.5

1

1.5

2

2.5

P [kgf]

Pórtico de Lee Estaiado ImperfeitoPórtico de Lee Imperfeito

Pontos da Deformada

A

B

C

D

E

(a)

0 20 40 60 80 100

u [cm]

-1

-0.5

0

0.5

1

1.5

2

2.5

P [kgf]

Pórtico de Lee Estaiado ImperfeitoPórtico de Lee Imperfeito

Pontos da Deformada

A

B

C

D E

(b)

Figura 4.14 – Trajetória de equilíbrio para Pórtico de Lee estaiado com imperfeições

geométricas iniciais: (a) P x v; (b) P x u.

Em um segundo exemplo, visando mostrar a influência das imperfeições iniciais,

procedeu-se a análise da coluna estaiada, submetendo-a agora, a uma deformação inicial

originada pela aplicação de tensionamento diferenciado nos cabos, sendo que, no cabo

esquerdo considerou-se 91 =T kN e no direito 102 =T kN.

P

uv

96 cm

θ

120 cm

59,129 cm

60,058

2,445cm

8,724 cm

24 cm

30 cm

1

6

11 12

16

21

80

A

B

CD

E

(a)

A

B

CD E

(b)

Figura 4.15 – Deformada do Pórtico de Lee estaiado com imperfeições geométricas iniciais:

(a) Sem cabo; (b) Com cabo.

Na prática, a diferença no tensionamento dos cabos acarreta no acomodamento da

estrutura, ficando esta, com uma imperfeição geométrica inicial que, neste exemplo, foi o da

coluna estaiada mostrada na Fig. 4.16, a qual possui um deslocamento na extremidade livre de

2,643 cm para a direita e um abaulamento de 5,903 cm para a esquerda, medido na metade da

altura.

As propriedades físicas e geométricas da coluna estaiada imperfeita são as

indicadas na Tabela 4.2. Os cabos possuem inclinação =α 60°, diâmetro 018,0=φ m e

módulo E = 1 x 10 +08 kN/m².

Figura 4.16 – Coluna estaiada imperfeita com dois cabos a 60°.

100 m

α α

M = 1 N/m

θv

u

P

50 m

0,05903 m

0,02643 m

81

A Fig. 4.17 mostra as trajetórias de equilíbrio da coluna estaiada sujeita à

tensionamento diferenciado, relacionando o parâmetro de carga vertical EIPL2 aos

deslocamentos verticais ( Lv ) e horizontais ( Lu ) da extremidade superior da coluna. Para

fins de comparação, apresenta-se também, nesta ilustração, a trajetória de equilíbrio da coluna

estaiada perfeita.

0 0.4 0.8 1.2

v/L

0

4

8

12

16

20

PL²/EI

T1 = T2T1 < T2

Pontos da Deformada

A

BC D E

(a)

-0.6 -0.4 -0.2 0 0.2 0.4 0.6

u/L

0

4

8

12

16

20

PL²/EI

T1 = T2T1 < T2

Pontos da Deformada

AB

C D E

F

(b)

Figura 4.17 – Trajetória de equilíbrio para coluna estaiada com geometria inicial deformada e

tensionamento diferenciado: (a) P x v; (b) P x u.

Comparando o caso perfeito com o imperfeito, observa-se que as trajetórias de

equilíbrio são complemente diferentes. O ganho de rigidez até o ponto limite é drasticamente

reduzido na coluna estaiada imperfeita e, posteriormente, verifica-se uma perda de rigidez

menos acentuada ao longo dos deslocamentos. Em ambos os casos a estrutura volta a ter

ganho de rigidez no trecho final das trajetórias que, na coluna perfeita, ocorre para valores de

deslocamento maiores. Pode-se observar ainda que, para os valores de tensionamento

aplicados, há uma inversão no sentido dos deslocamentos, como pode ser visto na Fig. 4.18.

Avançando no estudo das imperfeições, na torre estaiada da Fig. 4.19, foram

considerados também, os efeitos das imperfeições geométricas iniciais impostas por

tensionamentos diferenciados nos estais. Na Fig. 4.20 são apresentas as trajetórias de

equilíbrio da coluna em relação a seus deslocamentos verticais ( Lv ) e horizontais ( Lu ) no

topo da coluna. A estrutura foi solicitada por T1 = 10 kN e T2 = 10 kN, correspondendo ao

caso perfeito, e também para as combinações T1 = 10 kN e T2 = 9,8 kN, T1 = 10 kN e

T2 = 9,6 kN e T1 = 10 kN e T2 = 8,0 kN.

82

A

B

C

D

E

Figura 4.18 – Deformada coluna estaiada com tensionamento diferenciado.

Figura 4.19 – Coluna estaiada com tensionamentos diferentes.

0 0.4 0.8 1.2 1.6

v/L

0

4

8

12

16

20

PL²/EI

T2 = 10 kNT2 = 9,8 kNT2 = 9,6 kNT2 = 8,0 kN

(a)

-0.6 -0.4 -0.2 0

u/L

0

4

8

12

16

20

PL²/EI

T2 = 10 kNT2 = 9,8 kNT2 = 9,6 kNT2 = 8,0 kN

(b)

Figura 4.20 – Trajetória de equilíbrio para coluna estaiada: (a) P x v; (b) P x u.

100 m

α α

M = 1 N/m

θv

u

P

T1 T2

83

Verifica-se que, para T2 = 9,8 kN, a estrutura apresentou carga crítica maior que a

do caso perfeito; para T2 = 9,6 kN a carga crítica foi praticamente igual à do caso perfeito e

para T2 = 8,0 kN, foi menor. Na Fig. 4.21 é apresentada a curva de sensibilidade a

imperfeições geométricas iniciais da torre, considerando T1 = 10 kN e variando T2.

1 0.96 0.92 0.88 0.84 0.8

T2/T1

14

15

16

17

18

19

20

PL²/EI

Figura 4.21 – Curva de sensibilidade a imperfeições geométricas inicias para T1 = 10 kN e T2

variável.

Observa-se que no intervalo de 10 ≥ T2 ≥ 9,7 kN a estrutura apresentou ganho de

rigidez; para T2 = 9,6 kN, a rigidez aproxima-se do valor encontrado no caso perfeito; no

intervalo de 9,6 ≥ T2 ≥ 9,5 kN a estrutura perdeu drasticamente sua rigidez e a partir deste

valor até T2 = 8,0 kN a curva apresenta um patamar constante.

Procedeu-se ainda com o estudo das imperfeições em colunas estaiadas, avaliando

a influência de deslocamentos horizontais δ, aplicado no topo da estrutura, representado no

modelo esquemático da Fig. 4.22.

Figura 4.22 – Coluna estaiada imperfeita 0 ≤ δ ≤ 1,0 m.

100 m

α

M = 1 N/m

θv

u

P

α

δ

84

Na Fig. 4.23 são mostradas as trajetórias de equilíbrio da coluna para os

deslocamentos verticais ( Lv ) e horizontais ( Lu ) do topo da coluna, considerando δ = 0,0,

δ = 0,2 m e δ = 1,0 m.

Novamente, as propriedades físicas e geométricas da coluna são as indicadas na

Tabela 4.2. Os cabos possuem inclinação =α 60°, tensionamento constante 10=T kN,

diâmetro 018,0=φ m e módulo de elasticidade E = 1 x 10 +08 kN/m².

0 0.4 0.8 1.2 1.6

v/L

0

4

8

12

16

20

PL²/EI

δ = 0,0 mδ = 0,2 mδ = 1,0 m

(a)

-0.6 -0.5 -0.4 -0.3 -0.2 -0.1 0

u/L

0

4

8

12

16

20

PL²/EI

δ = 0,0 mδ = 0,2 mδ = 1,0 m

(b)

Figura 4.23 – Trajetória de equilíbrio para coluna estaiada: (a) P x v; (b) P x u.

Observa-se que as curvas para os casos imperfeitos são muito próximas à do caso

perfeito com variação máxima de 2% da carga crítica perfeita, apresentando a tendência de

aumentar com o aumento das imperfeições geométricas iniciais aplicadas, como mostra a

curva de sensibilidade da Fig. 4.24.

0 0.003 0.006 0.009

δ/L

18.7

18.8

18.9

19

19.1

19.2

19.3

PL²/EI

Figura 4.24 – Curva de sensibilidade a imperfeições geométricas inicias para 0 ≤ δ ≤ 0,9 m.

85

A curva de sensibilidade relaciona a carga crítica de flambagem com vários casos

de deslocamentos horizontais, δ ≤ 0,9 m. Verifica-se que, para as grandezas de deslocamentos

horizontais estudados, a estrutura apresenta ganho de rigidez quando comparada ao caso

perfeito.

4.3 ANÁLISE DINÂMICA LINEAR VIA MÉTODO DOS ELEMENTOS FINITOS

Nesta seção são determinadas as freqüências naturais e os modos de vibração de

colunas estaiadas perfeitas e imperfeitas.


Inicialmente, apresenta-se o comportamento da coluna engastada de seção

constante e sem força axial mostrada na Fig. 4.25. Essa coluna foi estudada anteriormente por

Orlando (2006). As propriedades físicas e geométricas, necessárias à análise, podem ser vistas

na Tabela 4.4.

Figura 4.25 – Coluna engastada de seção constante sem força axial.

Tabela 4.4 – Propriedades físicas e geométricas da coluna de seção constante e descarregada.

Propriedades Coluna

Área (A) 1,91 x 10 -02 m²

Momento de Inércia ( xI ) 5,69 x 10 -04 m4

Densidade ( ρ ) 150,25 kg/m

Módulo de Elasticidade (E) 1,18 x 10 +08 kN/m²

1m

v θ u

86

Na Tabela 4.5 são apresentadas as três primeiras freqüências encontradas por

Orlando (2006) e as obtidas neste trabalho. Como pode ser observado, os resultados obtidos

são bastante satisfatórios quando comparados aos da literatura. Na Fig. 4.26 são mostrados os

três primeiros modos de vibração da coluna.

Tabela 4.5 – Comparação das freqüências naturais.

Freqüências Naturais Orlando (2006) Presente Trabalho

rad/seg rad/seg

1ª 7,43328 x 10 +02 7,43320 x 10 +02

2ª 4,65835 x 10 +03 4,65831 x 10 +03

3ª 1,30435 x 10 +04 1,30436 x 10 +04

(a)

(b)

(c)

Figura 4.26 – Modos de Vibração Normalizados: (a) 1° Modo; (b) 2° Modo; (c) 3° Modo.

A Fig. 4.27 mostra dois modelos de colunas estaiadas sujeitas a dois tipos de

carregamento: uma carga axial P no topo da coluna e uma carga axial p a 7 m de sua base. Na

Fig. 4.27a cada cabo possui inclinação α e está fixado à barra a uma distância de 8 m do

apoio. Na Fig. 4.27b acrescentam-se dois cabos, ambos fixados à barra a uma distância de 6 m

do apoio. Consideram-se nestas colunas as propriedades físicas e geométricas indicadas na

Tabela 4.4. Os cabos têm tensionamento constante 10=T kN, diâmetro 018,0=φ m e

módulo de elasticidade E = 1 x 10 +08 kN/m².

A Fig. 4.28 exibe a variação da freqüência natural em função da inclinação dos

cabos. Para o modelo com dois cabos, a maior freqüência natural obtida foi para α = 40°,

enquanto que, para o modelo com quatro cabos, a maior freqüência natural foi para α = 45°.

Observa-se ainda que, para a coluna com quatro cabos, a freqüência natural foi

aproximadamente 30% maior.

87

(a)

(b)

Figura 4.27 – Coluna Estaiada: (a) Dois Cabos; (b) Quatro Cabos.

0 10 20 30 40 50 60 70 80 90

α [°]

0

0.4

0.8

1.2

1.6

2

2.4

2.8

3.2

3.6

4

f [H

z]

02 Cabos04 Cabos

Figura 4.28 – Variação da freqüência natural em função do ângulo α.

4.3.2 EFEITO DAS IMPERFEIÇÕES GEOMÉTRICAS INICIAIS NAS FREQÜÊNCIAS

NATURAIS

Na Tabela 4.6 são apresentadas as freqüências naturais para os casos imperfeitos

estudados na Fig. 4.19 e na Fig. 4.29 são mostrados os primeiros modos de vibração para os

respectivos casos. Pode-se observar que os valores das freqüências naturais são bastante

semelhantes. Quanto aos modos de vibração, observa-se a semelhança nas configurações do

primeiro modo de vibração das Figuras 4.29a e 4.29b. Já na Fig. 4.29c, pode-se observar uma

configuração bastante complexa e na Fig. 4.29d vê-se que o primeiro modo de vibração

possui as menores amplitudes de deslocamento quando comparado os demais.

θ

10 m

α α

M = 1 N/m

v u

P

p

6 m

θ

10 m

α α

M = 1 N/m

v u

P

p 8

m

7 m

88

Tabela 4.6 – Freqüência natural.

T2 (kN) rad/seg T2 (kN) rad/seg

10 0,35 9,6 0,35

9,8 0,36 8,0 0,37

(a) 1° Modo.

(b) 1° Modo.

(c) 1° Modo.

(d) 1° Modo.

Figura 4.29 – Modos de vibração normalizados para T1 = 10 kN e: (a) T2 = 10 kN; (b) T2 =

9,8 kN; (c) T2 = 9,6 kN; (d) T2 = 8,0 kN.

Na Fig. 4.30 é mostrado o primeiro modo de vibração da coluna perfeita e

imperfeita com δ = 0,2 m e δ = 1,0 m, mostrada na Fig. 4.22. Em todos os casos, a freqüência

natural foi de 0,35 rad/seg.

Observa-se que quase não há alteração na forma entre o primeiro modo de

vibração da Fig. 4.30a e da Fig.4.30b. Na Fig. 4.30c pode-se observar que o modo tem

amplitude máxima no topo da coluna.

Em resumo, pode-se afirmar que as imperfeições geométricas iniciais afetam a

forma do modo de vibração das colunas estudadas, porém sua influência depende do tipo e da

intensidade da imperfeição.

89

(a) 1° Modo.

(b) 1° Modo.

(c) 1° Modo.

Figura 4.30 – Modos de vibração normalizados: (a) δ = 0,0 m; (b) δ = 0,2 m; (c) δ = 1,0 m.

4.4 ANÁLISE DINÂMICA NÃO-LINEAR VIA MÉTODO DOS ELEMENTOS FINITOS

Nesta seção é analisado o comportamento da coluna estaiada quando submetida à

ação de uma carga axial súbita ou harmônica. É importante recordar que foi adotada para a

análise a razão de amortecimento ξ = 1% (Item 2.11).

4.4.1 INSTABILIDADE DINÂMICA PARA CARGA SÚBITA

Na Fig. 4.32 mostram-se seis respostas no tempo do topo da coluna engastada da

Fig. 4.5 (a se ver novamente na Fig. 4.31), quando submetida a valores incrementais da carga

súbita P. A carga crítica da coluna é de Pcr = 16,3517 kN.

Figura 4.31 – Coluna engastada com carregamento súbito P.

Observa-se que, para P ≤ 0,80 Pcr (Figuras 4.32a, 4.32b, 4.32c e 4.32d), a coluna

apresenta vibrações de pequena amplitude, com a resposta permanente tendendo à solução

P

vθ u

M = 1 N/m

L = 100 m

90

trivial (deslocamento igual a zero). A Fig. 4.32e exibe vibrações moderadas e, finalmente, na

Fig. 4.32f, verifica-se que a carga P = 1,00 Pcr provoca vibrações de grande amplitude,

ocasionando a perda de estabilidade da coluna, isto é, a resposta permanente não é mais

trivial.

0 1000 2000 3000 4000

Tempo [s]

-0.014

-0.012

-0.01

-0.008

-0.006

-0.004

-0.002

0

0.002

Des

loca

men

to (u

) [m

]

(a) P = 0,50 Pcr

0 1000 2000 3000 4000

Tempo [s]

-0.014

-0.012

-0.01

-0.008

-0.006

-0.004

-0.002

0

0.002

Des

loca

men

to (u

) [m

]

(b) P = 0,60 Pcr

0 1000 2000 3000 4000

Tempo [s]

-0.014

-0.012

-0.01

-0.008

-0.006

-0.004

-0.002

0

0.002

Des

loca

men

to (u

) [m

]

(c) P = 0,70 Pcr

0 1000 2000 3000 4000

Tempo [s]

-0.014

-0.012

-0.01

-0.008

-0.006

-0.004

-0.002

0

0.002

Des

loca

men

to (u

) [m

]

(d) P = 0,80 Pcr

0 1000 2000 3000 4000

Tempo [s]

-0.014

-0.012

-0.01

-0.008

-0.006

-0.004

-0.002

0

0.002

Des

loca

men

to (u

) [m

]

(e) P = 0,90 Pcr

0 1000 2000 3000 4000

Tempo [s]

-0.014

-0.012

-0.01

-0.008

-0.006

-0.004

-0.002

0

0.002

Des

loca

men

to (u

) [m

]

(f) P = 1,00 Pcr

Figura 4.32 – Resposta no tempo da coluna engastada para um carregamento súbito P.

A curva de Budianski é uma ferramenta muito útil na avaliação da instabilidade

dinâmica das estruturas. O critério de instabilidade dinâmica de Budianski diz que, resolvendo

as equações de movimento para alguns valores do parâmetro de carga, partindo de valores

pequenos e incrementando-os gradativamente, o sistema apresentará pequenas oscilações para

valores do parâmetro de carga também pequenos. A amplitude máxima da resposta aumentará

suavemente com o parâmetro de carga e a certo nível do parâmetro a amplitude da resposta

experimentará um grande salto. O valor em que este salto acontece é identificado como a

carga dinâmica crítica (SIMITSES e HODGES, 2006).

A Fig. 4.33 descreve bem esse fenômeno, relacionando o parâmetro de carga

crPP com a amplitude máxima da resposta no tempo da coluna engastada da Fig. 4.31. Com

o incremento do parâmetro de carga observa-se que a curva apresenta um grande salto para

9,0≅crPP , sendo esta a carga dinâmica crítica.

91

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

P / Pcr

0

-0.002

-0.004

-0.006

-0.008

-0.01

-0.012

-0.014

u máx

Figura 4.33 – Curva de Budianski para coluna engastada sob carga axial súbita P.

A Fig. 4.35 mostra a evolução da resposta no tempo dos deslocamentos

horizontais (u) do topo da coluna estaiada da Fig.4.8 (a se ver novamente na Fig. 4.34),

quando submetida a valores incrementais da carga súbita P. A carga crítica da coluna é de

Pcr = 126,01 kN, os cabos possuem inclinação igual a =α 60° e o tensionamento aplicado é

de 10=T kN.

Figura 4.34 – Coluna estaiada com dois cabos a 60°.

Verifica-se, nas Figuras 4.35a, 4.35b e 4.35c, que para o carregamento aplicado o

sistema volta ao estado fundamental. Isto é, após uma perturbação inicial, a amplitude da

resposta decresce rapidamente, convergindo para a solução trivial. Novamente, na Fig. 4.35d,

após uma perturbação inicial, a amplitude da resposta decresce rapidamente, porém,

apresentando vibrações de pequena amplitude no estado não fundamental. Nas Figuras 4.35e

e 4.35f, verifica-se que o carregamento aplicado ocasiona perda de estabilidade com

vibrações de grande amplitude.

L=100 m

α α

M = 1 N/m

θv

u

P

92

0 200 400 600 800 1000 1200

Tempo [s]

-0.008

-0.006

-0.004

-0.002

0

0.002

0.004D

eslo

cam

ento

(u) [

m]

(a) P = 0,92 Pcr

0 200 400 600 800 1000 1200

Tempo [s]

-0.008

-0.006

-0.004

-0.002

0

0.002

0.004

Des

loca

men

to (u

) [m

]

(b) P = 0,93 Pcr

0 200 400 600 800 1000 1200

Tempo [s]

-0.008

-0.006

-0.004

-0.002

0

0.002

0.004

Des

loca

men

to (u

) [m

]

(c) P = 0,94 Pcr

0 200 400 600 800 1000 1200

Tempo [s]

-0.008

-0.006

-0.004

-0.002

0

0.002

0.004

Des

loca

men

to (u

) [m

]

(d) P = 0,95 Pcr

0 200 400 600 800 1000 1200

Tempo [s]

-0.008

-0.006

-0.004

-0.002

0

0.002

0.004

Des

loca

men

to (u

) [m

]

(e) P = 0,96 Pcr

0 200 400 600 800 1000 1200

Tempo [s]

-0.008

-0.006

-0.004

-0.002

0

0.002

0.004

Des

loca

men

to (u

) [m

]

(f) P = 0,97 Pcr

Figura 4.35 – Resposta no tempo da coluna estaiada para um carregamento súbito P.

A Fig. 4.36 mostra a curva de Budianski desta coluna, onde observa-se que a

carga dinâmica crítica é igual à 95,0≅crPP .

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

P / Pcr

-0.003

-0.004

-0.005

-0.006

-0.007

-0.008

-0.009

u máx

Figura 4.36 – Curva de Budianski para coluna estaiada sob carga axial súbita P.

4.4.2 INSTABILIDADE DINÂMICA PARA CARGA HARMÔNICA

Na Fig. 4.38 ilustra-se a evolução da resposta no tempo dos deslocamentos

horizontais (u) do topo da coluna engastada da Fig. 4.37 quando sujeita a um carregamento

93

axial harmônico de )cos( tPP a Ω= , sendo aP uma parcela da carga crítica Pcr =16,3517 kN,

e Ω igual à freqüência natural da estrutura, que neste caso é 0,074332 rad/seg.

Figura 4.37 – Coluna engastada com carregamento harmônico )cos( tPP a Ω= .

0 5000 10000 15000

Tempo [s]

-0.004

-0.002

0

0.002

0.004

Des

loca

men

to (u

) [m

]

(a) aP = 0,01 Pcr

0 20000 40000 60000 80000 100000

Tempo [s]

-0.012

-0.008

-0.004

0

0.004

0.008

Des

loca

men

to (u

) [m

]

(b) aP = 0,02 Pcr

0 20000 40000 60000 80000 100000

Tempo [s]

-0.012

-0.008

-0.004

0

0.004

0.008D

eslo

cam

ento

(u) [

m]

(c) aP = 0,04 Pcr

0 20000 40000 60000 80000 100000

Tempo [s]

-0.2

-0.1

0

0.1

0.2

Des

loca

men

to (u

) [m

]

(d) aP = 0,06 Pcr

0 20000 40000 60000 80000 100000

Tempo [s]

-0.2

-0.1

0

0.1

0.2

Des

loca

men

to (u

) [m

]

(e) aP = 0,08 Pcr

0 20000 40000 60000 80000 100000

Tempo [s]

-0.2

-0.1

0

0.1

0.2

Des

loca

men

to (u

) [m

]

(f) aP = 0,10 Pcr

Figura 4.38 – Resposta no tempo da coluna engastada para um carregamento harmônico

)cos( tPP a Ω= .

P = Pa cos(Ωt)

vθ u

M = 1 N/m

L = 100 m

94

Observa-se que a amplitude da resposta no tempo do topo da coluna aumenta com

o incremento do parâmetro de carga, sendo mais significativo a partir da Fig. 4.38d. Observa-

se ainda que a amplitude da resposta exibe um crescimento exponencial ilimitado (Fig. 4.38f).

Na Fig. 4.39 é mostrada a curva de Budianski desta coluna, onde se verifica que a amplitude

máxima da resposta aumenta gradativamente com o incremento do parâmetro de carga.

0 0.02 0.04 0.06 0.08 0.1

Pa / Pcr

0

0.02

0.04

0.06

0.08

0.1

0.12

u máx

Figura 4.39 – Curva de Budianski para coluna engastada sob carga axial harmônica

)cos( tPP a Ω= .

Considera-se também o caso onde a coluna engastada da Fig. 4.40 está sujeita ao

carregamento axial harmônico )2cos( tPP a Ω= , sendo novamente, aP uma parcela da carga

crítica Pcr =16,3517 kN e Ω a freqüência natural da estrutura, isto é, 0,074332 rad/seg. Na

Fig. 4.41 mostra-se a evolução da resposta no tempo dos deslocamentos horizontais (u) do

topo da coluna.

Figura 4.40 – Coluna engastada com carregamento harmônico )2cos( tPP a Ω= .

P = Pa cos(2Ωt)

vθ u

M = 1 N/m

L = 100 m

95

0 2000 4000 6000 8000 10000

Tempo [s]

-6

-4

-2

0

2

4

6D

eslo

cam

ento

(u) [

m]

(a) aP = 0,05 Pcr

0 2000 4000 6000 8000 10000

Tempo [s]

-6

-4

-2

0

2

4

6

Des

loca

men

to (u

) [m

]

(b) aP = 0,06 Pcr

0 2000 4000 6000 8000 10000

Tempo [s]

-6

-4

-2

0

2

4

6

Des

loca

men

to (u

) [m

]

(c) aP = 0,07 Pcr

0 2000 4000 6000 8000 10000

Tempo [s]

-6

-4

-2

0

2

4

6

Des

loca

men

to (u

) [m

]

(d) aP = 0,08 Pcr

0 2000 4000 6000 8000 10000

Tempo [s]

-6

-4

-2

0

2

4

6

Des

loca

men

to (u

) [m

]

(e) aP = 0,09 Pcr

0 2000 4000 6000 8000 10000

Tempo [s]

-6

-4

-2

0

2

4

6

Des

loca

men

to (u

) [m

]

(f) aP = 0,10 Pcr

Figura 4.41 – Resposta no tempo da coluna engastada para um carregamento harmônico

)2cos( tPP a Ω= .

Observa-se que, inicialmente, a resposta apresenta amplitude de vibração

praticamente nula, mas, a partir de certo tempo, nota-se uma mudança nas características da

resposta, passando a apresentar um crescimento exponencial da amplitude. Com o aumento do

parâmetro de carga, o tempo decorrido para resposta apresentar essa mudança, diminui. Por

exemplo, na Fig. 4.41c o tempo é de 6500 segundos e na Fig. 4.41d é de 6000 segundos. Nas

Figuras 4.41e e 4.41f, observa-se que a resposta apresenta deslocamentos de grande

amplitude.

Verifica-se na curva de Budianski desta coluna (Fig. 4.42) que a amplitude

máxima da resposta aumenta gradualmente a partir de 04,0≅crPP . Com o incremento do

parâmetro de carga, identifica-se que a carga dinâmica crítica para 05,0≅crPP . Após o salto,

a curva apresenta um aumento discreto na amplitude da resposta.

A Fig. 4.44 apresenta seis respostas no tempo do topo da coluna estaiada mostrada

na Fig. 4.43, que está sob ação do carregamento axial harmônico )cos( tPP a Ω= , com aP

sendo a parcela da carga crítica Pcr =126,01 kN, e Ω a freqüência natural da estrutura, igual a

96

0,352215 rad/seg. Recordando que considera-se tensionamento constante 10=T kN e

inclinação =α 60° em ambos os cabos. Na Fig. 4.46 tem-se a resposta no tempo para o

mesmo sistema estrutural, considerando, no entanto, )2cos( tPP a Ω= .

0 0.02 0.04 0.06 0.08 0.1

Pa / Pcr

0

0.5

1

1.5

2

2.5

3

3.5

4

4.5

5

u máx

Figura 4.42 – Curva de Budianski para coluna engastada sob carga axial harmônica

)2cos( tPP a Ω= .

Figura 4.43 – Coluna estaiada com dois cabos a 60° e carga axial harmônica )cos( tPP a Ω= .

Nas Figuras 4.44a-d observa-se que, após uma perturbação inicial, a amplitude da

resposta decresce rapidamente, convergindo para a solução trivial. Na Fig. 4.44e observa-se

que a amplitude da resposta apresenta, após decrescer rapidamente, vibrações de pequena

amplitude. Na Fig. 4.44f nota-se uma mudança de comportamento, exibindo um crescimento

exponencial ilimitado.

L=100 m

α α

M = 1 N/m

θv

u

P = Pa cos(Ωt)

97

0 200 400 600 800 1000 1200

Tempo [s]

-0.004

-0.003

-0.002

-0.001

0

0.001

0.002

0.003

0.004D

eslo

cam

ento

(u) [

m]

(a) aP = 1,40 Pcr

0 200 400 600 800 1000 1200

Tempo [s]

-0.004

-0.003

-0.002

-0.001

0

0.001

0.002

0.003

0.004

Des

loca

men

to (u

) [m

]

(b) aP = 1,60 Pcr

0 200 400 600 800 1000 1200

Tempo [s]

-0.004

-0.003

-0.002

-0.001

0

0.001

0.002

0.003

0.004

Des

loca

men

to (u

) [m

]

(c) aP = 1,80 Pcr

0 200 400 600 800 1000 1200

Tempo [s]

-0.004

-0.003

-0.002

-0.001

0

0.001

0.002

0.003

0.004

Des

loca

men

to (u

) [m

]

(d) aP = 2,00 Pcr

0 200 400 600 800 1000 1200

Tempo [s]

-0.004

-0.003

-0.002

-0.001

0

0.001

0.002

0.003

0.004

Des

loca

men

to (u

) [m

]

(e) aP = 2,20 Pcr

0 200 400 600 800 1000 1200

Tempo [s]

-0.004

-0.003

-0.002

-0.001

0

0.001

0.002

0.003

0.004

Des

loca

men

to (u

) [m

]

(f) aP = 2,40 Pcr

Figura 4.44 – Resposta no tempo da coluna estaiada para um carregamento harmônico

)cos( tPP a Ω= .

Na curva de Budianski desta coluna (Fig. 4.45) verifica-se que a amplitude

máxima da resposta aumenta rapidamente a partir de 70,1≅crPP , seguido de um grande

salto, onde se identifica a carga dinâmica crítica em 30,2≅crPP .

0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 2.2 2.4

Pa / Pcr

0

0.0005

0.001

0.0015

0.002

0.0025

0.003

0.0035

u máx


)cos( tPP a Ω= .

98

No conjunto de respostas no tempo da Fig. 4.46 notam-se mudanças bruscas de

suas características à medida que se incrementa o parâmetro de carga. Nas Figuras 4.46a e

4.46b observa-se que, após uma perturbação inicial, a amplitude da resposta decresce

rapidamente, convergindo para solução trivial. Nas Figuras 4.46c e 4.46d a resposta apresenta,

após decrescer rapidamente, alternância entre ciclos de aumento e de redução na amplitude

das vibrações. Finalmente, as Figuras 4.46e e 4.46f exibem, a partir de certo tempo, um

crescimento abrupto da amplitude, convergindo para configuração de equilíbrio trivial.

0 300 600 900 1200 1500 1800

Tempo [s]

-0.008

-0.004

0

0.004

0.008

Des

loca

men

to (u

) [m

]

(a) aP = 0,50 Pcr

0 300 600 900 1200 1500 1800

Tempo [s]

-0.008

-0.004

0

0.004

0.008D

eslo

cam

ento

(u) [

m]

(b) aP = 0,60 Pcr

0 300 600 900 1200 1500 1800

Tempo [s]

-0.008

-0.004

0

0.004

0.008

Des

loca

men

to (u

) [m

]

(c) aP = 0,70 Pcr

0 300 600 900 1200 1500 1800

Tempo [s]

-0.008

-0.004

0

0.004

0.008

Des

loca

men

to (u

) [m

]

(d) aP = 0,71 Pcr

0 300 600 900 1200 1500 1800

Tempo [s]

-8

-4

0

4

8

Des

loca

men

to (u

) [m

]

(e) aP = 0,72 Pcr

0 300 600 900 1200 1500 1800

Tempo [s]

-8

-4

0

4

8

Des

loca

men

to (u

) [m

]

(f) aP = 0,73 Pcr

Figura 4.46 – Resposta no tempo da coluna estaiada para um carregamento harmônico

)2cos( tPP a Ω= .

Na curva de Budianski da Fig. 4.47 verifica-se que, após um longo trecho sem

apresentar variações nas máximas amplitudes da resposta, ocorre um salto brusco na curva

para 70,0≅crPP , seguido de um aumento gradual com o incremento do parâmetro de carga.

99

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

Pa / Pcr

0

0.25

0.5

0.75

1

1.25

1.5

1.75

2

2.25

2.5

2.75

3

3.25

u máx


)2cos( tPP a Ω= .

4.5 ANÁLISE ESTÁTICA LINEAR VIA MODELO DISCRETO


A Fig. 4.49 mostra os caminhos pós-crítico obtidos a partir da Eq. 2.80 para

diferentes valores de α1 da coluna estaidada apresentada na Fig. 4.48. A análise estática da

estabilidade de torres estaiadas, modeladas com molas lineares, foi aplicada anteriormente por

Pasquetti (2003) e os resultados apresentados estão em concordância com os obtidos em seu

trabalho. Pode-se observar que, dependendo do valor do ângulo α1, o caminho pós-crítico é

instável (Figuras 4.49a e 4.49d) ou terá ganho de rigidez inicial e posterior perda de

estabilidade (Figuras 4.49b e 4.49c).

Figura 4.48 – Coluna estaiada com duas molas lineares.

100

(a) º31,1,0,0 110 ==== αξξξ p

(b) º32,1,0,0 110 ==== αξξξ p

(c) º58,1,0,0 110 ==== αξξξ p

(d) º59,1,0,0 110 ==== αξξξ p

Figura 4.49 – Caminhos pós-críticos.

4.6 ANÁLISE DINÂMICA LINEAR VIA MODELO DISCRETO

O objetivo dessa seção é estudar o comportamento dinâmico não-linear da torre

estaiada quando submetida a um carregamento periódico. Para este fim, utiliza-se como

referência a coluna estaiada da Fig. 4.50, colocando no lugar dos cabos molas com

comportamento linear. Em ambas as molas considera-se um tensionamento aplicado igual a

10 kN e uma inclinação de º321 =α .

θ [rad.]

λP1

θ [rad.]

λP1

θ [rad.]

λP1

θ [rad.]

λP1

101

Figura 4.50 – Coluna estaiada com dois cabos a 32° e carga axial periódica )cos( tPP a Ω= .

O carregamento axial periódico é da forma )cos( tPP a Ω= , sendo Ω igual à

freqüência natural de vibração da estrutura, calculada indiretamente, através da solução

numérica das equações de movimento do sistema autônomo, isto é, na ausência do

amortecimento e das forças externas, em função da inclinação do cabo 1α . Para isso

consideram-se pequenas perturbações iniciais, integram-se as equações no tempo e mede-se o

período T. Logo, Tπω 20 = que, neste caso é 0,25234. Na Tabela 4.7 estão indicados os

parâmetros físicos e geométricos adotados na análise.

Tabela 4.7 – Parâmetros físicos e geométricos da coluna e das molas.

Propriedades Coluna Molas

Área (A) 1,914 x 10 -02 m² -

Densidade ( ρ ) 78,50 kN/m³ -

Constante de Rigidez (K) - 2,30 kN/m

4.6.1 INSTABILIDADE PARAMÉTRICA

Na Fig. 4.51 é mostrada a fronteira de instabilidade paramétrica no espaço de

controle (freqüência de excitação Ω - amplitude da excitação aP ) da coluna estaiada com

carga brusca. Para se obter este limite, integrou-se, para cada valor da freqüência de

102

excitação, as equações de movimento, aumentando-se a amplitude do carregamento e

verificando-se a cada incremento de carga as características da resposta permanente. O limite

inferior da fronteira de instabilidade paramétrica representa os parâmetros nos quais pequenas

perturbações levam o sistema de volta à solução trivial ou fundamental. O limite superior

corresponde ao escape do vale potencial de pré-flambagem (instabilidade dinâmica), isto é, à

solução não trivial. O vale mais profundo está associado à região principal de instabilidade,

com freqüência de excitação igual a Oω2 , enquanto o vale imediatamente a esquerda é a

região secundária de instabilidade, com freqüência de excitação igual a Oω , sendo Oω a

freqüência natural de vibração da torre. Os vales menores mais à esquerda estão relacionados

com as ressonâncias sub-harmônicas.

0.00 0.50 1.00 1.50 2.00 2.50 3.00

Ω

0.00

0.50

1.00

1.50

2.00

2.50

3.00

3.50

4.00

4.50

5.00

Pa

Ω = 2ωοΩ = ωο

2

5

6

7

4

1

3

Figura 4.51 – Fronteira de Instabilidade Paramétrica.

103

4.6.2 MECANISMOS DE ESCAPE

Com o objetivo de estudar a estabilidade das respostas do sistema no espaço dos

parâmetros da carga e as bifurcações associadas às fronteiras de estabilidade, foi realizado um

estudo das bifurcações que ocorrem nas duas regiões mais importantes de instabilidade

paramétrica. Para isto foram obtidos diagramas de bifurcação através do algoritmo de força

bruta, usando-se como parâmetro de controle a amplitude da força.

Na Fig. 4.52 é mostrado o diagrama de bifurcação da região secundária de

instabilidade paramétrica mostrada na Fig. 4.51 para 75,0=Ω (seção 1). Na Fig. 4.52a

observa-se que, incrementando a amplitude da força, a solução trivial em certo ponto perde a

estabilidade e salta, através de uma bifurcação do tipo sub-crítica, para uma solução não

trivial com grande amplitude de vibração e período 2T, tanto para frente quanto para traz,

como pode ser visto na Fig. 4.52b. O trecho construído com pontos na cor preta foi obtido

variando a amplitude de excitação da força de forma crescente (1,0 a 1,9), enquanto que, o

trecho construído com pontos na cor azul foi obtido variando a mesma amplitude de forma

decrescente (1,9 a 1,0).

1 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9

Pa

-100

0

100

200

300

400

500

600

700

θ

(a)

1.6 1.65 1.7 1.75 1.8 1.85 1.9

Pa

624.8

624.9

625

625.1

625.2

625.3

625.4

625.5

θ

(b)

Figura 4.52 – Mecanismo de escape para 75,0=Ω (seção 1): (a) Diagrama de bifurcação; (b)

Detalhe do diagrama de bifurcação.

O diagrama de bifurcação mostrado na Fig. 4.53a é típico do ramo ascendente da

fronteiras de estabilidade paramétrica secundária mostrada na Fig. 4.51 e refere-se à região

com 10,1=Ω (seção 2).

104


Resposta no tempo A; (c) Resposta no tempo B; (d) Plano fase e mapeamento de Poincaré B;

(e) Resposta no tempo C; (f) Plano fase e mapeamento de Poincaré C.

Analisando o diagrama de bifurcação da Fig. 4.53a observa-se que, incrementando

a amplitude da força, a solução trivial perde a estabilidade através de uma bifurcação do tipo

super-crítica e surgem duas soluções não triviais estáveis de período 1T (duplicação de

solução). Continuando a incrementar a amplitude da força, as vibrações crescem até ocorrer a

duplicação do período.

Para auxiliar o entendimento do mecanismo de escape nessa região, são mostradas

algumas respostas no tempo, os planos fase e as seções de Poincaré de pontos relativos ao

diagrama de bifurcação da Fig. 4.53a. A Fig. 4.53b confirma que a torre, após a perturbação

inicial, apresenta decaimento exponencial das amplitudes até retornar à solução trivial. Nas

Figuras 4.53c e 4.53d confirmam-se as soluções não triviais de período 1T e nas Figuras 4.53e

e 4.53f confirmam-se as soluções não triviais de período 2T.

O mecanismo de escape para 20,1=Ω (seção 3) da fronteiras de estabilidade

paramétrica da Fig. 4.51 é mostrado na Fig. 4.54. Verifica-se que essa, região possui o mesmo

comportamento apresentado para Fig.4.53, porém, sem apresentar duplicação de período.

(b)

(a)

2 2.1 2.2 2.3 2.4 2.5 2.6 2.7 2.8 2.9 3

Pa

-0.08

-0.06

-0.04

-0.02

0

0.02

0.04

0.06

0.08

θA

B

C

B'

C'

0 200 400 600 800 1000

t

-0.06

-0.04

-0.02

0

0.02

0.04

0.06

θ

-0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8

θ

-1.4-1.2

-1-0.8-0.6-0.4-0.2

00.20.40.60.8

11.21.4

θ

C

C'

-0.6 -0.4 -0.2 0 0.2 0.4 0.6

θ

-0.8

-0.4

0

0.4

0.8

θ

B'

B

0 200 400 600 800 1000

t

-0.8

-0.4

0

0.4

0.8

θ

0 200 400 600 800 1000

t

-0.6

-0.4

-0.2

0

0.2

0.4

0.6

θ

(d) (f)

(e) (c)

105


Resposta no tempo A; (c) Resposta no tempo B; (d) Plano fase e mapeamento de Poincaré B.

Na Fig. 4.55 é mostrado o diagrama de bifurcação do ramo descendente da região

principal de instabilidade paramétrica da Fig. 4.51 para 80,1=Ω (seção 4). Nessa região, ao

se incrementar a amplitude da força, observa-se uma bifurcação do tipo sub-crítica na qual, a

partir de certo ponto, a torre salta para uma região com grande amplitude de vibrações e

período 2T.

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8

Pa

-2

0

2

4

6

8

10

12

θ

Figura 4.55 – Mecanismo de escape para 80,1=Ω (seção 4): diagrama de bifurcação.

O diagrama de bifurcação mostrado na Fig. 4.56a é típico do ramo ascendente da

fronteira principal de estabilidade paramétrica mostrada na Fig. 4.51 e refere-se à região com

10,2=Ω (seção 5).

(a)

3 3.25 3.5 3.75 4 4.25

Pa

-0.006

-0.005

-0.004

-0.003

-0.002

-0.001

0

0.001

0.002

0.003

0.004

0.005

0.006

θA

B

B'

(b) (d)

(c)

-0.4 -0.2 0 0.2 0.4

θ

-0.5

-0.4

-0.3

-0.2

-0.1

0

0.1

0.2

0.3

0.4

0.5

θ

B

B'

0 200 400 600 800 1000

t

-0.4

-0.2

0

0.2

0.4

θ

0 200 400 600 800 1000

t

-0.04

-0.02

0

0.02

0.04

θ

106


Resposta no tempo do ponto A; (c) Resposta no tempo do ponto B; (d) Plano fase e

mapeamento de Poincaré do ponto B.

Nessa região, quando se incrementa a amplitude da força, a solução trivial perde a

estabilidade, a partir de certo ponto, pó meio de uma bifurcação do tipo super-crítica e surge

uma solução não trivial estável de período 2T cuja amplitude, inicialmente, cresce

rapidamente a partir do valor crítico e, posteriormente, decresce continuamente com

suavidade.

Para se visualizar melhor o escape nessa região, foram traçadas algumas respostas

no tempo, os planos fase e as seções de Poincaré para pontos relativos ao diagrama de

bifurcação da Fig. 4.56a. Na Fig. 4.56b fica evidente que, após uma perturbação inicial, a

amplitude da resposta decresce rapidamente, convergindo para a solução trivial. Na Fig.

4.56c, confirma-se a mudança de comportamento da estrutura, passando a apresentar

perturbações durante toda a resposta e na Fig. 4.56d o mapa de Poincaré novamente sugere a

presença de vibrações periódicas do tipo 2T.

Nas Figuras 4.57 e 4.58 são mostrados os mecanismos de escape para 30,2=Ω

(seção 6) e 60,2=Ω (seção 7) da fronteira principal de estabilidade paramétrica mostrada na

Fig. 4.51.

(a)

0 0.25 0.5 0.75 1 1.25 1.5

Pa

-0.05

-0.04

-0.03

-0.02

-0.01

0

0.01

0.02

0.03

0.04

0.05

θA

B

(b) (d)

(c)

-1.2 -0.8 -0.4 0 0.4 0.8 1.2

θ

-1.2

-0.8

-0.4

0

0.4

0.8

1.2

θ

0 200 400 600 800 1000

t

-2

-1

0

1

2

θ

0 200 400 600 800 1000

t

-0.015

-0.01

-0.005

0

0.005

0.01

0.015

θ

107







(a)

1 1.25 1.5 1.75 2 2.25 2.5

Pa

-0.015

-0.01

-0.005

0

0.005

0.01

0.015

θA

B

(b) (d)

(c)

-0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8

θ

-0.8

-0.6

-0.4

-0.2

0

0.2

0.4

0.6

0.8

θ

0 200 400 600 800 1000

t

-1.2

-0.8

-0.4

0

0.4

0.8

1.2

θ

0 200 400 600 800 1000

t

-0.04

-0.02

0

0.02

0.04

θ

(a)

0 0.25 0.5 0.75 1 1.25 1.5

Pa

-0.025

-0.02

-0.015

-0.01

-0.005

0

0.005

0.01

0.015

0.02

0.025

θA

B

(b) (d)

(c)

-1.2 -0.8 -0.4 0 0.4 0.8 1.2

θ

-1.2

-0.8

-0.4

0

0.4

0.8

1.2

θ

0 200 400 600 800 1000

t

-1.5

-1

-0.5

0

0.5

1

1.5

θ

0 200 400 600 800 1000

t

-0.02

-0.01

0

0.01

0.02

θ

108

Nessas regiões observa-se que o comportamento da estrutura é semelhante ao

visto no mecanismo de escape da Fig.4.56. Ainda assim, nos diagramas de bifurcação das

Figuras 4.57a e 4.58a verifica-se que a amplitude da solução não trivial de período 2T,

diferente do diagrama de bifurcação da Fig. 4.56a, apresenta crescimento contínuo ao longo

do caminho pós-crítico.

É importante destacar que um sistema onde o parâmetro de controle varia de

forma suave antes que o escape aconteça a solução de período 2T pode se tornar instável ou

seguir por uma cascata de duplicação de período.

Dos resultados obtidos pode-se afirmar que o ramo descendente da região

principal de instabilidade está associada à bifurcação sub-crítica e o ramo ascendente está

associada à bifurcação super-crítica. Observa-se ainda que, nos ramos ascendentes das regiões

principal e secundária de instabilidade paramétrica, à medida que se aumenta a freqüência de

excitação Ω a amplitude máxima de vibração diminui.

4.6.3 EVOLUÇÃO DA ESTABILIDADE GLOBAL

Os resultados até aqui apresentados não são suficientes para se projetar com

segurança uma estrutura sob cargas dinâmicas. Para isto, deve-se ter uma visão global do

comportamento da estrutura na presença de perturbações que podem ocorrer durante a

construção e vida útil da estrutura. O estudo da evolução das bacias de atração no espaço dos

parâmetros de controle tem-se mostrado uma ferramenta útil no estudo da estabilidade global

de sistemas dinâmicos (SANTEE, 1999; SOLIMAN, 1995). Esta tarefa relativamente fácil

quando se lida com sistemas com um grau de liberdade, torna-se bastante complexa para o

engenheiro e computacionalmente cara quando se trabalha com sistemas com vários graus de

liberdade. Entretanto, mesmo nestes casos, projeções criteriosas da bacia de atração podem

fornecer informações importantes sobre o grau de segurança da estrutura.

A bacia de atração de pontos fixos ou soluções periódicas está baseada no limite

quando ∞→t . Em termos numéricos, alcançar esse tempo é impossível. Por essa razão, usas-

se a chamada bacia de atração transiente (SOLIMAN, 1995), que é a bacia de atração onde a

integração numérica do sistema é truncada num múltiplo do período Ω= π2T da força

excitadora. À medida que o tempo de integração aumenta, a bacia de atração transiente

converge para a bacia de atração real.

109

Existem duas maneiras distintas de evolução da bacia de atração em função da

intensidade da força excitadora (SANTEE, 1999), a saber:

a) Redução gradativa da bacia de atração até seu completo desaparecimento

quando se atinge a fronteira de estabilidade.

b) Perda de integridade da bacia de atração, cuja fronteira torna-se subitamente

fractal, sendo este fenômeno seguido por um rápido processo de erosão e

estratificação da bacia.

Porém, em ambos os casos, a área da bacia de atração da solução periódica estável

tende a zero.

Na Fig. 4.59 são mostrados dois diagramas de bifurcação estudados

anteriormente. A seguir, apresenta-se a evolução das bacias de atração com o aumento da

amplitude de excitação para estes casos. As Figuras 4.60 e 4.61 correspondem ao plano •

−θθ .

As bacias foram geradas a partir de uma malha igualmente espaçada de 200x200, ou seja,

40000 células.

2 2.1 2.2 2.3 2.4 2.5 2.6 2.7 2.8 2.9 3

Pa

-0.08

-0.06

-0.04

-0.02

0

0.02

0.04

0.06

0.08

θA

B

(a) Ω =1,10

0 0.25 0.5 0.75 1 1.25 1.5

Pa

-0.05

-0.04

-0.03

-0.02

-0.01

0

0.01

0.02

0.03

0.04

0.05

θ A

(b) Ω =2,10

Figura 4.59 – Diagramas de bifurcação para estudo de erosão da bacia de atração.

Na Fig. 4.60 são vistas as bacias de atração correspondentes ao diagrama de

bifurcação da Fig. 4.59a. Nestas bacias observam-se quatro regiões distintas associadas a três

tipos distintos de comportamento, a saber: a) região azul – solução fundamental; b) região

verde e vermelha – solução pós-crítica estável e c) região branca – escape.

110

θ

θ

0.38

0.34

0.300.26

0.22

0.18

-0.02

0.140.10

0.06

0.02

-0.302 -0.004 0.024 0.052 0.080-0.060

(a) aP =2,35

θ

θ

0.38

0.34

0.300.26

0.22

0.18

-0.02

0.140.10

0.06

0.02

-0.302 -0.004 0.024 0.052 0.080-0.060

(b) aP =2,40

θ

θ

0.38

0.34

0.300.26

0.22

0.18

-0.02

0.140.10

0.06

0.02

-0.302 -0.004 0.024 0.052 0.080-0.060

(c) aP =2,55

θ

θ

0.38

0.34

0.300.26

0.22

0.18

-0.02

0.140.10

0.06

0.02

-0.302 -0.004 0.024 0.052 0.080-0.060

(d) aP =2,95

Figura 4.60 – Evolução da bacia de atração para o diagrama de bifurcação da Fig.4.59a.

As Figuras 4.60a mostra a resposta para um valor de força inferior ao ponto A na

Fig. 4.59a. Observa-se que todo o espaço das condições iniciais está associado à solução

trivial, estando envolvida por uma bacia contínua. Outras bacias são mostradas nas Figuras

4.60b-d. Após o ponto A, a bacia associadas às soluções triviais (azul) desaparece e aparecem

dois atratores não-triviais.

Na Fig. 4.61 são vistas as bacias de atração correspondentes ao diagrama de

bifurcação da Fig. 4.59b. Nestas projeções observam-se quatro regiões distintas associadas a

três tipos distintos de comportamento, a saber: a) região azul – solução fundamental; b) região

verde e vermelha – solução pós-crítica estável e c) região branca – escape.

θ

θ

-0.4 -0.3 -0.2 -0.1 0.0 0.1 0.2 0.3 0.4 0.5-0.5

1.20

0.96

0.720.48

0.24

0.00

-0.24

-0.48

-0.72

-0.96

-1.20

(a) aP =0,200

θ

θ

1.20

0.96

0.720.48

0.24

0.00

-0.24

-0.48

-0.72

-0.96

-1.20-0.4 -0.3 -0.2 -0.1 0.0 0.1 0.2 0.3 0.4 0.5-0.5

(b) aP =0,275

111

θ

θ

1.20

0.96

0.720.48

0.24

0.00

-0.24

-0.48

-0.72

-0.96

-1.20-0.4 -0.3 -0.2 -0.1 0.0 0.1 0.2 0.3 0.4 0.5-0.5

(c) aP =0,625

θ

θ

1.20

0.96

0.720.48

0.24

0.00

-0.24

-0.48

-0.72

-0.96

-1.20-0.4 -0.3 -0.2 -0.1 0.0 0.1 0.2 0.3 0.4 0.5-0.5

A

B

(d) aP =0,750

Figura 4.61 – Evolução da bacia de atração para o diagrama de bifurcação da Fig.4.59b.

As Figuras 4.61a mostra a bacia para valores de força inferiores ao ponto A da

Fig. 4.59b. Observa-se que todo o espaço das condições iniciais na Fig. 4.61a está associado à

solução trivial, sendo esta solução envolvida por uma bacia contínua. A partir do ponto A, a

região associada à solução trivial desaparece, surgindo as regiões associadas às soluções não

triviais (Figuras 4.61b-d), sendo que, à medida que crescem as excitações, a região decresce

paulatinamente conduzindo-se ao escape.

Para ilustrar a periodicidade das soluções não triviais, mostra-se nas Figuras 4.62

a resposta no tempo, o plano fase e o mapeamento de Poincaré de dois pontos na bacia de

atração da Fig. 4.61d.

0 200 400 600 800 1000

t

-1.2

-0.8

-0.4

0

0.4

0.8

1.2

θ

0 200 400 600 800 1000

t

-1.2

-0.8

-0.4

0

0.4

0.8

1.2

θ

-1.2 -0.8 -0.4 0 0.4 0.8 1.2

θ

-1.2

-0.8

-0.4

0

0.4

0.8

1.2

θ

(a)

-1.2 -0.8 -0.4 0 0.4 0.8 1.2

θ

-1.2

-0.8

-0.4

0

0.4

0.8

1.2

θ

(b)

Figura 4.62 – Resposta no tempo, plano fase e mapeamento de Poincaré do ponto A (a);

Resposta no tempo, plano fase e mapeamento de Poincaré do ponto B (b).

112

5 CONCLUSÕES E SUGESTÕES

Embora exista uma grande quantidade de resultados teóricos e experimentais

sobre o comportamento não-linear geométrico das estruturas estaiadas e sua perda de

estabilidade sob cargas axiais e súbitas, ainda se faz necessário conhecer mais seu

comportamento não-linear e sua perda de estabilidade sob cargas periódicas.

Neste trabalho foi apresentada uma formulação do método dos elementos finitos e

um modelo discreto para o estudo do comportamento não-linear geométrico de estruturas

estaiadas sob cargas estáticas e dinâmicas, dando-se ênfase aos efeitos das imperfeições

geométricas iniciais, à instabilidade dinâmica, à instabilidade paramétrica, aos mecanismos de

escape e à evolução das bacias de atração.

No método dos elementos finitos foi adotado o referencial Lagrangiano atualizado

onde os elementos de pórticos foram representados por elementos de viga/coluna modelados

utilizando a teoria de flexão de Euler-Bernoulli e os cabos foram representados por elementos

de cabo/treliça. Dentre os métodos numéricos, o método de Newton-Raphson foi utilizado

para obtenção das diversas configurações de equilíbrio. Sua associação com o método do

comprimento de arco permitiu a obtenção das trajetórias de equilíbrio avançadas. No

algoritmo de Newmark, usado para a integração das equações de movimento ao longo do

tempo, o acúmulo de resíduos a cada incremento, devido às simplificações e linearizações

inerentes à formulação do elemento finito, foi contornado com a adoção de um passo de

tempo extremamente reduzido. Assim, apesar de o método mostrar-se adequado para as torres

estaiadas analisadas neste trabalho, o processamento tornou-se dispendioso nas análises

envolvendo intervalos muito grandes de tempo. Por este método foi possível, ainda, avaliar a

instabilidade dinâmica das torres sob cargas harmônicas axiais pelo critério de Budianski.

113

No modelo discreto, foi adotado um esquema simplificado de torre estaiada,

modelando os cabos com molas com comportamento linear. Foram utilizadas ferramentas

numéricas para obter a fronteira de instabilidade no espaço dos parâmetros de controle e as

diversas bifurcações associadas às fronteiras de estabilidade. Através dos diagramas de

bifurcação e das respostas no tempo, foram estudados os mecanismos associados à perda de

estabilidade da torre estaiada. Pode-se verificar que, em cada região de instabilidade, o ramo

descendente está associado às bifurcações sub-críticas e o ramo ascendente às bifurcações

super-críticas. Além disso, verificou-se que, quando ocorrem bifurcações sub-críticas, a torre

sofre um salto com grandes deslocamentos, os quais, em muitos casos, podem causar danos e

até mesmo a ruína da estrutura.

Em ambos os métodos, o estudo das torres estaiadas evidenciou a forte influência

dos cabos no comportamento não-linear, havendo um forte acréscimo do valor da carga crítica

com sua inclusão. Em termos de características dinâmicas, os resultados obtidos mostraram

que o tensionamento dos cabos influencia as características dinâmicas da estrutura.

Quanto aos efeitos das imperfeições geométricas iniciais no comportamento

estático e dinâmico da torre, verificou-se que a torre, para algumas imperfeições, pode

apresentar freqüências naturais e/ou cargas críticas muito próximas ou coincidentes quando

comparados com a torre perfeita. Entretanto, dependendo da natureza da imperfeição, até

mesmo as de pequena magnitude podem provocar uma grande queda na capacidade de carga,

bem como variações na freqüência natural e nos modos de vibração da estrutura.

Comparando as duas formulações implementadas, constatou-se que a formulação

do método dos elementos finitos foi, computacionalmente, muito mais dispendiosa, assim, no

modelo discreto foi possível avaliar efeitos que no método dos elementos finitos dificilmente

seriam calculados.

Os resultados desse trabalho e a expansão das ferramentas computacionais

implementadas, podem ser usados para a análise de torres estaiadas reais sob outros

carregamentos dinâmicos como, por exemplo, na modelagem de problemas em três dimensões

sob a ação do vento, com amplitude e direção aleatórias.

Apesar de eficientes, faz-se necessário tornar os algoritmos mais amigáveis, no

que se refere à entrada de dados e, principalmente, à saída de resultados. Para isso, uma

interface gráfica torna-se desejável.

Para trabalhos futuros, também seria interessante estudar a influência do colapso

de um dado cabo nas vibrações da estrutura, aplicar técnicas de controle de vibração como,

114

por exemplo, amortecedores de massa sintonizado ou hibrido, otimizar o sistema

computacional e, no modelo discreto, utilizar molas não-lineares.

115

6 REFERÊNCIAS BIBLIOGRÁFICAS

ALVES, R. V. Instabilidade não-linear elástica de estruturas reticuladas espaciais. 1995. Tese (Doutorado) – COPPE, Universidade Federal do Rio de Janeiro, Rio de Janeiro, 1995.

AMARI, G. Ghodrati. Seismic sensitivity indicators for tall guyed telecommunication towers. Computers & Structures, nº 80, pp. 349-364, 2002.

ANDREU, A. et al. A new deformable catenary element for the analysis of cable net structures. Computers & Structures, n° 84, pp. 1882-1890, 2006.

BATHE, K. J. Finite elements procedures. Prentice Hall, 1996. 1037 f. BAZANT, Z. P. e CEDOLIN, L. Stabily of structures: elastic, inelastic, fracture, and

damage theories. Oxford University Press, 1991. 984 f. CAMPOS FILHO, Erlande da Costa. Análise do comportamento não-linear de estruturas

estaiadas planas. 2004. 109 f. Dissertação (Mestrado em Engenharia Civil) – Curso de Mestrado em Engenharia Civil, Universidade Federal de Goiás, Goiânia, 2004.

CHAN, S. L. et al. Stability analysis and parametric study of pre-stressed stayed columns. Engineering Structures, n° 24, pp. 115-124, 2002.

CHENG, J. et al. Advanced aerostatic stability analysis of cable-stayed bridges using finite-element method. Computers & Structures, nº 80, pp. 1145-1158, 2002.

CRISFIELD, M. A. Non-linear finite element analysis of solids and structures. John Willey & Sons Ltd, 1991. 345 f.

CROLL, J.G.A. e WALKER, A.C. Elementos de estabilidad estructural. Barcelona, Editora Reverté, 1975.

DEL PRADO, Z. J. G. N. Acoplamento e iteração modal na instabilidade de cascas cilíndricas. 2001. 222 f. Tese (Doutorado) – Pontifícia Universidade Católica do Rio de Janeiro, Rio de Janeiro, 2001.

FREIRE, A. M. S. et al. Geometrical nonlinearities on the static analysis of highly flexible steel cable-stayed bridges. Computers & Structures, nº 84, pp. 2128-2140, 2006.

GALVÃO, A. S. Formulações não-lineares de elementos finitos para análise de sistemas estruturais metálicos reticulados planos. 2000. 168 f. Dissertação (Mestrado) – Universidade Federal de Ouro Preto, Ouro Preto, 2000.

116

HARIKRISHNA, P. et al. Full scale measurements of the structural response of a 50 m guyed mast under wind loading. Engineering Structures, nº 25, pp. 859-867, 2003.

KAHLA, Nabil Ben. Nonlinear dynamic response of a guyed tower to a sudden guy rupture. Engineering Structures, v. 19, nº 11, pp. 879-890, 1997.

KAHLA, Nabil Ben. Response of a guyed tower to a guy rupture under no wind pressure. Engineering Structures, nº 22, pp. 699-706, 2000.

KAROUMI, Raid. Some modeling aspects in the nonlinear finite element analysis of cable supported bridges. Computers & Structures, nº 71, pp. 397-412, 1999.

LAW, S. S. et al. Time-varying wind load identification from structural responses. Engineering Structures, nº 27, pp. 1586-1598, 2005.

MACHADO, V. L. Bifurcações multiplas e comportamento não linear de sistemas dinâmicos. 1993. Tese (Doutorado) – Universidade Federal do Rio de Janeiro, COPPE, Rio de Janeiro, 1993.

MILLAR, Malcolm A. e BARGUIAN, Majid. Snap-through behaviour of cables in flexible structures. Computers & Structures, nº 77, pp. 361-369, 2000.

NAYFEH, A. H. e BALACHANDRAN, B. Applied nonlinear dynamics. Analitical, computational and experimental methods. New York, John Wiley & Sons, 1995.

NEVES, Francisco de Assis das. Vibrações de Estruturas Aporticadas Espaciais. 1990. 168 f. Tese (Doutorado) – Universidade Federal do Rio de Janeiro, COPPE, Rio de Janeiro, 1990.

OLIVEIRA, P. A. Análise estática não-linear de cabos suspensos utilizando o método dos elementos finitos. 2002. 93 f. Dissertação (Mestrado) – Universidade Federal do Paraná, Curitiba, 2002.

OLIVEIRA, P. A. et al. Análise estática não-linear de cabos utilizando o método dos elementos finitos. XXX Jornadas Sul-Americanas de Engenharia Estrutural, v. 1, pp. 1-19, 2002.

ORLANDO, Diego. Absorsor pendular para controle de vibrações de torres esbeltas. 2006. 168 f. Dissertação (Mestrado) – Pontifícia Universidade Católica do Rio de Janeiro, Rio de Janeiro, 2006.

PASQUETTI, Eduardo. Estabilidade estática e dinâmica de torres estaiadas. 2003. 99 f. Dissertação (Mestrado) – Pontifícia Universidade Católica do Rio de Janeiro, Rio de Janeiro, 2003.

PAZ, MARIO. Dinámica estructural. Barcelona, Editora Reverté, 1992. SANTEE, D. M. Vibrações não-lineares e instabilidades de elementos estruturais

sensíveis a imperfeições. 1999. Tese (Doutorado) – Pontifícia Universidade Católica do Rio de Janeiro, Rio de Janeiro, 1999.

SCHWEIZERHOF, K. H. e WRIGGERS, P. Consistent linearization for path following methods in nonlinear FE analysis. Computers Methods Applied Mechanics Engineering, v. 59, pp. 269-272, 1986.

SILVEIRA, R. A. M. Análise de elementos estruturais esbeltos com restrições unilaterais de contato. 1995. 212 f. Tese (Doutorado) – Pontifícia Universidade Católica do Rio de Janeiro, Rio de Janeiro, 1995.

SIMITSES, G. J. e HODGES, D. H. Fundamentals of structural stability. New York, Elsevier Science Publishing, 2006.

117

SOLIMAN, M. S. Fractal erosion of basins of attraction in coupled non-linear systems. Journal of Sound and Vibration, nº 182, v. 5, pp. 729-740, 1995.

THOMPSOM, J. M. T. e HUNT, G. W. Elastic instability phenomenal. New York, John Wiley & Sons, 1984.

THOMPSON, J. M. T. e STEWART, H. B. Nonlinear dynamics and chaos. Geometrical methods for engineers and scientists. Great Britain, John Wiley & Sons, 1993.

WAHBA, Y. M. F. et al. Dynamic response of guyed masts. Engineering Structures, v. 20, nº 12, pp. 1097-1101, 1998.b.

WAHBA, Y. M. F. et al. Evaluation of non-linear analysis of guyed antenna towers. Computers & Structures, nº 68, pp. 207-212, 1998.a.

WEI, P. et al. A catenary element for the analysis of cable structures. Applied Mathematics and Mechanics, v. 20, n° 5, pp. 532-534, 1999.

XU, Y. L. et al. Modal analysis of tower-cable system of Tsing Ma long suspension bridge. Engineering Structures, v. 19, nº 10, pp. 857-867, 1997.

YANG, Y. B e KUO, S. B. Theory and analysis of nonlinear framed structures. Prentice Hall, 1994. 580 f.

YANG, Y. B e SHIEH, M. S. Solution method for nonlinear problems with multiple critical points. AIAAJ, v. 28, n° 12, pp. 2110-2116, 1990.

YAN-LI, H. et al. Nonlinear discrete analysis method for random vibration of guyed masts under wind load. Journal of Wind Engineering and Industrial Aerodynamics, nº 91, pp. 513-525, 2003.

YAU, J. D. e YANG, Y. B. Vibration reduction for cable-stayed bridges traveled by high-speed trains. Finite Elements in Analysis and Design, nº 40, pp. 341-359, 2004.

118

7 APÊNDICE

Nesta seção apresenta-se a formulação do modelo discreto, desenvolvida neste

trabalho utilizando o software de análise simbólica MAPLE, a saber:

a) Energia potencial:

restart; with(linalg): with(student): with(plots): with(codegen,C):

Delta[pi]:=Delta[U]-Delta[W];

Delta[U]:=sum((1/2)*K[i]*Delta[L[i]]^2,i=1..nm);

Delta[L[i]]:=Delta[L[theta[i]]]+Delta[L[o[i]]];

Delta[L[o[i]]]:=F[o[i]]/K[i];

Delta[U]:=expand(Delta[U]);

Delta[W]:=p*yc*L*(1-cos(theta))+sum(P[i]*Gamma[i]*L*(1-cos(theta)),i=1..np);

# np -> n° de cargas externas concentradas ao longo da coluna;

# Gamma[i] -> Hi/L ... (posicao de apicacao da carga "i" na barra / comprimento da barra);

# y -> yb/L ... (cantro de gravidade da barra / comprimento da barra);

# p -> pesso proprio da barra;

# Pi -> carga externa concentrada "i" na barra;

Delta[pi];

nm:=2:

np:=1:

Delta[L[theta[1]]]:=gamma[1]*L*(sqrt(((1/tan(alpha[1]))+sin(theta))^2+cos(theta)^2)-

sqrt((1/tan(alpha[1])^2)+1));

# gamma[i] -> hi/L ... (posicao de fixacao da mola "i" na barra / comprmento da barra);

Delta[L[theta[2]]]:=gamma[1]*L*(sqrt(((1/tan(alpha[1]))-sin(theta))^2+cos(theta)^2)-

sqrt((1/tan(alpha[1])^2)+1));

# gamma[i] -> hi/L ... (posicao de fixacao da mola "i" na barra / comprmento da barra);

dPI:=diff(Delta[pi],theta);

119

b) Variáveis:

K[2]:=K[1]:

F[o[2]]:=F[o[1]]:

Rho:=7850:

Area:=1.914E-002:

Long:=100.0:

Fo1ca:=1000:

ha:=100:

peso:=0.0:

K1k:=230:

P1p:=P1p:

Hh:=100:

yy:=1/2:

Xi1:=(ha/Long)^2:

Xi0:=Fo1ca*(ha/Long)/(K1k*Long):

Xip:=peso*yy/(K1k*Long):

Lp1:=P1p*(Hh/Long)/(K1k*Long):

c) Equação do caminho pós-crítico dimensional:

solve(dPI,P[1]):

CPC_Dim:=expand(%);

Cpc:=powsubs(K[1]=K1k, L=Long, gamma[1]=ha/Long, Gamma[1]=Hh/Long, p=peso, yc=yy, F[o[1]]=Fo1ca,

alpha[1]=Pi*30/180, CPC_Dim):

plot(Cpc, theta=-Pi/9..Pi/9,thickness=2);

evalf(powsubs(theta=Pi*32/180,Cpc));

Load_cr:=limit(CPC_Dim, theta=0);

evalf(powsubs(K[1]=K1k, L=Long, gamma[1]=ha/Long, Gamma[1]=Hh/Long, p=peso, yc=yy, F[o[1]]=Fo1ca,

alpha[1]=Pi*32/180, Load_cr));

evalf(%*(Hh/Long)/(Long*K1k));

d) Equação do caminho pós-crítico adimensional:

expand(dPI/(L^2*K[1])):

powsubs(F[o[1]]/(L*K[1])=FoK, p/(L*K[1])=Pk, P[1]/(L*K[1])=Pa,%):

CPC_NDim:=solve(%, Pa);

evalf(powsubs( gamma[1]=ha/Long, Gamma[1]=Hh/Long, FoK=Fo1ca/(Long*K1k), Pk=peso/(Long*K1k),

alpha[1]=Pi*30/180, CPC_NDim));

plot(%, theta=-Pi/9..Pi/9,thickness=2);

e) Equação de movimento:

Eq_Dim:=expand(1/3*rho*A*L^3*diff(theta(t),t,t) + 2*1/3*rho*A*L^3*omega0*xi*diff(theta(t),t) + dPI);

120

f) Freqüência natural:

subs(xi=0, lambda_p1=Lp1, alpha1=Pi*32/180, A=Area, L=Long, rho=Rho, K1=K1k,

K[1]=K1k,gamma[1]=ha/Long, Gamma[1]=Hh/Long, F[o[1]]=0, P[1]=0, alpha[1]=Pi*32/180, p=0, yc=yy,

Eq_Dim):

powsubs(diff(theta(t),`$`(t,2))=uu, %):

evalf(isolate(%,uu)):

powsubs( theta=y1(t) , rhs(%)):

dsys1:=diff(y1(t),t)=y2(t), diff(y2(t),t)= % , y1(0)=10, y2(0)=0.0;

dsol1 := dsolve(dsys1, numeric,method=rkf45, output=listprocedure, maxfun=0, range=0..30):

odeplot(dsol1,[[t,y1(t)]],numpoints=1000);

dsol1y1 := subs(dsol1,y1(t)): dsol1y2 := subs(dsol1,y2(t)):

Period:=24.9:

dsol1y1(Period), dsol1y2(Period);

omega0:=evalf(2*Pi/Period);

g) Sistema de equações de primeira ordem (tempo adimensional):

expand(Eq_Dim/(L^2*K[1])):

powsubs(F[o[1]]/(L*K[1])=FoK, p/(L*K[1])=Pk, P[1]/(L*K[1])=Pa,%):

powsubs(diff(theta(t),`$`(t,2))=omega0^2*diff(theta(tau), tau, tau),

diff(theta(t),`$`(t,1))=omega0^2*diff(theta(tau), tau),%):

evalf(powsubs( gamma[1]=ha/Long, Gamma[1]=Hh/Long, FoK=Fo1ca/(Long*K1k), Pk=peso/(Long*K1k),

A=Area, L=Long, rho=Rho, K[1]=K1k, xi=0.01, alpha[1]=Pi*32/180, Pa=Ampl*cos(Omega*tau), %)):

eq_d1:=subs(tau=t, theta=delta1,%):

ee1:=powsubs(tt=diff(delta1(t),t$2), isolate(powsubs( diff(delta1(t),t$2)=tt,eq_d1),tt));

h) Sistema de equações em C++:

neq:=1: y:='y': count:=0: num:=0: xx:=1: temp:=0: d:=array(1..neq*2): for i from 1 to neq*2 do d[i]:=0 od:

for i from 1 to neq do

d[xx]:=F[count+1]=y[count+2]: xx:=xx+1:

dy:= coeff(powsubs(diff(delta||i(t),t)=uu, rhs(ee||i)),uu)*y[count+2]:

dir:=powsubs( diff(delta||i(t),t)=0,rhs(ee||i));

for j from 1 to nops(dir) do

temp:=temp + powsubs( seq(delta||k=y[2*k-1], k=1..neq), op(j, dir) ):

od;

d[xx]:=F[count+2]=dy+temp:

count:=count+2:

xx:=xx+1: temp:=0:

od:

C(d);

121

i) Solução do sistema utilizando Runge-Kutta:

Prim:=powsubs( diff(delta1(t),t)=y2(t) ,delta1=y1(t), Ampl=0.1, Omega=2.0, rhs(ee1) );

dsys1:=diff(y1(t),t)=y2(t), diff(y2(t),t)=Prim, y1(0)=0.01, y2(0)=0.0:

dsol1 := dsolve(dsys1, numeric,method=rkf45, output=listprocedure, maxfun=0, range=0..4000):

odeplot(dsol1,[[t,y1(t)]],numpoints=10000);

dsol1x := subs(dsol1,y1(t)): dsol1y := subs(dsol1,y2(t)):

dsol1x(20);

T_Total:=400.0: NDiv:=600: T_Step:=T_Total/NDiv;

j) Animação:

aa:=0:

for i from 1 to NDiv do

if (dsol1x(aa)>0) then

X[i]:=Long*cos(Pi/2-dsol1x(aa)): Y[i]:=Long*sin(Pi/2-dsol1x(aa)):

else

X[i]:=Long*cos(Pi/2-dsol1x(aa)): Y[i]:=Long*sin(Pi/2-dsol1x(aa)):

end if:

aa:=aa+T_Step: ti||i:=cat("t=",convert(aa,string)," seg"): od:

for jj from 1 to NDiv do

LL||jj:=PLOT( CURVES( [[0,0],[X[jj], Y[jj] ]]), THICKNESS(3), COLOR(RGB,1,0,0), TITLE(ti||jj) ):

Lc||jj:=PLOT( CURVES( [[60,0],[X[jj], Y[jj] ]]), CURVES( [[-60,0],[X[jj], Y[jj] ]]),THICKNESS(1),

COLOR(RGB,0,0,1) ):

Pl||jj:=display([LL||jj, Lc||jj]):

od:

display([seq(Pl||i, i=1..NDiv)], insequence=true);

Livros Grátis( http://www.livrosgratis.com.br )

Milhares de Livros para Download: Baixar livros de AdministraçãoBaixar livros de AgronomiaBaixar livros de ArquiteturaBaixar livros de ArtesBaixar livros de AstronomiaBaixar livros de Biologia GeralBaixar livros de Ciência da ComputaçãoBaixar livros de Ciência da InformaçãoBaixar livros de Ciência PolíticaBaixar livros de Ciências da SaúdeBaixar livros de ComunicaçãoBaixar livros do Conselho Nacional de Educação - CNEBaixar livros de Defesa civilBaixar livros de DireitoBaixar livros de Direitos humanosBaixar livros de EconomiaBaixar livros de Economia DomésticaBaixar livros de EducaçãoBaixar livros de Educação - TrânsitoBaixar livros de Educação FísicaBaixar livros de Engenharia AeroespacialBaixar livros de FarmáciaBaixar livros de FilosofiaBaixar livros de FísicaBaixar livros de GeociênciasBaixar livros de GeografiaBaixar livros de HistóriaBaixar livros de Línguas










http://www.livrosgratis.com.br/cat_1/administracao/1







http://www.livrosgratis.com.br/cat_2/agronomia/1







http://www.livrosgratis.com.br/cat_3/arquitetura/1







http://www.livrosgratis.com.br/cat_4/artes/1







http://www.livrosgratis.com.br/cat_5/astronomia/1







http://www.livrosgratis.com.br/cat_6/biologia_geral/1









http://www.livrosgratis.com.br/cat_8/ciencia_da_computacao/1











http://www.livrosgratis.com.br/cat_9/ciencia_da_informacao/1











http://www.livrosgratis.com.br/cat_7/ciencia_politica/1









http://www.livrosgratis.com.br/cat_10/ciencias_da_saude/1











http://www.livrosgratis.com.br/cat_11/comunicacao/1







http://www.livrosgratis.com.br/cat_12/conselho_nacional_de_educacao_-_cne/1















http://www.livrosgratis.com.br/cat_13/defesa_civil/1









http://www.livrosgratis.com.br/cat_14/direito/1







http://www.livrosgratis.com.br/cat_15/direitos_humanos/1









http://www.livrosgratis.com.br/cat_16/economia/1







http://www.livrosgratis.com.br/cat_17/economia_domestica/1









http://www.livrosgratis.com.br/cat_18/educacao/1







http://www.livrosgratis.com.br/cat_19/educacao_-_transito/1









http://www.livrosgratis.com.br/cat_20/educacao_fisica/1









http://www.livrosgratis.com.br/cat_21/engenharia_aeroespacial/1









http://www.livrosgratis.com.br/cat_22/farmacia/1







http://www.livrosgratis.com.br/cat_23/filosofia/1







http://www.livrosgratis.com.br/cat_24/fisica/1







http://www.livrosgratis.com.br/cat_25/geociencias/1







http://www.livrosgratis.com.br/cat_26/geografia/1







http://www.livrosgratis.com.br/cat_27/historia/1







http://www.livrosgratis.com.br/cat_31/linguas/1







Baixar livros de LiteraturaBaixar livros de Literatura de CordelBaixar livros de Literatura InfantilBaixar livros de MatemáticaBaixar livros de MedicinaBaixar livros de Medicina VeterináriaBaixar livros de Meio AmbienteBaixar livros de MeteorologiaBaixar Monografias e TCCBaixar livros MultidisciplinarBaixar livros de MúsicaBaixar livros de PsicologiaBaixar livros de QuímicaBaixar livros de Saúde ColetivaBaixar livros de Serviço SocialBaixar livros de SociologiaBaixar livros de TeologiaBaixar livros de TrabalhoBaixar livros de Turismo

http://www.livrosgratis.com.br/cat_28/literatura/1







http://www.livrosgratis.com.br/cat_30/literatura_de_cordel/1











http://www.livrosgratis.com.br/cat_29/literatura_infantil/1









http://www.livrosgratis.com.br/cat_32/matematica/1







http://www.livrosgratis.com.br/cat_33/medicina/1







http://www.livrosgratis.com.br/cat_34/medicina_veterinaria/1









http://www.livrosgratis.com.br/cat_35/meio_ambiente/1









http://www.livrosgratis.com.br/cat_36/meteorologia/1







http://www.livrosgratis.com.br/cat_45/monografias_e_tcc/1







http://www.livrosgratis.com.br/cat_37/multidisciplinar/1





http://www.livrosgratis.com.br/cat_38/musica/1







http://www.livrosgratis.com.br/cat_39/psicologia/1







http://www.livrosgratis.com.br/cat_40/quimica/1







http://www.livrosgratis.com.br/cat_41/saude_coletiva/1









http://www.livrosgratis.com.br/cat_42/servico_social/1









http://www.livrosgratis.com.br/cat_43/sociologia/1







http://www.livrosgratis.com.br/cat_44/teologia/1







http://www.livrosgratis.com.br/cat_46/trabalho/1







http://www.livrosgratis.com.br/cat_47/turismo/1







Documents

ANÁLISE DA INSTABILIDADE DINÂMICA DE ESTRUTURAS …livros01.livrosgratis.com.br/cp074018.pdf · CARVALHO, EULHER CHAVES (2008). Análise da Instabilidade Dinâmica de Estruturas