Embed Size (px)
Citation preview
UNIVERSIDADE ESTADUAL DO OESTE DO PARANÁ
CAMPUS DE FOZ DO IGUAÇU
PROGRAMA DE PÓS-GRADUAÇÃO EMENGENHARIA DE SISTEMAS DINÂMICOS E ENERGÉTICOS
DISSERTAÇÃO DE MESTRADO
ESTUDO NUMÉRICO E EXPERIMENTAL DADINÂMICA NÃO-LINEAR DE UM GIROSCÓPIO
ROSINEY DESIDÉRIO DA SILVA
FOZ DO IGUAÇU2012
Rosiney Desidério da Silva
Estudo Numérico e Experimental da Dinâmica Não-Linear deum Giroscópio
Dissertação de Mestrado apresentada ao Programa dePós-Graduação em Engenharia de Sistemas Dinâmi-cos e Energéticos como parte dos requisitos para ob-tenção do título de Mestre em Engenharia de Siste-mas Dinâmicos e Energéticos. Área de concentração:Sistemas Dinâmicos e Energéticos.
Orientador: Samuel da SilvaCo-orientador: Ubirajara Franco Moreno
Foz do Iguaçu2012
ii
Dados Internacionais de Catalogação na Publicação (CIP)
Biblioteca do Campus de Foz do Iguaçu – Unioeste Ficha catalográfica elaborada por Miriam Fenner R. Lucas - CRB-9/268
S586 Silva, Rosiney Desidério da Estudo numérico e experimental da dinâmica não-linear de um gi-
roscópio / Rosiney Desidério da Silva. - Foz do Iguaçu, 2012. 92 p. : il. : tab. : graf.
Orientador: Prof. Dr. Samuel da Silva. Co-orientador: Prof. Dr. Ubirajara Franco Moreno. Dissertação (Mestrado) – Programa de Pós-Graduação em Enge- nharia de Sistemas Dinâmicos e Energéticos - Universidade Estadual do Oeste do Paraná.
1. Giroscópio. 2. Sistemas dinâmicos não lineares. 3. Caos. 4. Expoentes de Lyapunov. 5. Séries temporais. I. Título.
CDU 629.1.054 517.93
iii
iv
Resumo
Este trabalho propõe um estudo da dinâmica de um giroscópio usando dados de simula-ção de um modelo analítico comparando com dados experimentais. Verifica-se a modelagemusando mecânica clássica, estudo de pontos de equilíbrio, estabilidade e verificação de regiõesonde o movimento do giroscópio pode ficar regular ou caótico. Os expoentes de Lyapunov sãoidentificados usando o método padrão, método de Eckmann-Ruelle, método de Wolf com sériestemporais e o teste 0-1. Os resultados alcançados nesta dissertação permitiram comparar asprincipais vantagens e desvantagens de cada um dos métodos e extrair informações qualitativase quantitativas sobre o movimento do giroscópio em estudo.
Palavras-chave: Giroscópio, Sistemas Dinâmicos Não-Lineares, Caos, Expoentes de Lyapu-nov, Séries Temporais.
v
Abstract
The present work proposes a study of the dynamics of a gyroscope using simulated dataof an analytical model by comparing with experimental data. Classical mechanical modelingapproaches are used to identify the equilibrium points, stability and verification of the regionswhere the motion equations of the gyroscope can present regular or chaotic behavior. The Lya-punov exponents are identified through the standard method, Eckmann-Ruelle Method, Wolfmethod with time series and the 0-1 test. The results achieved illustrate the main advantagesand drawbacks of each method and allow to observe qualitatively and quantitatively informationabout the motion of the gyroscope used.
Keywords: Gyroscope, Nonlinear Dynamical Systems, Chaos, Lyapunov Exponents, TimeSeries.
vi
Agradecimentos
Agradeço ao Professor Dr. Romeu Reginatto pelos conhecimentos repassados durante asaulas de Sistemas Dinâmicos Lineares e Sistemas Dinâmicos Não-Lineares, cujas discussõesforam de fundamental importância para o desenvolvimento desta dissertação.
Ao meu orientador, professor Dr. Samuel da Silva, agradeço pelas sugestões dadas du-rante o desenvolvimento do trabalho, sem as quais o mesmo não poderia ser finalizado.
Agradeço também as colaborações feitas pelo co-orientador, professor Dr. UbirajaraFranco Moreno, que contribuíram para elaboração do trabalho.
Agradeço aos professores Drs. Carlos Henrique Farias dos Santos e Daniel Iria Machado,pelas sugestões dadas durante o exame de qualificação.
Ao estudante de iniciação científica, Maurício Soares Leão, agradeço pela ajuda na reali-zação dos experimentos com o giroscópio que foi imprescindível para finalização do trabalho.
Agradeço às instituições:
Universidade Estadual do Oeste do Paraná (UNIOESTE), Campus de Foz do Iguaçu, portodo o apoio durante o período de realização do mestrado.
CAPES-PROAP por financiar a minha participação na Conferência Brasileira de Dinâ-mica e Controle (DINCON 2011, Águas de Lindóia, SP).
Fundação Parque Tecnológico Itaipu (FPTI) pela concessão da bolsa de estudos, possibi-litando assim dedicação exclusiva durante toda a realização da pesquisa.
Secretaria da Ciência, Tecnologia e Ensino Superior/Fundação Araucária (SETI-PR) pelacompra do giroscópio da Pascor utilizado neste trabalho.
E por fim, a todos aqueles que, de uma forma ou outra, contribuíram para a realizaçãodeste trabalho.
vii
viii
Sumário
Lista de Figuras xi
Lista de Tabelas xiv
Lista de Símbolos xv
1 Introdução 1
1.1 Objetivos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4
1.2 Contribuições do Trabalho . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4
1.3 Organização do Trabalho . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4
2 Revisão Bibliográfica 7
2.1 Modelagem e Solução Numérica . . . . . . . . . . . . . . . . . . . . . . . . . 7
2.2 Análise Não-Linear de Modelo de Giroscópios . . . . . . . . . . . . . . . . . 9
2.3 Análise Não-Linear Experimental . . . . . . . . . . . . . . . . . . . . . . . . 10
2.4 Considerações Finais . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
3 Modelagem do Giroscópio 15
3.1 Análise Cinemática . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
3.1.1 Sistema de Referência e Matrizes de Transformação de Coordenadas . 15
3.1.2 Vetores de Posição . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
3.1.3 Vetores de Velocidades e Acelerações Absolutas . . . . . . . . . . . . 20
3.2 Análise Dinâmica . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23
3.2.1 Diagrama de Corpo Livre . . . . . . . . . . . . . . . . . . . . . . . . 24
3.2.2 Tensor de Inércia . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28
ix
x
3.2.3 Torque Resultante . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31
3.2.4 Taxa de Variação do Momento Angular Resultante . . . . . . . . . . . 35
3.2.5 Equação de Euler . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38
3.2.6 Equação de Newton . . . . . . . . . . . . . . . . . . . . . . . . . . . 40
3.2.7 Equação da Trajetória na Base Inercial . . . . . . . . . . . . . . . . . . 41
3.3 Simulação e Comparação Experimental . . . . . . . . . . . . . . . . . . . . . 42
3.4 Considerações Finais . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46
4 Análise Dinâmica Não-Linear 49
4.1 Pontos de Equilíbrio e Estabilidade no Sentido de Lyapunov . . . . . . . . . . 49
4.2 Expoentes de Lyapunov . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50
4.2.1 Expoentes de Lyapunov do Modelo . . . . . . . . . . . . . . . . . . . 51
4.2.2 Expoentes de Lyapunov Experimental - Método de Wolf . . . . . . . . 55
4.3 Teste 0− 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57
4.4 Considerações Finais . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60
5 Análise da Dinâmica Não-Linear do Giroscópio 63
5.1 Pontos de equilíbrio e estabilidade . . . . . . . . . . . . . . . . . . . . . . . . 63
5.2 Expoentes de Lyapunov . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69
5.2.1 Expoentes de Lyapunov do modelo . . . . . . . . . . . . . . . . . . . 70
5.2.2 Expoentes de Lyapunov experimental . . . . . . . . . . . . . . . . . . 80
5.3 Teste 0− 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83
5.4 Conclusão . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84
6 Conclusões Finais e Trabalhos Futuros 87
Referências Bibliográficas 89
Apêndice 93
A Obtenção dos Parâmetros Indiretos do Modelo 93
xi
B Jacobiana da Equação Diferencial do Giroscópio 97
C Teste dos Programas de Análise Não-Linear 103
xii
Lista de Figuras
1.1 Giroscópio da Pascor, modelo ME-8960. . . . . . . . . . . . . . . . . . . . . 1
1.2 Movimentos de precessão, nutação e spin. . . . . . . . . . . . . . . . . . . . . 2
3.1 Sistema móvel B1. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
3.2 Sistema móvel B2. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
3.3 Ângulo de spin. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
3.4 Vetores de posição. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19
3.5 Força peso e força de reação. . . . . . . . . . . . . . . . . . . . . . . . . . . . 24
3.6 Força elástica. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26
3.7 Elementos para o cálculo do tensor de inércia. . . . . . . . . . . . . . . . . . . 29
3.8 Sensores de posição angular da Pascor, modelo CI-6625 (RMS). . . . . . . . . 43
3.9 Aquisição de dados. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44
3.10 Fio do sensor de nutação. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45
3.11 Posições angulares de precessão, nutação e spin. . . . . . . . . . . . . . . . . . 45
3.12 Velocidades angulares de precessão, nutação e spin. . . . . . . . . . . . . . . . 46
3.13 Trajetória do centro de massa do disco. . . . . . . . . . . . . . . . . . . . . . . 46
4.1 Exemplificação do método de Wolf. . . . . . . . . . . . . . . . . . . . . . . . 55
4.2 Gráfico p versus q para o mapa logístico. . . . . . . . . . . . . . . . . . . . . . 61
5.1 Exemplo de convergência de λ pelo método padrão. . . . . . . . . . . . . . . . 72
5.2 Exemplo de convergência de λ pelos métodos de Eckmann-Ruelle e Wolf. . . . 73
5.3 Expoentes de Lyapunov métodos padrão e de Eckmann-Ruelle variando bx. . . 74
5.4 Expoentes de Lyapunov métodos padrão e de Eckmann-Ruelle variando by. . . 74
xiii
xiv
5.5 Expoentes de Lyapunov métodos padrão e de Eckmann-Ruelle variando bz. . . 75
5.6 Soma dos expoentes de Lyapunov, métodos padrão e de Eckmann-Ruelle. . . . 75
5.7 Expoentes de Lyapunov com variação de h1. . . . . . . . . . . . . . . . . . . . 77
5.8 Gráfico p− q do teste 0− 1 com bx = 2,2× 10−4 N ·m · s/rad. . . . . . . . . 78
5.9 Gráfico p− q do teste 0− 1 com bx = 7,3× 10−4 N ·m · s/rad. . . . . . . . . 78
5.10 Espaço de fases de x2, x3 e x4 com bx = 2,2× 10−4 N ·m · s/rad. . . . . . . . 79
5.11 Espaço de fases de x2, x3 e x4 com bx = 7,3× 10−4 N ·m · s/rad. . . . . . . . 79
5.12 Exemplo de convergência de λ da série temporal experimental de x3. . . . . . . 81
5.13 Expoentes de Lyapunov da série temporal experimental de x3. . . . . . . . . . 81
5.14 Exemplo de convergência de λ da série temporal simulada de x3 . . . . . . . . 82
5.15 Expoentes de Lyapunov da série temporal simulada de x3. . . . . . . . . . . . 83
5.16 Gráfico p− q do teste 0− 1 com h1 = −0,1979 m. . . . . . . . . . . . . . . . 85
A.1 Comparação dos dados experimentais e simulados da precessão após ajuste. . . 95
A.2 Comparação dos dados experimentais e simulados da nutação após ajuste. . . . 95
C.1 Maior expoente de Lyapunov para o sistema de Lorenz. . . . . . . . . . . . . . 104
C.2 Atrator do sistema de Lorenz com comportamento regular. . . . . . . . . . . . 105
C.3 Atrator do sistema de Lorenz com comportamento caótico. . . . . . . . . . . . 106
Lista de Tabelas
3.1 Parâmetros dos contrapesos, da haste e do disco do giroscópio da Pascor. . . . 42
3.2 Condições iniciais para validação do modelo do giroscópio da Pascor. . . . . . 45
5.1 Sinais de γ dos pontos de equilíbrio do modelo. . . . . . . . . . . . . . . . . . 68
5.2 Condições iniciais para estimativa dos expoentes de Lyapunov. . . . . . . . . . 70
5.3 Resumo do teste 0− 1 aplicado a duas séries temporais simulada de x3. . . . . 78
5.4 Resumo da série temporal experimental de x3 e de λl e λw. . . . . . . . . . . . 80
5.5 Resumo da série temporal simulada de x3 e de λl e λw. . . . . . . . . . . . . . 82
5.6 Resumo da série temporal experimental de x3 e do teste 0− 1. . . . . . . . . . 84
5.7 Resumo da série temporal simulada de x3 e do teste 0− 1 . . . . . . . . . . . . 85
A.1 Parâmetros indiretos do giroscópio da Pascor. . . . . . . . . . . . . . . . . . . 96
C.1 Resumo da série temporal simulada de x e de λ do sistema de Lorentz. . . . . . 106
xv
xvi
Lista de Símbolos
PGESDE Programa de Pós-Graduação em Engenharia de Sistemas Dinâmicos e Ener-géticos
UNIOESTE Universidade Estadual do Oeste do ParanáEDO Equação Diferencial Ordinária
IS, B1S e B2S Representação geral de vetores do tipo posição, força e torque nas basesinercial, móvel B1 e móvel B2 respectivamente
X, Y e Z Eixos do sistema inercialX1, Y1 e Z1 Eixos do sistema móvel B1
X2, Y2 e Z2 Eixos do sistema móvel B2
O, O1 e O2 Origens dos sistemas inercial, móvel B1 e móvel B2 respectivamenteT φ e T T
φ Matriz de transformação da base I para a base B1 e transformação inversarespectivamente
T θ e T Tθ Matriz de transformação da base B1 para a base B2 e transformação inversa
respectivamente
IΦ, B1Φ e B2Φ Velocidades angulares de precessão nos sistemas inercial, móvel B1 e móvelB2 respectivamente
B1Θ e B2Θ Velocidades angulares de nutação nos sistemas móveis B1 e B2
B2Ψ Velocidade angular de spin na base móvel B2
X Vetor de estados da EDO de dimensão n cujas componentes são xi com i
variando de 1 a nF Função da EDOXe Ponto de equilíbrio da EDO FJ Matriz Jacobiana da EDO Fφ, x1 Posição angular de precessãoθ, x3 Posição angular de nutaçãoψ, x5 Posição angular de spinφ, x2 Velocidade angular de precessãoθ, x4 Velocidade angular de nutaçãoψ, x6 Velocidade angular de spin
B2Lr Taxa de variação do momento angular resultante
xvii
xviii
B2τ r Torque resultanteτx, τy e τz Componentes do torque resultanteτey (θ) Componente do torque elástico na direção de 2m1, m2, md e mh Massas do contrapeso maior, contrapeso menor, disco e haste respectiva-
menteh1, h2, hd e hh Posições do centro de massa do contrapeso maior, contrapeso menor, disco
e haste respectivamenteR1, R2, Rd e Rh Raios do contrapeso maior, contrapeso menor, disco e haste respectivamentel1, l2 e mh Larguras do contrapeso maior, contrapeso menor e haste respectivamenteM Massa total do sistemaH Componente na direção ı2 da posição do centro de massa do sistemag Aceleração da gravidadeθi, θf Limites inferior e superior do ângulo de nutação respectivamenteTx, Ty e Tz Coeficientes do torque cinéticobx, by e bz Coeficientes do torque viscosok Coeficiente do torque elásticoIxx, Iyy e Izz Componentes diagonais do tensor de inércia do discoIx, Iy e Iz Componentes diagonais do tensor de inércia resultanteγ Autovalorλ Expoente de Lyapunovxi, x[i], x(iδt) Representa a série temporalδt Tempo de amostragem de uma série temporalN Número de pontos de uma série temporalN
′ Número de pontos de um atrator reconstruído∆t Número de evoluções (tempo discreto) ou passo (tempo contínuo)m Dimensão de imersão do atrator reconstruídoτ Tempo de atraso de Takens (tempo discreto)ν Número total de passos de substituição
Capítulo 1
Introdução
O giroscópio mecânico é um dispositivo composto de um rotor ou mais, que gira emtorno de um eixo de simetria com liberdade de mudança de direção. Um giroscópio em equilí-brio apresenta a propriedade de manter a direção do eixo de rotação com base no princípio deconservação do momento angular. Quando se aplica uma força externa no giroscópio, causandoum torque, ele tende a se mover de forma não-intuitiva em uma direção perpendicular, tanto aoeixo de rotação do disco, quanto à direção dessa força. Um exemplo de giroscópio mecânicousado em ensino ou demonstração pode ser visto na figura 1.1, que é o modelo usado nestetrabalho.
Figura 1.1: Giroscópio da Pascor, modelo ME-8960.
Os créditos pela descoberta do efeito giroscópico geralmente são dados a Foucault, pela
1
2
investigação do movimento de rotação da Terra em 1852 (Titterton and Weston, 2004; Acarand Shkel, 2008). O efeito giroscópico, pode aparecer em dispositivos mecânicos diversos, porexemplo pião e roda de bicicleta.
No giroscópio há três movimentos de rotação, que são exemplificados na figura 1.2 econhecidos como (Crabtree, 1909):
Spin: Há um eixo pelo qual o rotor do giroscópio pode girar. O movimento de rotação emtorno desse eixo chama-se spin ou rotação própria (R na figura);
Precessão: Movimento causado quando existe um torque aplicado que causa uma variação dadireção do vetor momento angular. O vetor momento angular sempre terá a direção dovetor velocidade angular e neste caso, pela equação do movimento de Euler, uma reaçãocausará o movimento de precessão (P na figura);
Nutação: No movimento de precessão, pode ocorrer um movimento oscilatório, tal movimentoé a nutação (N na figura).
Figura 1.2: Movimentos de precessão (P), nutação (N) e spin (S). Fonte: (Angular, 2012).
Existem outros dispositivos que possuem o mesmo comportamento do giroscópio mecâ-nico descrito acima, mas com um princípio de funcionamento diferente. Apesar disso, estes dis-positivos também são chamados de giroscópios. Como exemplo de giroscópios temos (Tittertonand Weston, 2004): de fibra ótica, com anel laser e MEMS (MicroElectroMechanical Sys-tems) (Acar and Shkel, 2008), entre muitos outros.
Algumas aplicações práticas do efeito giroscópico em sistemas inerciais incluem (Tittertonand Weston, 2004):
3
Estabilização: utilizada para evitar, por exemplo, que as imagens gravadas por uma filmadorafiquem trêmulas;
Piloto automático: utilizado para manter, por exemplo, um navio ou avião em determinadarota;
Navegação: substituindo a bússola como referência de direção.
Na navegação por inércia, giroscópios e acelerômetros (sensores de inércia) são usa-dos para detectar o movimento de rotação e de translação em relação a um referencial iner-cial (Titterton and Weston, 2004; Lawrence, 1998).
Equipamentos que possuem partes girantes, por exemplo, discos excêntricos, mecanismobiela-manivela entre outros, também podem apresentar o efeito giroscópico na forma de movi-mento, como o de precessão. Muitas vezes tais movimentos não são desejados e nem levadosem consideração durante o projeto. Porém, em muitas situações, tais efeitos podem originarforças excessivas e não previstas inicialmente, provocando desgastes que levem à falha doscomponentes do sistema.
Assim, o entendimento do comportamento dinâmico do efeito giroscópico em sistemasé desejado e necessário. Entre as ferramentas para se fazer isto, a análise não-linear de umgiroscópio é útil para a verificação da estabilidade do sistema e também verificar se o sistemapode ter um comportamento caótico. Um sistema apresenta comportamento caótico quandoé sensível às condições iniciais, ou seja, para uma condição inicial com erro, a solução dosistema para um dado tempo apresenta um valor completamente diferente daquele, caso nãohouvesse erro. Isso faz com o comportamento do sistema, para longos períodos de tempo, setorne imprevisível.
De forma muito resumida a Teoria de Sistemas Dinâmicos teve basicamente seu iní-cio com Newton através de suas leis da mecânica, juntamente com o desenvolvimento docálculo diferencial e integral, desenvolvido também de forma independente por G. W. Leib-niz. Com essas ferramentas é possível obter matematicamente as trajetórias de corpos no es-paço (Monteiro, 2006; Strogatz, 1994).
Uma parte importante da análise de equações diferenciais não-lineares é o estudo de es-tabilidade. Contribuições a respeito de estabilidade de sistemas dinâmicos foram feitas pormuitos cientistas, podendo-se citar: Laplace, Lagrange, Poisson, Haret, Lyapunov, Poincaré,Kolmogorov, Arnold, Moser, dentre outros (Monteiro, 2006; Strogatz, 1994).
Um marco no estudo de sistemas não-lineares foi a publicação de Lorenz em 1963, deum artigo intitulado Deterministic Nonperiodic Flow (Lorenz, 1963), com um sistema simpli-ficado com apenas três equações diferenciais não-lineares de primeira ordem que apresentavacomportamento caótico. Este trabalho comprovou muitas das ideias lançadas por Poincaré anosantes, mas não estudadas na época por falta de recursos computacionais que permitissem a so-
4
lução numérica de equações diferenciais para uma gama ampla de condições iniciais. A partirdaí, o interesse em estudos caóticos e suas implicações nos ramos diversos da ciência têm sidoamplamente abordada na literatura.
1.1 Objetivos
O objetivo desta dissertação é estudar e analisar a estabilidade dos pontos de equilíbrio eos expoentes de Lyapunov de um modelo de giroscópio bem como o uso de séries temporais.O modelo matemático do giroscópio utilizado neste trabalho é obtido e validado com dadosexperimentais. O trabalho apresenta uma comparação entre métodos de obtenção de expoentesde Lyapunov com relação ao desempenho e complexidade para o caso do giroscópio em estudo.
1.2 Contribuições do Trabalho
Com base na realização deste trabalho, observa-se que esta dissertação contribuirá comos tópicos:
Análise não-linear do modelo: os trabalhos anteriores que modelaram giroscópios equivalen-tes ao da Pascor (Aballay and Avilés, 2002), (Kostov and Hammer, 2010) e (McGlynn,2007) assumem muitas simplificações, como desconsiderar o atrito, que podem compro-meter a análise do movimento. Esta dissertação, ao contrário, emprega um modelo usandoforças de atrito diversas, e assumindo várias hipóteses que permitem que regimes dinâmi-cos bem mais ricos e complexos sejam representados e comprovados experimentalmente.A partir deste modelo proposto, uma análise completa envolvendo estudo de estabilidadee expoentes de Lyapunov, assim como a utilização do teste 0− 1, é feita. Algumas destasanálises numéricas foram também correlacionadas com os testes experimentais;
Análise não-linear experimental: poucos trabalhos realizaram análise não-linear usando da-dos experimentais de operação. Como este trabalho irá empregar o giroscópio Pascor
existente no laboratório e medir sinais temporais de posição de precessão e nutação,espera-se com os resultados contribuir com análise quantitativa destes aspectos a partirdos sinais. Em especial, com a obtenção dos expoentes de Lyapunov experimentalmente,assim como realização do teste 0−1, e verificar regiões de movimento caótico ou regular.
1.3 Organização do Trabalho
Este trabalho está organizado em seis capítulos e três apêndices:
5
Capítulo 2 - Revisão Bibliográfica: Uma revisão bibliográfica simples é feita mostrando oatual estágio de desenvolvimento do tema;
Capítulo 3 - Modelagem do Giroscópio: Neste capítulo são feitas todas as considerações ne-cessárias para a modelagem do giroscópio de forma mais completa visando a obtençãodas equações diferenciais do movimento e a sua validação a partir de dados experimen-tais;
Capítulo 4 - Análise Dinâmica Não-Linear: Nesta parte é feita uma revisão dos conceitos bá-sicos de sistemas dinâmicos não-lineares, com ênfase nos métodos para estimativa dosexpoentes de Lyapunov;
Capítulo 5 - Análise da Dinâmica Não-Linear do Giroscópio: Tendo como base os concei-tos do capítulo 4 é feita uma análise não-linear tanto do modelo matemático quanto dosresultados experimentais do giroscópio da Pascor;
Capítulo 6 - Conclusões Finais e Trabalhos Futuros: As considerações finais com sugestõespara trabalhos futuros são apresentados neste capítulo;
Apêndice A - Obtenção dos Parâmetros Indiretos do Modelo: Os parâmetros das equaçõesdiferenciais do giroscópio obtidas no capítulo 3 são estimados experimentalmente nesteapêndice;
Apêndice B - Jacobiana da Equação Diferencial do Giroscópio: Neste apêndice a matriz ja-cobiana do sistema giroscópico é obtida. Esta matriz é usada no estudo de estabilidadedo sistema e também é utilizada nos métodos de obtenção de Expoentes de Lyapunovrelacionados a modelos;
Apêndice C - Teste dos Programas de Análise Não-Linear: Neste apêndice os programas con-feccionados e usados no capítulo 5 são testados para mostrar que estão corretos. No caso,o teste é feito com o sistema não-linear de Lorentz por ser clássico, bem conhecido e defácil comparação.
6
Capítulo 2
Revisão Bibliográfica
O estudo de sistemas físicos pode ser dividido de forma didática em quatro partes:
• Modelagem usando leis fundamentais;
• Solução do modelo analítico, quando possível ou numérico;
• Análise do modelo (linear ou não-linear);
• Análise experimental (linear ou não-linear).
Um bom modelo de um sistema físico, em tese, é capaz de predizer os resultados expe-rimentais, e, normalmente, quanto mais elaborado for o modelo, melhor são as predições. Poroutro lado, o modelo poderia se tornar demasiadamente complicado. Portanto, deve-se esco-lher um meio termo, ou seja, procurar ser simples, mas representativo. Neste sentido, muitossistemas não-lineares podem ser linearizados para funcionarem em uma certa faixa de operação.
Este trabalho trata do estudo numérico e experimental do comportamento dinâmico deum giroscópio, considerando aspectos não-lineares. Assim, esta revisão busca apresentar umpanorama da modelagem e análise de sistemas giroscópicos, além de formas gerais para carac-terização de sistemas não-lineares.
As subseções a seguir foram feitas de forma a separar os conteúdos dos diversos trabalhosusados nessa revisão, em três partes:
• Modelagem e solução numérica dos modelos de giroscópios (analítica quando possível);
• Análise não-linear de modelo;
• Análise não-linear experimental.
2.1 Modelagem e Solução Numérica
Muitos modelos simplificados de giroscópios mostram as relações entre as grandezas en-volvidas sem a utilização de equações diferenciais. Já outros modelos, mais sofisticados, podem
7
8
levar em consideração efeitos dissipativos. Esta seção mostra a descrição de alguns estudos fei-tos sobre giroscópios diversos.
Nos artigos (Aballay and Avilés, 2002), (Kostov and Hammer, 2010) e (McGlynn, 2007)são obtidos modelos de um giroscópio com as mesmas características do giroscópio Pascor
ME-8960, usado neste trabalho. Este giroscópio é composto por uma haste com dois graus deliberdade (precessão e nutação) e um disco que pode girar em torno dessa haste (spin). A mo-delagem feita por (Aballay and Avilés, 2002) usou conceitos de física básica desconsiderandoo movimento de nutação e assumindo apenas as relações de velocidade entre spin e precessão,portanto, não foi obtida a equação diferencial do movimento. Já no artigo (Kostov and Ham-mer, 2010) considerou-se a velocidade de nutação. Nestes dois casos foram feitos experimentosque validam todas as equações deduzidas. Já (McGlynn, 2007) faz uma modelagem mais elabo-rada, mas ainda assim com várias restrições, como desconsiderar o atrito e supor velocidades deprecessão e spin constantes. Em consequência das hipóteses assumidas e utilizando o métodode Newton-Euler chegou a uma solução analítica para o problema, da qual obteve o resultadoque o ângulo de nutação é constante.
Com o objetivo de estudar um giroscópio não axis-simétrico conservativo, com três grausde liberdade de rotação, com enfase no movimento com oscilações de grandes amplitude, foiusado no estudo em questão ângulos cardânicos e o conceito de quatérnions (Carrera and Weber,2009).
Um caso mais sofisticado que envolveu a modelagem de um giroscópio simétrico, mon-tado sobre uma base vibratória com movimento harmônico, com forças dissipativas propor-cionais a velocidade, foram feitas nos trabalhos (Chen and Ge, 2005) e (Ge et al., 1996). Omovimento foi descrito por ângulos de Euler e o modelo foi deduzido usando a formulaçãoLagrangiana da mecânica. No primeiro artigo, (Chen and Ge, 2005), o sistema foi consideradocom dois graus de liberdade, com um amortecedor de vibrações mecânicas fixado no interiorao longo do eixo de rotação central, sob a forma de um sistema massa-mola-amortecedor. Jáem (Ge et al., 1996) o sistema foi considerado com apenas um grau de liberdade e a forçadissipativa foi viscosa.
Outro estudo de giroscópio que leva em consideração o efeito do atrito é apresentadoem (Carrera et al., 2010). As forças de atrito foram causadas pelos rolamentos que compunhamas partes móveis do giroscópio. Foi assumido por hipótese que o atrito era viscoso e no modeloforam considerados os momentos de inércia dos balancins e do rotor. Para obter as equações demovimento foram utilizadas as leis de Newton-Euler e empregados ângulos cardânicos. Parafacilitar a análise do modelo usou-se o processo de normalização com o objetivo de reduzir umparâmetro na equação.
Uma etapa importante na análise de sistemas dinâmicos é a solução numérica das equa-ções do movimento. Para obter a solução numérica de uma equação diferencial pode-se em-pregar algum método da família de Runge-Kutta, como utilizado nos artigos (Ge et al., 1996),
9
(Ramasubramanian and Sriram, 2000). Uma referência de como o método de Runge-Kutta fun-ciona, além de vários outros métodos numéricos de resolução de equações diferenciais, podeser visto em (Hairer et al., 1993) e (Hairer and Wanner, 1991).
2.2 Análise Não-Linear de Modelo de Giroscópios
Diversas ferramentas matemáticas podem ser usadas para o estudo de sistemas não-lineares.Nos trabalhos (Chen and Ge, 2005), (Ge et al., 1996) utilizou-se: diagramas de bifurcação, re-tratos de fases, mapas de Poincaré e o espectro dos expoentes de Lyapunov usando o métodopara fluxos que é aplicado quando se tem conhecimento das equações diferenciais do sistema, ométodo para fluxos é descrito em (Wolf et al., 1985). Com a ajuda destas ferramentas foi verifi-cado que os sistemas giroscópicos em questão apresentavam tanto movimentos regulares quantocaóticos conforme os valores dos parâmetros. Em (Chen and Ge, 2005) foi obtido também asdimensões de Lyapunov e foi feito um estudo qualitativo usando o teorema da variedade centrale o teorema da forma normal. Já em (Ge et al., 1996) a estabilidade do sistema foi verificadacom a teoria de Mathieu e o método direto de Lyapunov; como há um distúrbio externo (basevibrando) é usado o método analítico de Melnikov que serve para provar a existência de bifur-cações homoclínicas ou heteroclínicas ajudando verificar a existência de caos no movimento; oespectro de potência média é encontrado e também usado para verificar a existência de caos nomovimento.
Os expoentes de Lyapunov fornecem uma caracterização qualitativa e quantitativa docomportamento dinâmico de um sistema e estão relacionados diretamente com as taxas médiasexponenciais de divergência ou convergência de órbitas próximas no espaço de fases. Qualquersistema com pelo menos um expoente de Lyapunov positivo é definido como sendo caótico.Para um sistema dinâmico de dimensão n, têm-se n expoentes de Lyapunov. Essas proprieda-des são descritas em (Wolf et al., 1985) e (Ramasubramanian and Sriram, 2000). Para sistemasem que se conhece as equações diferenciais do movimento, pode ser utilizado para o cálculodo espectro de Lyapunov o método descrito em (Wolf et al., 1985), conhecido como métodopadrão.
Em (Ramasubramanian and Sriram, 2000) foi feito uma comparação numérica entre trêsalgoritmos para o cálculo do espectro completo de Lyapunov. Estes algoritmos foram testadospara alguns sistemas não-lineares com 2, 3 e 4 variáveis no espaço de estados usando o mé-todo de Runge-Kutta com passo variável para resolver as equações diferenciais. Os algoritmostestados foram: o método padrão, uma formulação diferencial do método padrão, e um novoalgoritmo que não exigia a ortonormalização. Houve uma concordância razoável entre os resul-tados obtidos utilizando os três algoritmos na maioria dos casos. Porém, considerando o tempode processamento para o cálculo dos espectros de Lyapunov, o método padrão foi mais eficienteseguido pelo novo método e pela versão diferencial do método padrão. Apesar da formulação
10
diferencial do método padrão e o novo algoritmo serem menos eficientes, eles ainda são úteiscomo algoritmos alternativos para o cálculo do espectro de Lyapunov.
No método padrão, que é descrito em (Wolf et al., 1985), os expoentes são obtidos atravésda comparação entre a trajetória obtida pela solução da equação diferencial não-linear e a traje-tória obtida pela solução linearizada da equação original. Desta forma, são feitas aproximaçõessucessivas e, após um tempo de evolução da solução do sistema, é feito um processo de ortonor-malização das soluções linearizadas. No método padrão é necessária a matriz Jacobiana, que éusada para obter as versões linearizadas da EDO original.
Outro método de obtenção dos expoentes de Lyapunov, que utiliza as versões lineariza-das da EDO original é descrito no método de Eckmann-Ruelle (Eckmann and Ruelle, 1985).As aproximações para os expoentes de Lyapunov são feitas através dos módulos dos autovalo-res dos produtos da matriz Jacobiana e, por isso, este método acaba exigindo mais tempo deprocessamento em relação ao método padrão.
Um método para obter os expoentes de Lyapunov que não utiliza linearização das equa-ções diferenciais originais do sistema pode ser visto em (Fazanaro et al., 2010). Tal método,que é conhecido como o método via dinâmicas clonadas, utiliza, no lugar das equações linea-rizadas, cópias (clones) da equação original. O algoritmo utilizado para obter os expoentes deLyapunov, para o método em questão, é o mesmo do método padrão com a diferença de queas n equações linearizadas são substituídas por cópias da equação diferencial original, cujas ncondições iniciais ortonormais diferem de ε da condição inicial da trajetória fiducial, que é atrajetória original obtida a partir da solução da equação diferencial não-linear a partir das con-dições iniciais. O método via dinâmicas clonadas tem vantagens na aplicação em sistemas nãosuaves, já que neste caso não é necessário a linearização que, para tais sistemas, pode apre-sentar complicações devido a descontinuidade das equações linearizadas. No teste do métodoem questão é utilizado um sistema clássico, no caso o oscilador de Chua, mostrando resultadossatisfatórios.
2.3 Análise Não-Linear Experimental
Para um dado sistema dinâmico, nem sempre é possível ter acesso aos dados experimen-tais de todas as grandezas físicas que o caracterizam de forma simultânea. Porém, mesmoassim, em geral é possível, através de algumas dessas grandezas, na forma de uma série tempo-ral 1, caracterizar o sistema em sua totalidade. No geral, a série temporal de uma das grandezascontém informações das outras grandezas que, de alguma forma, não se tem acesso. Portanto,a partir de uma série temporal de uma grandeza pode ser possível reconstruir a trajetória noespaço de fases do sistema completo. Isto significa, por exemplo, reconstruir o atrator, a partir
1Definida como sucessivos valores de um dado estado com intervalo regular de tempo.
11
do que podem ser obtidas certas quantidades invariantes, como os expoentes de Lyapunov e adimensão do atrator.
O atrator reconstruído é topologicamente equivalente ao atrator que seria produzido casose conhecesse a dinâmica do sistema através da equação diferencial. Também vale ressaltarque a dimensão e os expoentes de Lyapunov do atrator reconstruído são aproximadamente osmesmos do atrator original. Um método bastante popular para a reconstrução de atratores é ométodo das coordenadas de atraso temporal onde os pontos do espaço de fases estão na forma
Si =x(ti) x(ti + τ) . . . x(ti + (m− 1)τ)
T, onde x(ti) representa a série temporal, τ é
o tempo de atraso de Takens e m é a dimensão de imersão do atrator reconstruído (Campanharoet al., 2005).
Há alguns métodos para estimar o tempo de atraso τ , por exemplo, a função de auto-correlação e o conceito de informação mútua (Perc, 2006; Kodba et al., 2005; Small, 2005;Addison, 1997; Kantz and Schreiber, 2004). Para séries advindas de um sistema caótico afunção de auto-correlação não é o método mais adequado para encontrar o tempo de atraso,pois este não leva em conta as correlações não-lineares. Já o mínimo da informação mútua éum critério mais adequado em relação a função de auto-correlação quando aplicado a sériestemporais não-lineares (Fraser and Swinney, 1986; Perc, 2006).
Para a determinação do valor mínimo da dimensão de imersão m pode-se usar o métododos falsos vizinhos próximos (Perc, 2006; Kodba et al., 2005; Small, 2005; Kantz and Schreiber,2004).
Em (Kodba et al., 2005) é feito um estudo de dados experimentais de um circuito resso-nante resistor-indutor-diodo em série, alimentado por uma tensão de entrada senoidal, que é oparâmetro de estudo, e a saída do circuito é a tensão sobre o indutor-diodo, que forma os dadosda série temporal. Já em (Perc, 2006) o estudo feito envolve uma série temporal obtida de umadas variáveis de estado do sistema de Lorenz simulado, ou seja, pela resolução da EDO. Nostrabalhos (Kodba et al., 2005; Perc, 2006) a partir das séries temporais os atratores são recons-truídos através do método das coordenadas de atraso temporal. Primeiramente o tempo de atrasoé encontrado pelo primeiro mínimo da informação mútua, em seguida a dimensão de imersão éencontrada pelo método dos falsos vizinhos próximos. Com os atratores reconstruídos, verifica-se se os sistemas são determinísticos. Caso não sejam, os sistemas são aleatórios. Neste caso osinvariantes obtidos do atrator reconstruído não apresentariam informações válidas. Na sequên-cia, calculou-se o maior expoente de Lyapunov para cada sistema, através do método de Wolf,onde foram encontrados expoentes positivos indicando comportamento caótico.
No método de Wolf (Wolf et al., 1985), os expoentes de Lyapunov são obtidos da análisedo atrator reconstruído a partir dos dados da série temporal de uma das variáveis de estadodo sistema em estudo. Tal método é baseado nas ideias do método padrão, que é utilizadoquando se conhece as equações diferenciais do sistema. A abordagem do método de Wolfimplementa o conceito de expoentes de Lyapunov de forma simples e bastante direta (Kodba
12
et al., 2005; Perc, 2006).
Em (Gottwald and Melbourne, 2004) foi apresentado um novo teste para verificar seum sistema dinâmico determinístico apresenta comportamento regular (periódico ou quase-periódico) ou caótico, o qual foi denominado “teste 0 − 1”. O método se aplica diretamentea dados de uma série temporal, logo não tem-se a necessidade de reconstrução do espaço defases. Como resultado do método, tem-se o valor K = 0 se o sistema dinâmico é regular ouK = 1 se for caótico, ou seja, a saída do método é binária. Na prática os valores de saída dométodo são próximos de 0 ou de 1, pois a convergência depende, por exemplo, da quantidadede dados da série temporal. O teste 0− 1 tem implementação computacional fácil e o tempo deexecução baixo.
Enquanto (Gottwald and Melbourne, 2004) trata de sistemas dinâmicos determinísticossem ruído, em (Gottwald and Melbourne, 2005) o ruído, de forma moderada, foi levado emconsideração na aplicação do teste 0− 1. No artigo em questão foram feitas algumas modifica-ções no teste original e comparações com os expoentes de Lyapunov foram feitas por meio desimulações numéricas de séries temporais de alguns sistemas determinísticos clássicos.
Em (Falconer et al., 2007), para verificar a eficacia do teste 0 − 1 em relação a dadosexperimentais, executa-se um experimento com base em um motor com dois pólos. Foi veri-ficado para um conjunto de dados do motor a uma frequência de alimentação ω = 0,9 Hz queao aplicar o teste o resultado obtido foi de K = 0,02 indicando dinâmica regular já para o con-junto de dados cuja frequência ω = 0,6 Hz tem-se K = 0,92 que apresenta dinâmica caótica.Obteve-se um comportamento melhor com o aumento de N . Concluiu-se que quanto menorfosse a quantidade de ruído da série temporal, maior seria a convergência do valor de K. Combase nos valores de K obtidos, verificou-se que o teste é robusto à contaminação dos sinais peloruído.
As justificativas teóricas do teste 0 − 1 podem ser vista em (Gottwald and Melbourne,2009b). No trabalho citado uma versão simplificada do teste foi verificada de forma rigorosa e,com a análise do teste, propostas de melhorias foram feitas.
Uma discussão sobre a implementação do teste 0 − 1 foi feita em (Gottwald and Mel-bourne, 2009a) com base nas melhorias feitas nos artigos anteriores, como o artigo sobre ajustificativa teórica do teste (Gottwald and Melbourne, 2009b).
Em um artigo mais recente (Ke-Hui et al., 2010), foram feitas simulações numéricas dediversos sistemas dinâmicos não-lineares, inclusive de sistemas de ordem fracionária 2, e apre-sentada uma comparação do teste 0 − 1 com o valor máximo do expoente de Lyapunov. Ossistemas dinâmicos usados no estudo foram o mapa de Hénon, sistema de Lorentz simplificado,de ordem inteira e de ordem fracionária. A conclusão foi que o teste funciona bem apesar de,
2Um sistema dinâmico é de ordem fracionária quando a equação diferencial que o representa contém derivadasde ordem fracionária.
13
eventualmente, ter zonas cinzentas 3, que podem ocorrer, por exemplo, caso não se tenha umasérie temporal com tamanho suficiente.
2.4 Considerações Finais
Para o modelo de giroscópio da Pascor, adotado neste trabalho, as modelagens matemá-ticas encontradas na literatura foram simples, não levando em consideração certas propriedadescomo, por exemplo, forças de atrito, que causam a diminuição dos movimentos de precessão,nutação e spin. Nestes casos nenhuma análise não-linear foi feita.
Para os modelos matemáticos mais elaborados dos giroscópios, que são diferentes aomodelo da Pascor, foram feitas diversas análises não-lineares com o objetivo de verificar aestabilidade e regiões onde o comportamento é regular ou caótico. Não há estudos de resultadosexperimentais em relação a análise não-linear, ou seja, foi dado ênfase na análise dos modelos.
Na análise não-linear dos modelos matemáticos dos giroscópios foram vistos: estudode estabilidade de pontos de equilíbrio através do teorema da variedade central, teorema daforma normal, teoria de Mathieu e método direto de Lyapunov; retratos de fases; diagramasde bifurcação; mapas de Poincaré; dimensões de Lyapunov; espectro de potência; espectrodos expoentes de Lyapunov. Já na análise não-linear experimental que utiliza séries temporais,foram vistos expoentes de Lyapunov pelo método de Wolf e o teste 0− 1.
Pretende-se nessa dissertação fazer uma modelagem mais completa, por exemplo, levandoem consideração os efeitos das diversas forças de atrito, desta forma pode-se prever a dinâmicado giroscópio real. Com isso pretende-se fazer um estudo não-linear do modelo verificandoa estabilidade no sentido de Lyapunov dos pontos de equilíbrio, outras formas de verificar aestabilidade não serão vistas, pois a ênfase é encontrar os expoentes de Lyapunov usando, paraisso, algumas abordagens clássicas, no caso: o método padrão, o método de Eckmann-Ruelle eo método de Wolf. Também será usado o teste 0− 1 por ter uma implementação rápida e fácil.
3Considera-se como zona cinza os valores de K que não estejam próximos o suficiente de zero ou de um.
14
Capítulo 3
Modelagem do Giroscópio
Neste capítulo o giroscópio utilizado neste trabalho é modelado visando obter as equaçõesdiferenciais ordinárias (EDOS) do movimento usando o método de Newton-Euler. O giroscópioconsiderado é fabricado pela Pascor, modelo ME-8960. As equações do movimento são obti-das em uma base móvel onde as equações são mais simples do que em um sistema de referênciainercial, uma vez que o tensor de inércia na base móvel é invariante com o tempo. O modeloresultante é composto por três equações diferenciais ordinárias não-lineares de segunda ordemque são transformadas para um sistema de seis equações diferenciais de primeira ordem. Assoluções deste sistema de equações compõe o espaço de estados (espaço de fases) do problemaem estudo.
A validação do modelo é feita comparando a solução numérica das equações diferenciaisordinárias, resolvidas pelo método de Runge-Kutta de oitava ordem, com resultados experimen-tais para uma dada condição inicial.
3.1 Análise Cinemática
Nesta seção são descritos os sistemas de referência e as grandezas de movimento, comoposições, velocidades e acelerações, que foram usadas para descrever o movimento do giroscó-pio. Já as causas do movimento são detalhadas na próxima seção.
3.1.1 Sistema de Referência e Matrizes de Transformação de Coordena-das
Os sistemas de referência usados são:
Sistema inercial I: Esse sistema de referência tem origem em O e seus eixos são X, Y e Z.Um vetor pode ser descrito neste sistema usando os versores ı, e k;
15
16
Sistema móvel B1: A origem deste sistema é O1 (O1 = O) contendo os eixos X1, Y1 e Z1 e osversores ı1, 1 e k1. Este sistema gira em torno do eixo Z (Z1 = Z) com um ângulo deprecessão φ e velocidade angular IΦ (figura 3.1):
IΦ =
0
0
φ
(3.1)
Sistema móvel B2: A origem deste sistema é O2 (O2 = O1) contendo os eixos X2, Y2 e Z2 eos versores ı2, 2 e k2. Este sistema gira em torno do eixo Y1 (Y2 = Y1) com um ângulode nutação θ 1 e velocidade angular B1Θ (figura 3.2):
B1Θ =
0
θ
0
(3.2)
Figura 3.1: Sistema móvel B1.
O disco gira na base B2 em torno de X2 com um ângulo de spin ψ e com velocidadeangular B2Ψ (figura 3.3):
B2Ψ =
ψ
0
0
(3.3)
1Este ângulo é limitado. Veja a seção 3.2.1 para mais detalhes.
17
Figura 3.2: Sistema móvel B2.
Figura 3.3: Ângulo de spin.
18
A mudança de sistemas de referências ocorre através das matrizes de transformação 2 (Neto,2004; Santos, 2001):
De I para B1:
T φ =
cos(φ) sin(φ) 0
− sin(φ) cos(φ) 0
0 0 1
(3.4)
B1S = T φ IS (3.5)
De B1 para I:
T Tφ =
cos(φ) − sin(φ) 0
sin(φ) cos(φ) 0
0 0 1
(3.6)
IS = T Tφ B1S (3.7)
De B1 para B2:
T θ =
cos(θ) 0 − sin(θ)
0 1 0
sin(θ) 0 cos(θ)
(3.8)
B2S = T θ B1S (3.9)
De B2 para B1:
T Tθ =
cos(θ) 0 sin(θ)
0 1 0
− sin(θ) 0 cos(θ)
(3.10)
B1S = T Tθ B2S (3.11)
sendo IS, B1S e B2S os vetores de cada sistema de referência.
Ressalta-se que a obtenção das equações diferenciais do giroscópio é mais simples em umsistema de referência móvel do que no sistema de referência inercial.
3.1.2 Vetores de Posição
Aqui serão consideradas as posições dos contrapesos e do disco na haste que os suporta,representados na base B2, assim como a posição do centro de massa, as quais são mostrados na
2As matrizes de transformação são ortogonais, isto é, têm a propriedade de que a inversa da matriz é igual asua transposta.
19
figura 3.4. O eixo da haste coincide com o eixo X2, logo as posições dos contrapesos e do discoestão sobre o eixo X2. Todos estes vetores posição têm origem no ponto de apoio O2.
Figura 3.4: Vetores de posição.
O contrapeso maior tem massa m1 e sua posição pode ser mudada para que o centro demassa do sistema haste-contrapesos-disco varie. A componente do vetor posição deste contra-peso é h1 < 0:
B2r1 =
h1
0
0
(3.12)
O contrapeso menor tem massa m2 e sua posição também pode ser mudada para que umajuste fino seja feito na posição do centro de massa. A componente do vetor posição destecontrapeso é h2 < 0:
B2r2 =
h2
0
0
(3.13)
Há também um contrapeso extra de massa m3. Neste trabalho a massa m3 não foi mode-lada, já que experimentalmente ela não foi usada.
O disco tem massa md e sua posição é fixa. A componente do vetor posição do disco é
20
hd > 0:
B2rd =
hd
0
0
(3.14)
Como a haste é assimétrica, em relação ao ponto de apoio, sua massa influencia na posiçãodo centro de massa do sistema, por isso a massa da haste mh e a posição do centro de massa damesma hh são levados em consideração.
Para a posição do centro de massa tem-se a relação (Arnold, 1987; Marion and Thornton,1995; Neto, 2004; Symon, 1996):
B2rc =
H
0
0
(3.15)
sendo:M = m1 +m2 +md +mh (3.16)
H =m1h1 +m2h2 +mdhd +mhhh
M(3.17)
3.1.3 Vetores de Velocidades e Acelerações Absolutas
A velocidade angular absoluta do sistema móvel B2 é dada pela soma das velocidadesangulares de precessão e nutação na referida base:
B2Ω = B2Φ + B2Θ (3.18)
sendo a velocidade angular de precessão descrita por:
B2Φ = T θ T φ IΦ
=
cos(θ) 0 − sin(θ)
0 1 0
sin(θ) 0 cos(θ)
cos(φ) sin(φ) 0
− sin(φ) cos(φ) 0
0 0 1
0
0
φ
=
−φ sin(θ)
0
φ cos(θ)
(3.19)
21
e a velocidade angular de nutação:
B2Θ = T θ B1Θ
=
cos(θ) 0 − sin(θ)
0 1 0
sin(θ) 0 cos(θ)
0
θ
0
=
0
θ
0
(3.20)
Substituindo as eqs. (3.19) e (3.20) na eq. (3.18) obtém-se:
B2Ω =
Ωx
Ωy
Ωz
=
−φ sin(θ)
θ
φ cos(θ)
(3.21)
A velocidade angular absoluta do disco no sistema móvel B2 é dada pela soma das velo-cidades de precessão, nutação e spin na referida base:
B2ω = B2Ω + B2Ψ (3.22)
Substituindo as eqs. (3.21) e (3.3) na eq. (3.22) obtém-se:
B2ω =
ωx
ωy
ωz
=
−φ sin(θ) + ψ
θ
φ cos(θ)
(3.23)
A aceleração angular absoluta do disco no sistema móvelB2 é descrita por (Santos, 2001):
B2ω =d
dt(B2ω) + B2Ω× B2ω (3.24)
sendo a derivada da eq. (3.23) dada por:
d
dt(B2ω) =
−φ sin(θ) + ψ − φθ cos(θ)
θ
φ cos(θ)− φθ sin(θ)
(3.25)
22
o produto vetorial da eq. (3.21) com a eq. (3.23) é:
B2Ω× B2ω =
∣∣∣∣∣∣∣ı2 2 k2
−φ sin(θ) θ φ cos(θ)
−φ sin(θ) + ψ θ φ cos(θ)
∣∣∣∣∣∣∣=
0
φψ cos(θ)
−θψ
(3.26)
Substituindo as eqs. (3.25) e (3.26) na eq. (3.24) obtém-se a aceleração angular absoluta dodisco:
B2ω =
−φ sin(θ) + ψ − φθ cos(θ)
θ + φψ cos(θ)
φ cos(θ)− φθ sin(θ)− θψ
(3.27)
A velocidade linear absoluta do centro de massa é descrita por (Santos, 2001):
B2V c = B2V O + B2Ω× B2rc + B2V rel (3.28)
sendo B2V O a velocidade da origem, ponto O2, B2V rel a velocidade relativa entre o ponto O2
e o centro de massa em B2. Como o ponto de apoio está na origem de B2 (O2) e o centro demassa, na baseB2, está fixo temos que B2V O = 0 e B2V rel = 0. Resolvendo o produto vetorialentre as eqs. (3.21) e (3.15) obtém-se:
B2Ω× B2rc =
∣∣∣∣∣∣∣ı2 2 k2
−φ sin(θ) θ φ cos(θ)
H 0 0
∣∣∣∣∣∣∣=
0
Hφ cos(θ)
−Hθ
(3.29)
Com isto encontra-se a velocidade linear absoluta do centro de massa:
B2V c =
0
Hφ cos(θ)
−Hθ
(3.30)
A aceleração linear absoluta do centro de massa na base B2 é dada por (Santos, 2001):
B2Ac = B2AO + B2Ω× B2rc + B2Ω× (B2Ω× B2rc) + 2B2Ω× B2V rel + B2Arel (3.31)
sendo B2AO a aceleração da origem e B2Arel é a aceleração relativa do centro de massa em B2.
23
Como B2V O = 0 e B2V rel = 0, tem-se que B2AO = 0, B2Ω × B2V rel = 0 e B2Arel = 0. Aderivada da eq. (3.21) em relação ao tempo é dada por:
B2Ω =
−φ sin(θ)− φθ cos(θ)
θ
φ cos(θ)− φθ sin(θ)
(3.32)
Resolvendo o produto vetorial entre as eqs. (3.32) e (3.15) encontra-se a aceleração tangencial:
B2Ω× B2rc =
∣∣∣∣∣∣∣ı2 2 k2
−φ sin(θ)− φθ cos(θ) θ φ cos(θ)− φθ sin(θ)
H 0 0
∣∣∣∣∣∣∣=
0
H[φ cos(θ)− φθ sin(θ)]
−Hθ
(3.33)
O produto vetorial entre as eqs. (3.21) e (3.29) resulta na aceleração normal:
B2Ω× (B2Ω× B2rc) =
∣∣∣∣∣∣∣ı2 2 k2
−φ sin(θ) θ φ cos(θ)
0 Hφ cos(θ) −Hθ
∣∣∣∣∣∣∣=
−H[φ2 cos2(θ) + θ2]
−Hφθ sin(θ)
−Hφ2 sin(θ) cos(θ)
(3.34)
Substituindo as eqs. (3.33) e (3.34) na eq. (3.31) encontra-se a aceleração do centro de massa:
B2Ac =
−H[φ2 cos2(θ) + θ2]
H[φ cos(θ)− 2φθ sin(θ)]
−H[θ + φ2 sin(θ) cos(θ)]
(3.35)
Com isto todas as grandezas cinemáticas necessárias para a obtenção de um modelo dadinâmica do giroscópio estão estabelecidas.
3.2 Análise Dinâmica
Nesta seção são mostradas as causas do movimento do giroscópio, isto é, as forças etorques que atuam e causam os movimentos no giroscópio. Conceitos de Mecânica Clássicaforam empregados, como a segunda lei de Newton, que relaciona força resultante e aceleraçãodo centro de massa, e também a equação de Euler, que relaciona o torque resultante, em relaçãoa origem de B2, com a variação do momento angular resultante em relação ao mesmo ponto.
24
3.2.1 Diagrama de Corpo Livre
Como já discutido anteriormente, o giroscópio da Pascor é composto por um disco commassa md, dois contrapesos (com massas m1 e m2), uma base e uma haste com massa mh naqual o disco e os contrapesos estão ligados.
No estudo do giroscópio foi feita a análise da parte do giroscópio que fica em movimento,ou seja, os contrapesos, disco e a haste na qual esses componentes estão ligados. Assume-se a hipótese de que a haste é rígida, não sofrendo deformações. Consideram-se cinco forçasatuando no giroscópio, a saber: força peso (B2F g), força de reação no ponto de apoio (B2F re),força elástica (B2F e), força de atrito cinético (B2F k) e força de atrito viscoso (B2F v).
A força peso atua em todos os pontos do corpo do giroscópio, mas é equivalente e muitomais fácil considerar a sua atuação no centro de massa. No sistema inercial o peso é escritocomo (figura 3.5):
Figura 3.5: Força peso e força de reação.
IF g =
0
0
−Mg
(3.36)
sendo g a aceleração da gravidade e M a massa do giroscópico. Como está sendo utilizadaa base móvel B2, a força peso será transformada para esta base. Primeiramente, é feita uma
25
transformação para a base móvel B1:
B1F g = T φ IF g
=
cos(φ) sin(φ) 0
− sin(φ) cos(φ) 0
0 0 1
0
0
−Mg
=
0
0
−Mg
(3.37)
e por fim, é feita a transformação da base B1 para a base B2:
B2F g = T θ B1F g
=
cos(θ) 0 − sin(θ)
0 1 0
sin(θ) 0 cos(θ)
0
0
−Mg
=
Mg sin(θ)
0
−Mg cos(θ)
(3.38)
A força de reação age no ponto de apoio, sendo responsável por manter o ponto O emequilíbrio (figura 3.5) e é descrita na base inercial como:
IF re =
Fx
Fy
Fz
(3.39)
Caso exista a necessidade de igualar essa força com outra que esteja na base B2, aplica-se atransformação:
B2F re = T θ T φ IF re (3.40)
Uma força elástica é necessária para poder limitar o ângulo de nutação 3 θ, definido nabase B2, agindo perpendicularmente à haste na posição B2re (veja figura 3.6 para maiores de-talhes) e ocorre somente em um certo intervalo de ângulo, como descrito abaixo (Neto, 2004):
θi − ε ≤ θ ≤ θf + ε (3.41)
comθi = −50π/180 rad; θf = 35π/180 rad (3.42)
3Testes experimentais preliminares mostraram a necessidade de incluir esta força para correta correlação nu-mérica e o experimental.
26
e ε a deformação elástica máxima. A força elástica tem a forma:
B2F e =
0
0
fez (θ)
(3.43)
fez (θ) =
−ke(θ − θi) : θ < θi
0 : θi ≤ θ ≤ θf
ke(θ − θf ) : θ > θf
(3.44)
onde ke a constante de elasticidade ajustada experimentalmente.
Figura 3.6: Força elástica.
Para as forças de atrito cinético que atuam nos rolamentos, pode-se verificar que o soma-tório destas é nulo pois são forças internas ao sistema, ou seja:
B2F k = 0 (3.45)
A força devida ao atrito da haste 4 e disco com o ar (viscoso) é dada por (Goldstein, 1981;Lemos, 2007; Marion and Thornton, 1995; Neto, 2004):
B2F v = −B2BB2vv (3.46)
4Haste mais os contrapesos.
27
sendo:
B2B =
bbx 0 0
0 bby 0
0 0 bbz
(3.47)
e as componentes de B2B são consideradas constantes, dependentes da geometria da haste e/oudo disco e B2vv é a velocidade linear da haste e/ou disco em um certa posição.
Assumindo que a força de atrito viscoso, que atua na superfície de um corpo, possa sersubstituída por uma força equivalente, que atua em uma linha do corpo, os pontos dessa linhasão representados por B2rv e a força atuando no corpo é obtida integrando as forças que atuamna linha.
A força de atrito viscoso que atua no disco é dada por:
B2F vd(ψ) =Fvx Fvy Fvz
T= −B2BB2ω × B2r(ψ) (3.48)
com
B2r(ψ) =
rx
ry
rz
=
hd
Rd cos(ψ)
Rd sin(ψ)
(3.49)
B2F vd(ψ) = −
∣∣∣∣∣∣∣ı2 2 k2
bbxωx bbyωy bbzωz
hd ry rz
∣∣∣∣∣∣∣=
−bbyωyrz + bbzωzry
−bbzωzhd + bbxωxrz
−bbxωxry + bbyωyhd
(3.50)
B2F vd =
∫ 2π
0B2F vd(ψ)dψ (3.51)
ωx, ωy e ωz não dependem de ψ, logo
B2F vd = −
bbyωyRd
∫ 2π
0sin(ψ)dψ − bbzωzRd
∫ 2π
0cos(ψ)dψ
bbzωzhd∫ 2π
0dψ − bbxωxRd
∫ 2π
0sin(ψ)dψ
bbxωxRd
∫ 2π
0cos(ψ)dψ − bbyωyhd
∫ 2π
0dψ
(3.52)
lembrando que ∫ 2π
0
cos(ψ)dψ =
∫ 2π
0
sin(ψ)dψ = 0 (3.53)
28
B2F vd =
0
2πbbzhdωz
2πbbyhdωy
=
0
b′vyωz
b′vzωy
(3.54)
sendo b′vy = 2πbbzhd e b′vz = 2πbbyhd constantes.
De forma análoga para a haste tem-se:
B2F vh =
0
b′′vyΩz
b′′vzΩy
(3.55)
e a força de atrito viscoso resultante é:
B2F v = B2F vd + B2F vh (3.56)
logo
B2F v =
0
bvyφ cos(θ)
bvz θ
(3.57)
3.2.2 Tensor de Inércia
O tensor de inércia calculado em relação ao ponto de apoio pode ser dividido em doistensores distintos representados na base móvel B2
5: um tensor relacionado com a haste e aoscontrapesos que estão presos à haste (B2If ), com velocidade B2Ω, e outro relacionado ao disco(B2Id), cuja velocidade é B2ω.
Os tensores de inércia são diagonais, pois os contrapesos, haste e o disco são simétricos.As grandezas necessárias para o cálculo dos tensores de inércia encontram-se na figura 3.7.
O tensor de inércia do disco em relação ao ponto de apoio é descrito por:
B2Id =
Ixx 0 0
0 Iyy 0
0 0 Izz
(3.58)
Para o disco, os momentos de inércia de massa, em relação ao centro de massa do disco, sãodados por (Santos, 2001):
Icxx =1
2mdR
2d (3.59)
5Na base móvel solidária ao movimento do corpo, o tensor de inércia é invariante com o tempo (Neto, 2004;Santos, 2001).
29
Figura 3.7: Elementos para o cálculo do tensor de inércia.
Icyy =1
4mdR
2d (3.60)
Iczz = Icyy (3.61)
As componentes em relação ao ponto de apoio são obtidas com a aplicação do teorema doseixos paralelos (Arnold, 1987; Lemos, 2007; Santos, 2001):
Ixx = Icxx (3.62)
Iyy = Icyy +mdh2d (3.63)
Izz = Icyy +mdh2d (3.64)
Logo, as componentes do tensor de inércia do disco, eq. (3.58), são dadas por:
Ixx =1
2mdR
2d (3.65)
Iyy =1
4mdR
2d +mdh
2d (3.66)
30
Izz = Iyy (3.67)
sendo Rd o raio do disco e hd a posição do centro de massa do disco em relação à origem O2.
Já o tensor de inércia das partes fixas é dado por:
B2If =
Ifxx 0 0
0 Ifyy 0
0 0 Ifzz
(3.68)
onde as componentes Ixx, Iyy e Izz são compostas pela soma das componentes dos tensores doscontrapesos e da haste:
Ifxx = I1xx + I2xx + Ihxx (3.69)
Ifyy = I1yy + I2yy + Ihyy (3.70)
Ifzz = I1zz + I2zz + Ihzz (3.71)
O contrapeso de massa m1 tem o formato de um cilindro oco, já que é preso à haste e seumomento de inércia de massa é:
I1cxx =1
2m1
(R2
1 +R2h
)(3.72)
I1cyy =1
4m1
(R2
1 +R2h
)+
1
12m1l
21 (3.73)
I1czz = I1cyy (3.74)
Usando o teorema dos eixos paralelos, tem-se:
I1xx =1
2m1
(R2
1 +R2h
)(3.75)
I1yy =1
4m1
(R2
1 +R2h
)+
1
12m1l
21 +m1h
21 (3.76)
I1zz = I1yy (3.77)
O contrapeso de massa m2 também tem formato de um cilindro oco e a haste mh é cilín-drica. Assim para m2:
I2xx =1
2m2
(R2
2 +R2h
)(3.78)
31
I2yy =1
4m2
(R2
2 +R2h
)+
1
12m2l
22 +m2h
22 (3.79)
I2zz = I2yy (3.80)
e para mh:
Ihxx =1
2mhR
2h (3.81)
Ihyy =1
4mhR
2h +
1
12mhl
2h +mhh
2h (3.82)
Ihzz = Ihyy (3.83)
sendo R1, R2 e Rh os raios; h1, h2 e hh as posições dos centros de massa das peças em relaçãoà origem O2; e l1, l2 e lh as larguras.
Substituindo as eqs. (3.75), (3.78) e (3.81) na eq. (3.69), obtêm-se Ifxx:
Ifxx =1
2
[m1R
21 +m2R
22 + (m1 +m2 +mh)R
2h
](3.84)
Substituindo as eqs. (3.76), (3.79) e (3.82) na eq. (3.70), encontram-se Ifyy:
Ifyy =1
4
[m1R
21 +m2R
22 + (m1 +m2 +mh)R
2h
]+
1
12
[m1l
21 +m2l
22 +mhl
2h
]+m1h
21 +m2h
22 +mhh
2h (3.85)
E substituindo as eqs. (3.77), (3.80) e (3.83) na eq. (3.71), encontram-se Ifzz:
Ifzz = Ifyy (3.86)
3.2.3 Torque Resultante
O torque resultante no giroscópio B2τ r, em relação ao ponto de apoio, é composto pelosseguintes torques: torque da força peso B2τ g, torque do atrito cinético B2τ k, torque do atritoviscoso B2τ v e torque da força elástica B2τ e. Ou seja:
B2τ r = B2τ g + B2τ k + B2τ v + B2τ e (3.87)
O torque causado pela força peso é dado por (Goldstein, 1981; Lemos, 2007; Neto, 2004;Takwale and Puranik, 1979):
B2τ g = B2rc × B2F g (3.88)
32
Usando as eqs. (3.15) e (3.38) na eq. (3.88), encontra-se o torque da força peso:
B2τ g =
∣∣∣∣∣∣∣ı2 2 k2
H 0 0
Mg sin(θ) 0 −Mg cos(θ)
∣∣∣∣∣∣∣ =
0
MgH cos(θ)
0
(3.89)
O torque devido à força de atrito cinético é:
B2τ k = B2τ k1 + B2τ k2 + B2τ k3 (3.90)
Este é causado por uma força de atrito no eixo de rotação que se opõe ao movimento do gi-roscópio, ou seja, o torque é constante e tem sentido contrário ao da velocidade angular doeixo. Este torque também é composto por três componentes um para cada eixo de rotação, ouseja, torque no eixo do disco (B2τ k1), torque no eixo de nutação (B2τ k2) e torque no eixo deprecessão (B2τ k3). O torque no eixo do disco, que se encontra na base B2, é descrito por:
B2τ k1 =
−Tx sgn (ψ)
0
0
(3.91)
sendo a função sgn (x) definida como:
sgn (x) =
−1 : x < 0
0 : x = 0
+1 : x > 0
(3.92)
O torque no eixo de nutação na base B1 é dado por:
B1τ k2 =
0
−Ty sgn (θ)
0
(3.93)
e na base B2 é dado por:
B2τ k2 = T θ B1τ k2
=
cos(θ) 0 − sin(θ)
0 1 0
sin(θ) 0 cos(θ)
0
−Ty sgn (θ)
0
=
0
−Ty sgn (θ)
0
(3.94)
33
Já o torque no eixo de precessão na base I é descrito como:
Iτ k3 =
0
0
−Tz sgn (φ)
(3.95)
e na base B2 fica:
B2τ k3 = T θ T φ Iτ k3 (3.96)
como:
T φ Iτ k3 =
cos(φ) sin(φ) 0
− sin(φ) cos(φ) 0
0 0 1
0
0
−Tz sgn (φ)
=
0
0
−Tz sgn (φ)
(3.97)
Substituindo a eq. (3.97) na eq. (3.96):
B2τ k3 =
cos(θ) 0 − sin(θ)
0 1 0
sin(θ) 0 cos(θ)
0
0
−Tz sgn (φ)
=
Tz sgn (φ) sin(θ)
0
−Tz sgn (φ) cos(θ)
(3.98)
Portanto, somando os torques nestes três eixos, eqs. (3.91), (3.94) e (3.98), na eq. (3.90)obtém-se:
B2τ k =
−Tx sgn (ψ) + Tz sgn(φ) sin(θ)
−Ty sgn (θ)
−Tz sgn (φ) cos(θ)
(3.99)
sendo Tx, Ty e Tz os coeficientes do torque cinético.
O torque devido ao atrito viscoso que atua no giroscópio é dividido em duas partes, umaque atua no disco (B2τ vd) e outra que atua na haste do disco (B2τ vh) que inclui os contrapesos:
B2τ v = B2τ vd + B2τ vh (3.100)
De forma geral o torque devido ao atrito viscoso é dado por:
B2τ v = B2rv × B2F v (3.101)
34
Agora, o torque do atrito viscoso que atua no disco, é descrito pelas equações:
B2τ vd(ψ) = B2r(ψ)× B2F vd(ψ) (3.102)
Note que na eq. (3.102) há dependência da posição de spin, logo:
B2τ vd(ψ) =
∣∣∣∣∣∣∣ı2 2 k2
hd ry rz
Fvx Fvy Fvz
∣∣∣∣∣∣∣ =
−ryFvz + rzFvy−rzFvx + hdFvz−hdFvy + ryFvx
(3.103)
B2τ vd(ψ) =
−bbxωx(r2y + r2z) + bbyωyhdry + bbzωzhdrz
−bbyωy(r2z + h2d) + bbzωzryrz + bbxhdωxry
−bbzωz(h2d + r2y) + bbxhdωxrz + bbyωyryrz
(3.104)
B2τ vd(ψ) =
−bbxR2
dωx + hdRd[bbyωy cos(ψ) + bbzωz sin(ψ)]
−bbyωy[R2d sin2(ψ) + h2d] + bbzR
2dωz cos(ψ) sin(ψ) + bbxhdRdωx cos(ψ)
−bbzωz[h2d +R2d cos2(ψ)] + bbxhdRdωx sin(ψ) + bbyR
2dωy cos(ψ) sin(ψ)
(3.105)
O torque viscoso resultante que atua no disco é dado por:
B2τ vd =
∫ 2π
0B2τ vd(ψ)dψ (3.106)
Lembrando que ∫ 2π
0
cos(ψ) sin(ψ)dψ = 0 (3.107)
∫ 2π
0
cos2(ψ)dψ =
∫ 2π
0
sin2(ψ)dψ = π (3.108)
Chega-se a
B2τ vd = −
2πbbxR
2dωx
πbby(R2d + 2h2d)ωy
πbbz(2h2d +R2
d)ωz
= −
b′xωx
b′yωy
b′zωz
(3.109)
onde b′x = 2πbbxR2d, b
′y = πbby(R
2d + 2h2d) e b′z = πbbz(2h
2d + R2
d). De forma análoga tem-se,para a haste; o torque:
B2τ vh = −
b′′xΩx
b′′yΩy
b′′zΩz
(3.110)
que somando ao torque do disco leva ao torque viscoso.
35
Usando as eqs. (3.109) e (3.110) na eq. (3.100) obtêm-se:
B2τ v = −
−bx2φ sin(θ) + bx1ψ
byθ
bzφ cos(θ)
(3.111)
sendo b′x + b′′x = bx2 , b′′x = bx1 , b′y + b
′′y = by e b′z + b
′′z = bz constantes com bx2 > bx1 . Por
simplificação no restante desse trabalho considera-se bx1 = bx2 = bx.
O torque devido à força elástica é:
B2τ e = B2re × B2F e =
0 τey (θ) 0T
(3.112)
sendo:
τey (θ) =
−k(θ − θi) : θ < θi
0 : θi ≤ θ ≤ θf
−k(θ − θf ) : θ > θf
(3.113)
Substituindo as eqs. (3.89), (3.99), (3.111) e (3.112) em (3.87), obtém-se o torque resul-tante:
B2τ r =
τx
τy
τz
=
−Tx sgn (ψ) +
[Tz sgn (φ) + bxφ
]sin(θ)− bxψ
MgH cos(θ)− Ty sgn (θ)− byθ + τey (θ)
−[Tz sgn (φ) + bzφ] cos(θ)
(3.114)
3.2.4 Taxa de Variação do Momento Angular Resultante
A taxa de variação do momento angular resultante, calculado em relação ao ponto deapoio O2 fixo, tem uma contribuição das partes fixas em relação à haste móvel do giroscópio(B2Lf ) e uma contribuição do disco que gira em torno desta haste (B2Ld), ou seja:
B2Lr = B2Lf + B2Ld (3.115)
A taxa de variação do momento angular para um corpo rígido na base B2 é dada por (Santos,2001):
B2L = B2Id
dt(B2ω) + B2Ω× (B2I B2ω) (3.116)
Para encontrar a variação do momento angular das partes fixas é necessário realizar asseguintes substituições na eq. (3.116): B2I = B2If e B2ω = B2Ω. Assim:
B2Lf = B2Ifd
dt(B2Ω) + B2Ω× (B2If B2Ω) (3.117)
36
Utilizando as eqs. (3.68) e (3.32) chega-se a:
B2Ifd
dt(B2Ω) =
Ifxx 0 0
0 Ifyy 0
0 0 Ifzz
−φ sin(θ)− φθ cos(θ)
θ
φ cos(θ)− φθ sin(θ)
=
−Ifxx[φ sin(θ) + φθ cos(θ)]
Ifyyθ
Ifzz[φ cos(θ)− φθ sin(θ)]
(3.118)
Similarmente, usando as eqs. (3.68) e (3.21), pode-se obter:
B2If B2Ω =
Ifxx 0 0
0 Ifyy 0
0 0 Ifzz
−φ sin(θ)
θ
φ cos(θ)
=
−Ifxxφ sin(θ)
Ifyyθ
Ifzzφ cos(θ)
(3.119)
Calculando o produto vetorial entre as eqs. (3.21) e (3.119), resulta:
B2Ω× (B2If B2Ω) =
∣∣∣∣∣∣∣ı2 2 k2
−φ sin(θ) θ φ cos(θ)
−Ifxxφ sin(θ) Ifyyθ Ifzzφ cos(θ)
∣∣∣∣∣∣∣=
−(Ifyy − Ifzz)φθ cos(θ)
(−Ifxx + Ifzz)φ2 sin(θ) cos(θ)
(Ifxx − Ifyy)φθ sin(θ)
(3.120)
Por fim, substituindo as eqs. (3.118) e (3.120) na eq. (3.117), obtém-se a taxa de variação domomento angular das partes fixas:
B2Lf =
−Ifxxφ sin(θ)− (Ifxx + Ifyy − Ifzz)φθ cos(θ)
Ifyyθ + (−Ifxx + Ifzz)φ2 sin(θ) cos(θ)
Ifzzφ cos(θ) + (Ifxx − Ifyy − Ifzz)φθ sin(θ)
(3.121)
Para encontrar a variação do momento angular do disco considera-se B2I = B2Id. Assim
B2Ld = B2Idd
dt(B2ω) + B2Ω× (B2Id B2ω) (3.122)
37
e, substituindo a eq. (3.22) na eq. (3.122),
B2Ld = [B2Idd
dt(B2Ω) + B2Ω× (B2Id B2Ω)] +
[B2Idd
dt(B2Ψ) + B2Ω× (B2Id B2Ψ)] (3.123)
A primeira parte da eq. (3.123) corresponde ao resultado anterior (3.121) mudando apenas ascomponentes do tensor de inércia. Já para a segunda parte da equação tem-se o desenvolvimentoa seguir.
Utilizando a eq. (3.58) e a derivada da eq. (3.3) obtém-se
B2Idd
dt(B2Ψ) =
Ixx 0 0
0 Iyy 0
0 0 Izz
ψ
0
0
=
Ixxψ
0
0
(3.124)
De forma análoga:
B2Id B2Ψ =
Ixxψ
0
0
(3.125)
Com o resultado da eq. (3.125):
B2Ω× (B2Id B2Ψ) =
∣∣∣∣∣∣∣ı2 2 k2
−φ sin(θ) θ φ cos(θ)
Ixxψ 0 0
∣∣∣∣∣∣∣=
0
Ixxφψ cos(θ)
−Ixxθψ
(3.126)
Por fim, substituindo as eqs. (3.124) e (3.126) na eq. (3.123), obtém-se a forma final da taxa devariação do momento angular do disco:
B2Ld =
−Ixxφ sin(θ)− (Ixx + Iyy − Izz)φθ cos(θ)
Iyyθ + (−Ixx + Izz)φ2 sin(θ) cos(θ)
Izzφ cos(θ) + (Ixx − Iyy − Izz)φθ sin(θ)
+
Ixxψ
Ixxφψ cos(θ)
−Ixxθψ
=
Ixx(−φ sin(θ) + ψ)− (Ixx + Iyy − Izz)φθ cos(θ)
Iyyθ + [(−Ixx + Izz)φ sin(θ) + Ixxψ]φ cos(θ)
Izzφ cos(θ)− [(−Ixx + Iyy + Izz)φ sin(θ) + Ixxψ]θ
(3.127)
38
Substituindo as eqs. (3.121) e (3.127) na eq. (3.115), finalmente encontra-se a taxa devariação do momento angular resultante:
B2Lr =
−Ixφ sin(θ) + Ixxψ − (Ix + Iy − Iz)φθ cos(θ)
Iyθ + [(−Ix + Iz)φ sin(θ) + Ixxψ]φ cos(θ)
Izφ cos(θ)− [(−Ix + Iy + Iz)φ sin(θ) + Ixxψ]θ
(3.128)
sendo:Ix = Ifxx + Ixx, Iy = Ifyy + Iyy, Iz = Ifzz + Izz (3.129)
3.2.5 Equação de Euler
Para obter as equações de movimento do giroscópio, resolve-se a equação de Euler nabase do sistema móvel B2 (Santos, 2001):
B2τ r = B2Lr (3.130)
Utilizando as eqs. (3.114) e (3.128) obtém-se um sistema de equações diferenciais:τx
τy
τz
=
−Ixφ sin(θ) + Ixxψ − (Ix + Iy − Iz)φθ cos(θ)
Iyθ + [(−Ix + Iz)φ sin(θ) + Ixxψ]φ cos(θ)
Izφ cos(θ)− [(−Ix + Iy + Iz)φ sin(θ) + Ixxψ]θ
(3.131)
Resolvendo o sistema de eqs. (3.131) nas variáveis φ, θ e ψ obtém-se:
φ =τz + [(−Ix + Iy + Iz)φ sin(θ) + Ixxψ]θ
Iz cos(θ)(3.132)
θ =τy − Ixxφψ cos(θ) + (Ix − Iz)φ2 cos(θ) sin(θ)
Iy(3.133)
ψ =τx + (Ix + Iy − Iz)φθ cos(θ) + Ix
Izτz + [(−Ix + Iy + Iz)φ sin(θ) + Ixxψ]θ tan(θ)
Ixx(3.134)
As equações diferenciais (3.132) e (3.134) foram obtidas levando em consideração que cos(θ) 6=0 e pela expressão (3.41) tal condição é satisfeita.
Usando o conceito de variáveis de estado, pode-se converter o sistema de três equaçõesdiferenciais de segunda ordem para um sistema de seis equações de primeira ordem:
X = F (X); X ∈ <6 (3.135)
39
sendoX o vetor de estados definido por:
X =x1 x2 x3 x4 x5 x6
T=φ φ θ θ ψ ψ
T(3.136)
Portanto, as equações diferenciais no espaço de estado são dadas por:
x1 = x2 (3.137)
x2 =τz + [(−Ix + Iy + Iz)x2 sin(x3) + Ixxx6]x4
Iz cos(x3)(3.138)
x3 = x4 (3.139)
x4 =τy − Ixxx2x6 cos(x3) + (Ix − Iz)x22 cos(x3) sin(x3)
Iy(3.140)
x5 = x6 (3.141)
x6 = τx+(Ix+Iy−Iz)x2x4 cos(x3)Ixx
+ Ixτz+[(−Ix+Iy+Iz)x2 sin(x3)+Ixxx6]x4
IxxIztan(x3) (3.142)
onde τx, τy e τz são obtidos das eqs. (3.113) e (3.114), utilizando a definição dos estados (3.136)na forma
B2τ r =
τx
τy
τz
=
−Tx sgn (x6) + [Tz sgn (x2) + bxx2] sin(x3)− bxx6MgH cos(x3)− Ty sgn (x4)− byx4 + τey (x3)
−[Tz sgn (x2) + bzx2] cos(x3)
(3.143)
com
τey (x3) =
−k(x3 − θi) : x3 < θi
0 : θi ≤ x3 ≤ θf
−k(x3 − θf ) : x3 > θf
(3.144)
O sistema no espaço de estados dado pelas eqs. (3.137) até (3.142) pode ser solucionadonumericamente por meio da aplicação de uma condição inicial X0 e da implementação dealgum método numérico como Runge-Kutta de oitava ordem de passo variável, usado nessetrabalho.
40
3.2.6 Equação de Newton
Para obter as equações das reações dinâmicas no ponto de apoio, aplica-se a segunda leide Newton (Marion and Thornton, 1995; Santos, 2001):
B2F r = MB2Ac (3.145)
sendo a força resultante dada por:
B2F r = B2F g + B2F re + B2F e + B2F k + B2F v (3.146)
Substituindo a eq. (3.145) na eq. (3.146) e isolando B2F re, obtém-se:
B2F re = MB2Ac − [B2F g + B2F e + B2F k + B2F v] (3.147)
Com as eqs. (3.35), (3.38), (3.43), (3.45) e (3.57) aplicadas na eq. (3.147) obtêm-se:
B2F re =
−MH[φ2 cos2(θ) + θ2] + g sin(θ)
MH[φ cos(θ)− 2φθ sin(θ)]− bvzφ cos(θ)
−MH[θ + φ2 sin(θ) cos(θ)]− g cos(θ)+ fez (θ) + bvyθ
(3.148)
Para expressar as forças de reação na base inercial 6, usam-se as eqs. (3.39), (3.40) e isola-se
IF re: Fx
Fy
Fz
= T Tφ T
Tθ B2F re (3.149)
O produto das matrizes de transformação na eq. (3.149) é dado por:
T Tφ T
Tθ =
cos(φ) − sin(φ) 0
sin(φ) cos(φ) 0
0 0 1
cos(θ) 0 sin(θ)
0 1 0
− sin(θ) 0 cos(θ)
=
cos(φ) cos(θ) − sin(φ) cos(φ) sin(θ)
sin(φ) cos(θ) cos(φ) sin(φ) sin(θ)
− sin(θ) 0 cos(θ)
(3.150)
Com as eqs. (3.150) e (3.148) aplicadas na eq. (3.149), obtêm-se as forças de reação na
6Assumindo que os sensores de força sejam instalados na base inercial.
41
base inercial:
Fx = −M [Hθ +Hφ2 cos(θ) sin(θ)− g cos(θ)] + bvyθ + fez (θ) cos(φ) sin(θ)−−M [Hφ2 cos2(θ) +Hθ2 + g sin(θ)] cos(φ) cos(θ) +
+[(−MHφ+ bvzφ) cos(θ) + 2MHφθ sin(θ)] sin(φ) (3.151)
Fy = −M [Hθ +Hφ2 cos(θ) sin(θ)− g cos(θ)] + bvyθ + fez (θ) sin(φ) sin(θ)−−M [Hφ2 cos2(θ) +Hθ2 + g sin(θ)] sin(φ) cos(θ) +
+[(MHφ− bvzφ) cos(θ)− 2MHφθ sin(θ)] cos(φ) (3.152)
Fz = −M [Hθ +Hφ2 cos(θ) sin(θ)− g cos(θ)] + bvyθ + fez (θ) cos(θ) +
+M [Hφ2 cos2(θ) +Hθ2 + g sin(θ)] sin(θ) (3.153)
Expressando em termos das variáveis de estado, eq. (3.136), obtém-se:
Fx = −M [Hx4 +Hx22 cos(x3) sin(x3)− g cos(x3)] + bvyx4 + fez (x3) cos(x1) sin(x3)−−M [Hx22 cos2(x3) +Hx24 + g sin(x3)] cos(x1) cos(x3) +
+[(−MHx2 + bvzx2) cos(x3) + 2MHx2x4 sin(x3)] sin(x1) (3.154)
Fy = −M [Hx4 +Hx22 cos(x3) sin(x3)− g cos(x3)] + bvyx4 + fez (x3) sin(x1) sin(x3)−−M [Hx22 cos2(x3) +Hx24 + g sin(x3)] sin(x1) cos(x3) +
+[(MHx2 − bvzx2) cos(x3)− 2MHx2x4 sin(x3)] cos(x1) (3.155)
Fz = −M [Hx4 +Hx22 cos(x3) sin(x3)− g cos(x3)] + bvyx4 + fez (x3) cos(x3) +
+M [Hx22 cos2(x3) +Hx24 + g sin(x3)] sin(x3) (3.156)
3.2.7 Equação da Trajetória na Base Inercial
Para saber a trajetória de um ponto da haste do giroscópio no sistema inercial, bastatransformar um vetor na base móvel B2 para a base inercial I via matrizes de transformação.
42
Para a posição do centro de massa do disco temos:
B1rd = T Tθ B2rd
=
cos(θ) 0 sin(θ)
0 1 0
− sin(θ) 0 cos(θ)
hd
0
0
=
hd cos(θ)
0
−hd sin(θ)
(3.157)
Usando a expressão (3.157) pode-se encontrar a trajetória do centro de massa do disco:
Ird = T Tφ B1rd
=
cos(φ) − sin(φ) 0
sin(φ) cos(φ) 0
0 0 1
hd cos(θ)
0
−hd sin(θ)
=
hd cos(φ) cos(θ)
hd sin(φ) cos(θ)
−hd sin(θ)
(3.158)
3.3 Simulação e Comparação Experimental
A simulação numérica do giroscópio é feita encontrando-se uma solução usando o métodode Runge-Kutta de oitava ordem de passo variável para o sistema de equações diferenciais parauma dada condição inicial e um certo intervalo de tempo.
Antes de fazer a simulação é necessário obter os parâmetros da equação diferencial. Noapêndice A é mostrado o procedimento usado para estimar experimentalmente os parâmetrosobtidos por medidas indiretas, como os relacionados ao torque viscoso (bx, by, bz), torque ci-nético (Tx, Ty, Tz), torque elástico k e a aceleração da gravidade g. Os parâmetros obtidos pormedidas diretas são mostrados na tabela 3.1.
Tabela 3.1: Parâmetros dos contrapesos, da haste e do disco do giroscópio da Pascor. Sendom as massas, h as posições, l as larguras e R os raios. As massas m1 e m2 foram colocadas emuma posição tal que o giroscópio está em equilíbrio.
m( kg) h( m) l( m) R( m)
md 1,740 0,1465 0,127m1 0,9 −0,2579 0,0315 0,035m2 0,03 −0,1720 0,019 0,0225mh 0,2714 −0,065 0,574 0,0063
43
O giroscópio possuí dois sensores de posição angular que fazem as leituras dos movi-mentos de precessão e nutação 7 (figura 3.8). A aquisição de dados é feita com o programaDataStudio (figura 3.9(a)) com o auxílio da interface de aquisição de dados da Pascor, mo-delo CI-6859C 8 (figura 3.9(b)). No programa de aquisição de dados foi usada uma taxa deamostragem de 1000 Hz.
(a) Sensor de precessão. (b) Sensor de nutação.
Figura 3.8: Sensores de posição angular da Pascor, modelo CI-6625 (RMS).
Como não há sensores para velocidades angulares, uma alternativa para obter tais veloci-dades no tempo seria por um processo de derivação numérica dos dados medidos de desloca-mentos angulares, porem a derivação numérica se mostrou ruidosa e pouco precisa e, por isto,foi omitido desta dissertação.
O experimento realizado foi feito com o giroscópio desequilibrado, estando o contrapesomenor encostado próximo ao ponto de apoio na posição h2 = −0,0827 m, e o contrapeso maiorna posição h1 = −0,1779 m. As distâncias entre as faces mais próximas dos contrapesos é de7,0 cm.
Antes de iniciar o experimento a velocidade de spin é estimada, primeiramente usa-se umtacômetro (figura 3.9(c)) que mede a velocidade do disco através de uma marca reflexiva feitano mesmo para que a velocidade possa ser lida, não pode haver movimento de precessão nem denutação. O disco é posto a girar através de um cordão que é enrolado na base do disco e atravésde um sistema de roldanas o movimento de queda de uma esfera é transferido para o disco.O valor encontrado para a velocidade de spin foi 68,49 (rad/s), posteriormente no apêndice Aobteve-se um valor melhorado.
7Veja a referência (PASCO, 1996) para mais detalhes.8Veja a referência (PASCO, 1997) para mais detalhes.
44
(a) Programa DataStudio. (b) Interface de aquisição de dados.
(c) Tacômetro.
Figura 3.9: Aquisição de dados.
Na realização do experimento admite-se que a velocidade de spin será a mesma da que foiestimada anteriormente já que a condição adotada é a mesma. Para o experimento mantém-se ogiroscópio parado manualmente, com a posição de nutação x3 = 0,0 rad. Isso é feito com umtransferidor preso ao giroscópio. Simultaneamente o disco é posto a girar e, logo em seguida,o giroscópio é solto e inicia-se a aquisição de dados. O tempo de duração dos experimentosnão pode ser muito longo, pois o cabo do sensor da posição angular de nutação induz umamortecimento e prende o movimento depois de poucas voltas, dificultando uma aquisição dedados mais longa e confiável (figura 3.10). Para o experimento em questão o tempo total domovimento foi de tf = 27,54 s.
As condições iniciais foram extraídas dos dados experimentais, no caso as velocidadesde precessão e de nutação foram obtidas por uma derivação numérica simples. Como pode-se ver na tabela 3.2 apesar de se querer que o giroscópio ficasse parado, não foi o que ocorreupois segurando o giroscópio manualmente acaba-se introduzindo alguma velocidade no sistema,sendo assim não se tem controle da condição inicialX0.
Na simulação são usadas as mesmas condições iniciais e o mesmo intervalo de tempoque foi usado no experimento para poder compará-los. No método de Runge-Kutta foi adotadoum erro de 1,0 × 10−12 tanto para o erro absoluto quanto para o erro relativo. A comparação
45
Figura 3.10: Fio do sensor de nutação.
Tabela 3.2: Condições iniciais para validação do modelo do giroscópio da Pascor. Com h1 =−0,1779 m e h2 = −0,0827 m.
x1 (rad) x2 (rad/s) x3 (rad) x4 (rad/s) x5 (rad) x6 (rad/s)
X0 0,1920 2,0854 0,2090 0,7708 0,0000 60,777
das posições obtidas da simulação (linha azul) e do experimento (linha vermelha) encontra-se nafigura 3.11. As velocidades simuladas são mostradas na figura 3.12. A comparação da trajetóriasimulada e experimental do disco pode ser visto na figura 3.13.
0 5 10 15 20 25 30−5
0
5
10
15
20
25
30
35
tempo [s]
x 1 [rad
]
SimuladoExperimental
(a) Precessão.
0 5 10 15 20 25 30−0.1
0
0.1
0.2
0.3
0.4
0.5
0.6
tempo [s]
x 3 [rad
]
SimuladoExperimental
(b) Nutação.
0 5 10 15 20 25 300
200
400
600
800
1000
1200
1400
1600
tempo [s]
x 5 [ra
d]
(c) Spin simulado.
Figura 3.11: Posições angulares de precessão, nutação e spin.
Pode-se verificar dos gráficos que os resultados experimentais e simulados têm tendênciasparecidas. Para o caso do movimento de spin, o giroscópio não tem um sensor específico paramedição e, portanto não se apresenta a comparação com dados experimentais. O mesmo valepara os sinais de velocidade angular que não são medidos em nenhum dos movimentos.
46
0 5 10 15 20 25 30−1
−0.5
0
0.5
1
1.5
2
2.5
tempo [s]
x 2 [rad
/s]
(a) Precessão.
0 5 10 15 20 25 30−1.5
−1
−0.5
0
0.5
1
1.5
tempo [s]
x 4 [rad
/s]
(b) Nutação.
0 5 10 15 20 25 3051
52
53
54
55
56
57
58
59
60
61
tempo [s]
x 6 [ra
d/s
]
(c) Spin simulado.
Figura 3.12: Velocidades angulares de precessão, nutação e spin.
−0.2
−0.1
0
0.1
0.2
−0.2
−0.1
0
0.1
0.2
−0.04
−0.03
−0.02
−0.01
0
0.01
X [m]Y [m]
Z [m
]
(a) Simulado.
−0.2
−0.1
0
0.1
0.2
−0.2
−0.1
0
0.1
0.2
−0.04
−0.03
−0.02
−0.01
0
0.01
X [m]Y [m]
Z [m
]
(b) Experimental.
Figura 3.13: Trajetória do centro de massa do disco.
3.4 Considerações Finais
Este capítulo detalhou a modelagem do giroscópio da Pascor utilizando conceitos damecânica de Newton-Euler. Nesta modelagem, foi levado em consideração efeitos elásticosque limitam o movimento de nutação, assim como efeitos dissipativos causados pelos atritosdas partes girantes (rolamentos) e do movimento do giroscópio com o ar. O resultado destemodelo foi um conjunto de equações diferenciais não-lineares que descreve o movimento deprecessão, nutação e spin do giroscópio.
As equações diferenciais foram resolvidas numericamente com o método de Runge-Kuttade oitava ordem com passo variável e as soluções foram comparadas com os dados experimen-tais que puderam ser medidos, como deslocamento angular de precessão e nutação. Algunspontos são destacados a seguir:
• A condição inicial experimental só pode ser verificada numericamente com o modelo deforma aproximada, por não ser possível garanti-la com precisão;
47
• Não se teve acesso com o sistema experimental aos dados de velocidades angulares, queportanto, não foram comparadas com as velocidades da simulação. Para a velocidadeinicial de spin foi possível estima-la com o uso de um tacômetro, mas não foi possívelsaber tal velocidade no tempo pois para se obter esta informação, é necessário que o discodo giroscópio esteja parado em frente ao tacômetro;
• O procedimento de obtenção dos parâmetros indiretos como os relacionados aos torquesdissipativos e ao torque elástico foram feitos experimentalmente e, naturalmente estãosujeitos a um certo grau de incerteza.
Apesar disto, os resultados mostraram que o modelo proposto apresenta uma boa concor-dância com o movimento real do giroscópio e pode ser utilizado para as finalidades de análisenão-linear que são mostradas nos próximos capítulos desta dissertação.
48
Capítulo 4
Análise Dinâmica Não-Linear
Apresenta-se neste capítulo um resumo dos conceitos de sistemas dinâmicos não-linearesque são usados no capítulo 5 deste trabalho. Os conceitos abordados são pontos de equilíbrio,e estabilidade no sentido de Lyapunov. Não serão vistos outros tipos de estabilidade, pois aênfase do trabalho são os expoentes de Lyapunov. Para completar é apresentado o teste 0 − 1
que tem uma implementação fácil e que permite verificar se o sistema é regular ou apresentacaós ou caós fraco. Os métodos para o cálculo dos expoentes de Lyapunov aqui abordados são:o método padrão, o método de Eckmann-Ruelle e o método de Wolf.
Para muitos dos conceitos abordados a seguir considera-se um sistema dinâmico autô-nomo formado pelo sistema de equações diferenciais de ordem n visto abaixo:
X = F (X) (4.1)
sendoX =
x1 x2 . . . xn
T(4.2)
eF (X) =
f1(X) f2(X) . . . fn(X)
T(4.3)
4.1 Pontos de Equilíbrio e Estabilidade no Sentido de Lyapu-nov
Os pontos de equilíbrio Xe correspondem às soluções da eq. (4.1) que não mudam como passar do tempo, ou seja:
F (Xe) = 0 (4.4)
A estabilidade dos pontos de equilíbrio pode ser classificada como (Fiedler-Ferrara anddo Prado, 1994; Wiggins, 2003; Monteiro, 2006):
49
50
Neutralmente ou marginalmente estável - dado ε > 0, ∃δ = δ(ε) > 0 tal que se |X0 −Xe| < δ então |X(t) − Xe| < ε, ∀t > 0. Se o δ existir, o ponto de equilíbrio serámarginalmente estável;
Instável - usando a definição de marginalmente estável, caso não exista um δ que satisfaça osrequisitos da definição, o ponto de equilíbrio será instável;
Assintoticamente estável - Para que um ponto de equilíbrio seja assintoticamente estável énecessário que ele seja marginalmente estável e que também existam condições ini-ciais que satisfazem a relação |X0 − Xe| < δ, para algum δ > 0 e que a relaçãolimt→∞X(t) = Xe seja satisfeita. Tal ponto é denominado atrator.
Para muitos sistemas não-lineares, a estabilidade dos pontos de equilíbrio pode ser verifi-cada através da análise da linearização do sistema, ou seja verifica-se os sinais dos autovaloresobtidos da matriz Jacobiana do sistema linearizado. A estabilidade do sistema linearizado im-plica na estabilidade do sistema não-linear desde que o ponto de equilíbrio seja hiperbólico, ouseja não deve existir autovalor com parte real nula no sistema linearizado. Caso o autovalor te-nha parte real nula o ponto de equilíbrio é dito elíptico. No caso de pontos de equilíbrio elípticosa análise do sistema linear não é conclusiva e para saber a estabilidade do ponto de equilíbriodo sistema não-linear é necessário a utilização de outros meios como o teorema da variedadecentral e o método direto de Lyapunov (Fiedler-Ferrara and do Prado, 1994; Monteiro, 2006).
Caso a parte real de todos os autovalores relacionados a um ponto de equilíbrio sejamnegativas, o ponto de equilíbrio é estável. Caso exista ao menos um autovalor, relacionado aum ponto de equilíbrio, com parte real positiva, o ponto de equilíbrio é instável.
4.2 Expoentes de Lyapunov
Expoentes de Lyapunov são úteis para verificar se um sistema dinâmico não-linear, for-mado por n equações diferenciais ordinárias de primeira ordem, apresenta sensibilidade às con-dições iniciais, ou seja, identificar se o sistema tem comportamento caótico. A determinaçãoda sensibilidade às condições iniciais é feita considerando uma hiperesfera centrada em umacondição inicial X(t0). Na superfície da hiperesfera há n condições iniciais, uma em cada di-reção principal que forma, juntamente com X(t0), os raios iniciais ri(t0). Conforme o tempopassa esses n raios ri(t) evoluem no tempo deformando a hiperesfera. Assumindo que essaevolução seja exponencial 1 tem-se a seguinte relação que dá uma ideia inicial sobre expoentesde Lyapunov (Wolf et al., 1985):
ri(t) = ri(t0)eλi(t−t0); i = 1, 2, · · · , n (4.5)
1Pode-se usar outras bases além da base definida pela exponencial e.
51
Isolando os expoentes λi e considerando um tempo de evolução longo:
λi = limt→∞
ln(ri(t)ri(t0)
)t− t0
(4.6)
Os expoentes λi são conhecidos como expoentes de Lyapunov e representam a médiadas expansões ou contrações das soluções próximas no espaço de fases. Estes expoentes têmalgumas propriedades (Wolf et al., 1985; Monteiro, 2006). Para um sistema conservativo, porexemplo:
n∑i=1
λi = 0, (4.7)
já para um sistema dissipativo:n∑i=1
λi < 0. (4.8)
Em relação às propriedades individuais dos expoentes λi, pode-se destacar que (Wolfet al., 1985; Monteiro, 2006):
• Se λi = 0, na média as distâncias entre duas condições iniciais se mantém, sendo umasolução periódica;
• Para λi < 0, tem-se que o raio ri(t) diminui, ou seja, as trajetórias no espaço de fases seaproximam (contração de volume no espaço de fases);
• No caso de λi > 0, há divergência exponencial das trajetórias e, portanto, sensibilidadeàs condições iniciais, apresentando comportamento caótico.
Para o sistema ser caótico basta que exista pelo menos um λi > 0. Caso nenhum dosexpoentes seja positivo, o comportamento do sistema é regular. Com uma divergência exponen-cial, um erro nas condições iniciais faz com que a trajetória se afaste rapidamente da trajetóriaque o sistema teria caso não houvesse este erro. Em um sistema dissipativo, quando existecaos, o atrator no espaço de fases é dito estranho. Tal atrator tem um volume limitado, apesarda divergência exponencial causada pelo expoente positivo. O volume finito do atrator podeser explicado pelos repetidos processos de alongamentos e dobras (Wolf et al., 1985; Mon-teiro, 2006; Savi, 2006). Os alongamentos ocorrem nas direções em que os expoentes de Lya-punov são positivos. Nas direções em que os expoentes são negativos há contração do volumedo atrator. As dobras ocorrem modificando as direções das expansões e contrações.
4.2.1 Expoentes de Lyapunov do Modelo
Nesta seção são descritos o método padrão e o método de Eckmann-Ruelle que são uti-lizados para encontrar os expoentes de Lyapunov de um sistema não-linear autônomo com n
52
equações diferenciais, eq. (4.1).
Método Padrão
Para o cálculo dos expoentes de Lyapunov com o método padrão é usada uma condiçãoinicialX(t0) e a partir disso o sistema evolui no tempo, conforme a solução do sistema de equa-ções diferenciais não-lineares. A solução desse sistema de equações diferenciais origina a tra-jetória fiducial ou trajetória de referência. No mesmo instante, n condições iniciais, orientadasnas n direções canônicas em torno da condição inicial da trajetória fiducial, formam raios ri(t0)unitários. A evolução no tempo do sistema partindo das n condições iniciais ocorre conforme asolução da equação diferencial linearizada η em torno da trajetória fiducial X = X(t) (Wolfet al., 1985):
η = J(X)η, J(X) = [aij]n×n, fi = fi(X), aij =∂fi∂xj
, (4.9)
sendo J(X) a matriz jacobiana obtida da equação diferencial não-linear original. Com isso épossível observar se as trajetórias linearizadas se afastam ou se aproximam da trajetória fiducial.A resolução da equação diferencial não-linear é feita simultaneamente com as n equações line-arizadas, ou seja, é resolvida a equação diferencial expandida, com condição inicial expandida.Para n direções ortogonais (condições iniciais) tem-se:[
η1 η2 · · · ηn
]= J(X)
[η1 η2 · · · ηn
](4.10)
ηi =η1i η2i . . . ηni
T(4.11)
sendo a equação diferencial expandida da forma:
XE = F E(XE) (4.12)
onde:XE =
X η1 . . . ηn
T(4.13)
e
F E(XE) =
F (X)
J(X)η1...
J(X)ηn
(4.14)
ondeXE e F E são vetores de dimensão n(n+ 1).
Com as soluções ηi obtêm-se aproximações para os expoentes de Lyapunov (log2 |ηi|/t).Em determinados instantes de tempo é necessário fazer uma ortonormalização, pois, os raios
53
ri(t) podem crescer rapidamente, além da base formada pela evolução temporal das condiçõesiniciais deixarem de ser ortogonal. Para fazer a ortonormalização pode-se usar o processo deGram-Schmidt (Wolf et al., 1985; Ramasubramanian and Sriram, 2000; Fiedler-Ferrara anddo Prado, 1994; Monteiro, 2006; Savi, 2006), descrito na sequência deste texto.
Para uma base[η1 · · · ηn
]uma nova base pode ser construída
[η′1 · · · η′j · · · η′n
]η′1 =
η1
|η1|(4.15)
η′j = ηj −j−1∑i=1
(ηj · η′i)η′i (4.16)
η′j =η′j|η′j|
(4.17)
Pode-se usar um tempo (passo) entre ortonormalizações fixo ∆t e, na sequência, verificara aproximação dos expoentes de Lyapunov. O valor |ηi| é armazenado antes do processo deortonormalização:
λik =
∑kj=1 log2 |ηi(j∆t)|
k∆t. (4.18)
λi = limk→∞
λik (4.19)
sendo λi os expoentes de Lyapunov.
Pode-se adotar como critério de convergência um valor alto de k (consequentementet). Tomando os últimos m λik e calculando a média, ou seja, λi = (λik + λi(k−1) + · · · +λi(k−m+1))/m e o erro é o desvio padrão da média σi = σi/
√m. Em vez de usar um k alto,
pode-se empregar um critério de parada baseado em σi.
Método de Eckmann-Ruelle
No método de Eckmann-Ruelle são feitas algumas alterações na equação diferencial ori-ginal do sistema de forma a poder usar a definição de expoentes de Lyapunov para mapasmultidimensionais. Dado um mapa de dimensão m tem-se:
xn+1 = F (xn). (4.20)
54
Os expoentes de Lyapunov para o mapa (Eckmann and Ruelle, 1985; Fiedler-Ferrara and do Prado,1994) são definidos como 2
λj(x0) = limN→∞
1
Nlog2 |Λj,N |, j = 1, 2, . . . ,m (4.21)
sendo |Λj,N | o módulo da j-ésima componente dos autovalores da matrizM
M =N−1∏i=0
Γ(xi), Γ(xi) =∂F
∂x
∣∣∣∣x=xi
, (4.22)
no caso Γ(xi) é a matriz Jacobiana no ponto xi.
Com a definição dos expoentes de Lyapunov dada pela eq. (4.21), pode-se agora alterar aequação diferencial continua no tempo (4.1) para poder utilizar esta definição. Deve-se lembrarque os expoentes de Lyapunov são invariantes em relação às transformações na EDO que serãofeitas.
Assume-se a EDO representativa do sistema dinâmico de interesse:
dx
dt= F (x(t)) = U (t), (4.23)
logo:dU
dt=
dF
dx
dx
dt= J(x)U (t) (4.24)
Assumindo que a solução de (4.23) seja x(t) = f(x0, t), assim:
U(t) =dx
dt=
df(x0, t)
dx0
dx0
dt= j0(t)U 0 (4.25)
edU
dt= U 0
dj0(t)
dt(4.26)
igualando as eqs. (4.24) e (4.26) e usando a eq. (4.25) tem-se:
dj0(t)
dt= J(x)j0(t) (4.27)
Ao obter a solução de j0(t) numericamente, usando um “passo” de integração τ , as so-luções são como um mapa com ti = iτ . Considerando que j0(ti) = Γ(xi) e utilizando aseqs. (4.21) e (4.22), encontra-se os expoentes de Lyapunov pelo método de Eckmann-Ruelle.
2Será usado log2 em vez de ln, por questão de compatibilidade com o método padrão desenvolvido anterior-mente.
55
4.2.2 Expoentes de Lyapunov Experimental - Método de Wolf
O método de Wolf para séries temporais é utilizado para se obter os expoentes de Lya-punov de um sistema experimental. Com o método de Wolf é possível obter os dois maioresexpoentes de Lyapunov, porém, este trabalho irá tratar somente do algoritmo para obtenção domaior expoente com tempo de evolução fixo, ou seja, o cálculo das aproximações dos expoentesde Lyapunov é feito com variação de tempo fixo.
O método de Wolf é baseado nas ideias do método padrão que é utilizado quando seconhece a equação diferencial do sistema. No método padrão tem-se a trajetória fiducial, que éobtida pela solução numérica da equação diferencial não-linear. No método de Wolf a trajetóriafiducial é o espaço de estados reconstruído através da série temporal experimental.
No método padrão as n soluções linearizadas, que são obtidas da equação diferencial ori-ginal, são comparadas com a trajetória fiducial para se obter as aproximações sucessivas dosexpoentes de Lyapunov. Para cada aproximação sucessiva é realizado o processo de ortonor-malização de Gram-Schmidt e no caso o primeiro vetor vai tendendo à uma direção que tendea maximizar o expoente de Lyapunov correspondente, enquanto o primeiro vetor não sofre mo-dificações além de ficar unitário após o processo (Wolf et al., 1985).
Para o método de Wolf, a comparação com a trajetória fiducial, para posterior cálculo doexpoente de Lyapunov, é feita com uma trajetória próxima desta trajetória fiducial. Esta traje-tória deve ser um trecho próximo o suficiente e faz parte da trajetória original, mas se compor-tando como se fosse uma outra solução com condição inicial próxima. Com estas aproximaçõessucessivas, é feito um procedimento que tem o mesmo objetivo do processo de Gram-Schmidt,que é usado no método padrão. O procedimento de Wolf consiste em obter um ponto o maispróximo da trajetória fiducial vigente e que seja mantida a orientação do vetor formado pelascondições iniciais, o mais próximo possível do vetor anterior (veja o ângulo θ1 na figura 4.1).Para se empregar o método de Wolf há duas etapas: a primeira consiste em obter o atratorreconstruído, e a segunda usar o método propriamente dito.
Figura 4.1: Exemplificação do método de Wolf. Adaptado de (Wolf et al., 1985)
56
Em experimentos, geralmente não se tem ou mesmo se conhece todas as variáveis de es-tado do sistema em estudo. Normalmente tem-se uma série temporal de dados experimentais deum dos estados do sistema. Para fazer um estudo do sistema é preciso ter o atrator do mesmo.Tal atrator pode ser reconstruído a partir da série temporal através do método dos atrasos tem-porais de Takens (Monteiro, 2006). O atrator reconstruído é topologicamente equivalente aoatrator original, tendo os mesmos invariantes, como por exemplo os expoentes de Lyapunov.
Dada a série temporal xj = x[j] = x(tj) = x(jδt), com j = 0, . . . , N − 1, sendoN o número de pontos e δt a taxa de amostragem da série temporal. No método dos atrasostemporais de Takens os pontos do atrator reconstruído tem a forma:
Si =xi xi+τ . . . xi+(m−1)τ
T, (4.28)
sendo i a evolução temporal ou tempo discreto, m a dimensão de imersão e τ o tempo de atraso.O número de pontos do atrator reconstruído é N ′
= N − (m − 1)τ . O tempo de atraso éobtido pelo primeiro mínimo da informação mútua (Fraser and Swinney, 1986), e a dimensãode imersão é obtida pelo método falsos vizinhos próximos (Perc, 2006; Kodba et al., 2005).
Com o atrator reconstruído, pode-se calcular o maior expoente de Lyapunov do sistemaem estudo usando a seguinte expressão:
λm =1
te
ν−1∑j=0
log2
(L
′j
Lj
)(4.29)
as variáveis Lj , L′j , ν e te são descritas a seguir.
No caso ν é o número total de passos de substituição:
ν = (N′ − 1)/∆t (4.30)
e como ν é discreto a divisão que aparece na eq. (4.30) é inteira, e assim, o tempo total deevolução te é:
te = νδt∆t (4.31)
sendo ∆t o número de evoluções (δt∆t é o tempo de uma evolução). Idealmente ∆t > τ e ∆t
não deve ser muito maior que mτ , isso é feito para que a análise da evolução do sistema nãopasse por uma região de dobragem no atrator reconstruído (Kodba et al., 2005; Perc, 2006).
Inicialmente tem-se uma trajetória próxima da trajetória fiducial, cuja separação Lj é:
Lj = |S(j∆t)− S(κj)| (4.32)
e após um tempo de evolução do sistema, as trajetórias ficam distantes entre si deL′j (figura 4.1):
57
L′
j = |S((j + 1)∆t)− S(κj + ∆t)| (4.33)
Caso as trajetórias se afastem L′j > Lj e se elas se aproximem L
′j < Lj .
Para se encontrar o valor do maior expoente de Lyapunov é necessário achar os valores deκj . O valor de κj representa o tempo da posição da trajetória, que é representado pela funçãovetorial S(.). Encontrando κj , encontra-se as trajetórias pelas quais se deve comparar com atrajetória fiducial.
Para κ0, procura-se um ponto da trajetória que comparado com S(0) tenha a menor dis-tância além de κ0 ≥ 10. Com esta última imposição evita-se uma separação temporal pequenaque contribuiria para encontrar um expoente de Lyapunov zero ao invés do maior expoente.Com isso, κ0 deve minimizar L0 da eq. (4.32), com a ressalva de que L0 < υmax, sendo υmax(υmim) a máxima (mínima 3) distância permitida entre as “duas” trajetórias.
Agora para κj , com j = 1, . . ., busca-se |κj − j| ≥ 10 e procura-se um κj que minimizeLj com υmim < Lj < pυmax, sendo p uma ponderação e que em condições normais vale p = 1;é necessário também que obedeça a seguinte relação:
arccos
(Lj ·Lj−1LjLj−1
)< θmin = 0,3 rad (4.34)
Esta relação é equivalente ao processo de ortonormalização de Gram-Schmidt que é usadono método padrão, mantendo o vetor Lj com quase a mesma orientação do vetor anterior Lj−1com uma diferença máxima de 0,3 rad.
Caso não seja encontrado um κj que satisfaça os critérios acima, muda-se o valor da pon-deração (p = 2) e repete-se os passos para obter κj e se mesmo assim o κj não for encontrado ovalor de p vai sendo alterado de unidade em unidade até um o valor máximo de p = 5. Depoisdisso se κj ainda não for encontrado, p volta a ter o valor p = 1 e θmin tem o seu valor dobradoe o processo se repete até que seja encontrado o respectivo κj ou até que o θmin > 3,14 rad
situação esta que faz com que κj = j.
4.3 Teste 0− 1
A descrição do funcionamento do teste 0 − 1 é feita com base no artigo (Gottwald andMelbourne, 2004), que apresenta o teste na sua forma original, e nos artigos (Gottwald andMelbourne, 2005), (Gottwald and Melbourne, 2009b) e (Gottwald and Melbourne, 2009a) queintroduziram algumas melhorias no teste original e discussões práticas sobre a implementaçãodo teste.
3O valor de υmim é da ordem do erro experimental.
58
O teste 0−1 é usado para verificar se um sistema dinâmico determinístico apresenta com-portamento regular (periódico ou quase-periódico) ou caótico. O método é aplicado diretamenteem dados de uma série temporal, logo não tem-se a necessidade de reconstrução do espaço defases. Como resultado do método, tem-se o valor K ≈ 0 se o sistema dinâmico é regular ouK ≈ 1 se for caótico.
A formulação do teste 0 − 1 é feita para tempo discreto, embora os dados experimentaisde um sistema de tempo continuo sejam discretizados através da amostragem do sinal, certoscuidados devem ser tomados ao utilizar a formulação do teste 0− 1. Pode haver para os dadosde um sistema contínuo no tempo Φ(t), um problema de sobreamostragem.
Para uma série de tempo discreta proveniente de um sistema contínuo no tempo, tem-se x(j) = Φ(tj) = Φ(jδt), sendo δt > 0 é o tempo de amostragem. A sobreamostragemocorre quando δt é muito pequeno. Isto pode levar a resultados do teste 0 − 1 (K) errados,pois a convergência de K pode se tornar lenta e como o tempo da série é finito isso se torna umproblema. Um bom método para se evitar a sobreamostragem é discretizar a série de tal formaque δt seja igual ao valor do primeiro mínimo da informação mútua da série em questão.
Dados os elementos de uma série temporal x(n) para n = 1, . . . , N temos as variáveis detranslação pc e qc:
pc(n) =n∑j=1
x(j) cos(jc), (4.35)
qc(n) =n∑j=1
x(j) sin(jc), (4.36)
para c ∈ (0, π).
As função pc(n) e qc(n) são limitadas se o sistema dinâmico for regular e caso a dinâmicaseja caótica, pc(n) e qc(n) se comportam assintoticamente com um movimento browniano 4.O crescimento das funções pc(n) e qc(n) pode ser determinado pelo deslocamento quadráticomédio (deslocamento médio quadrático) Mc(n). Para um sistema dinâmico regular pc(n) eqc(n) são limitadas assim como Mc(n). Já para o caso caótico Mc(n) cresce linearmente notempo. O deslocamento quadrático médio Mc(n) é dado por:
Mc(n) = limN→∞
1
N
N∑j=1
[pc(j + n)− pc(j)]2 + [qc(j + n)− qc(j)]2, (4.37)
que é valida quando n N ou n ≤ ncorte onde ncorte N . Na prática concluí-se que
4O movimento browniano é a movimentação de partículas macroscópicas que ocorre tipicamente em um fluidode forma aleatória em um espaço n-dimensional. Tal movimento tem comportamento caótico. Veja (Mörterset al., 2010; Addison, 1997).
59
ncorte = N/10 pode produzir bons resultados. De forma geral Mc(n) tem o aspecto abaixo:
Mc(n) = V (c)n+ Vosc(c, n) + e(c, n), (4.38)
sendo que a parcela do erro tem a propriedade limn→∞ e(c, n)/n = 0 no domínio (0, π) de c eo termo oscilatório de Mc(n) é dado por:
Vosc(c, n) = (Ex)21− cos(nc)
1− cos(c), (4.39)
com
Ex = limN→∞
1
N
N∑j=1
x(j). (4.40)
Uma melhoria na estimativa do deslocamento quadrático médio pode ser feita retirando-se o termo oscilatório. Deve ficar claro que isto não afeta as taxas de crescimento assintóticase com a retirada de Vosc(c, n), o comportamento linear de Mc(n) é regularizado. Com isso odeslocamento quadrático médio melhorado é descrito abaixo:
Dc(n) = Mc(n)− Vosc(c, n). (4.41)
Para verificar se um sistema dinâmico tem comportamento regular ou caótico, no teste0−1 emprega-se a taxa de crescimento assintóticaKc do deslocamento quadrático médioMc(n)
ou Dc(n). Um valor de Kc ≈ 0 indica dinâmica regular, e Kc ≈ 1 indica dinâmica caótica.
O valor de Kc é calculado, de forma eficiente, pelo método de correlação linear. Este mé-todo é descrito em termos deDc(n) por este apresentar melhores resultados (Gottwald and Mel-
bourne, 2009b; Gottwald and Melbourne, 2009a). Dados os vetores ξ =
1 2 . . . ncorte
Te ∆ =
Dc(1) Dc(2) . . . Dc(ncorte)
T, tem-se:
Kc = corr (ξ,∆) =cov (ξ,∆)√
var (ξ) var (∆)∈ [−1, 1], (4.42)
A covariância de dois vetores y, z de comprimento q pode ser descrita por:
cov (y, z) =1
q
q∑j=1
(yj − y)(zj − z), (4.43)
sendo a variância dada porvar (y) = cov (y,y), (4.44)
60
e a média definida por
y =1
q
q∑j=1
yj. (4.45)
Para certos valores de c no cálculo de Kc, pode ocorrer um fenômeno de ressonância, quefaz com que Kc ≈ 1. Isso ocorre independente da dinâmica ser regular ou caótica, logo, parauma dinâmica regular, tal resultado conduz a erros. Para contornar esse problema escolhe-seum conjunto de Nc valores de c de forma aleatória dentro do seu domínio e, para cada valor dec, calcula-se um Kc correspondente.
O valor de K é obtido pelo cálculo da mediana dos Nc valores de Kc ou seja:
K = mediana (Ξ), (4.46)
comΞ =
Kc(1) Kc(2) . . . Kc(Nc)
T. (4.47)
É preferível usar a mediana do que a média aritmética, pois a primeira é robusta em relaçãoaos valores de Kc relacionados à ressonância. Para reduzir ainda mais os efeitos da ressonânciao domínio de c é alterado para c ∈ (π/5, 4π/5) e foi constatado que um conjunto Nc = 100
valores distintos de c é adequado (Gottwald and Melbourne, 2009a).
No cálculo de K a série temporal precisa ter tamanho suficiente para que o valor de Kpossa convergir satisfatoriamente. No caso de “caos fraco” a série temporal deve ser mais longa,pois a convergência é mais lenta. Nas situações onde a série não seja longa o suficiente, o valorde K pode não ser o correto, pois ainda não convergiu. Neste caso, pode-se distinguir entre adinâmica regular e o caos fraco de duas formas:
• Inspeção visual do gráfico no plano p − q, na figura 4.2(a) tem-se o caso de dinâmicaregular, na figura 4.2(b) tem-se o caos fraco e na figura 4.2(c) o caso caótico. Quantomais dados a série temporal tiver, mais próximo, para o caso de caos fraco, a figura 4.2(b)estará da figura 4.2(c);
• Observar o gráfico de K versus N . Segundo (Gottwald and Melbourne, 2009a) se adinâmica é regular, K vai convergindo pra zero, já quando tem-se o caos fraco o valor deK vai aumentando lentamente.
4.4 Considerações Finais
Neste capítulo foi visto o conceito de pontos de equilíbrio e o conceito de estabilidade nosentido de Lyapunov utilizando a versão linearizada da EDO original que é valida para pontosde equilíbrio hiperbólicos.
61
−1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6−0.4
−0.2
0
0.2
0.4
0.6
0.8
1
1.2
p
q
(a) Regular, µ = 3,569945672.
−1.5 −1 −0.5 0 0.5 1−0.5
0
0.5
1
1.5
p
q
(b) Caos fraco, µ = 3,569945672 +0,001.
−50 −40 −30 −20 −10 0 10−20
−15
−10
−5
0
5
p
q
(c) Caos, µ = 3,9.
Figura 4.2: Gráfico p versus q para o mapa logístico xn+1 = µxn(1− xn) com N=5000.
Como foi visto, a ênfase maior foi dada aos métodos de obtenção dos expoentes de Lya-punov que serve para caracterizar o sistema como regular ou caótico, sendo vistos dois métodosclássicos aplicados quando se conhece a EDO do sistema, no caso os métodos vistos foram:o método padrão, que é relativamente simples e também eficiente, e o método de Eckmann-Ruelle, que serve de contraponto ao método padrão, e é um pouco mais difícil de entender eapresenta um tempo de processamento maior.
Também foram vistos métodos de análise não-linear que se aplicam a dados experimen-tais. No caso foi visto um método de obtenção do maior expoente de Lyapunov que é baseadonas ideias do método padrão, há muitos detalhes para implementar esse método como obter oatrator reconstruído. Outro método visto que se aplica a dados experimentais na forma de umasérie temporal é o teste 0 − 1 que apresenta a vantagem de não ser necessário reconstruir oatrator do sistema, sendo aplicado diretamente aos dados da série temporal, o resultado do testeé K ≈ 0 caso o sistema seja regular ou K ≈ 1 caso o sistema seja caótico.
62
Capítulo 5
Análise da Dinâmica Não-Linear doGiroscópio
Este capítulo apresenta uma análise não-linear do modelo matemático do giroscópio, as-sim como dos resultados experimentais usando séries temporais. Inicialmente é realizada umaanálise dos pontos de equilíbrio e de sua estabilidade. Em seguida, três métodos clássicos de ob-tenção de expoentes de Lyapunov são utilizados e seus resultados comparados. Estes métodossão: o método padrão, o método de Eckmann-Ruelle e o método de Wolf com séries temporais.Destaca-se que o método Padrão e método de Eckmann-Ruelle são baseados no uso da equa-ções diferencial do modelo e da matriz Jacobiana. Já o método de Wolf se baseia em uma sérietemporal que pode ser experimental ou proveniente da solução das equações do movimento. Porfim, este capítulo também utiliza o teste 0− 1. A partir destes testes é possível ter uma ideia docomportamento do sistema, obtendo informações sobre estabilidade e as condições que podemlevar a comportamento regular ou caótico.
5.1 Pontos de equilíbrio e estabilidade
Os pontos de equilíbrio são soluções no espaço de fases que se mantém constantes notempo, ou seja para uma solução da equação diferencial cuja condição inicial seja um dessespontos de equilíbrio, a solução é constante no tempo.
Para encontrar os pontos de equilíbrio do sistema deve-se resolver a equação a seguir:
F (Xe) = 0 (5.1)
com a equação diferencial do sistema na forma
X = F (X); X ∈ <6 (5.2)
63
64
sendoXe o vetor que representa os pontos de equilíbrio do sistema.
Xe =xe1 xe2 xe3 xe4 xe5 xe6
T(5.3)
Como pode ser visto na eq. (5.1), para obter os pontos de equilíbrio é necessária a equaçãodiferencial do sistema no espaço de estados. Para facilitar o cálculo dos pontos de equilíbrio,esta equação, que foi obtida no capítulo 3, é representada abaixo.
x1 = x2 (5.4)
x2 =τz + [(−Ix + Iy + Iz)x2 sin(x3) + Ixxx6]x4
Iz cos(x3)(5.5)
x3 = x4 (5.6)
x4 =τy − Ixxx2x6 cos(x3) + (Ix − Iz)x22 cos(x3) sin(x3)
Iy(5.7)
x5 = x6 (5.8)
x6 = τx+(Ix+Iy−Iz)x2x4 cos(x3)Ixx
+ Ixτz+[(−Ix+Iy+Iz)x2 sin(x3)+Ixxx6]x4
IxxIztan(x3) (5.9)
Resolvendo a eq. (5.1), de imediato com as eqs. (5.4), (5.6) e (5.8), obtêm-se:
xe2 = xe4 = xe6 = 0 (5.10)
Das eqs. (5.5), (5.7) e (5.9), juntamente com o resultado da eq. (5.10), obtém-se o restante dasolução de equilíbrio.
Da eq. (5.5) tem-se:τz
Iz cos(xe3)= 0 (5.11)
ou seja,τz = 0 (5.12)
Agora da eq. (5.7) obtém-se:τyIy
= 0 (5.13)
logoτy = 0 (5.14)
65
Finalmente, da eq. (5.9):τxIxx
+Ixτz tan(xe3)
IxxIz= 0 (5.15)
Usando o resultado da eq. (5.12) em (5.15), resulta
τx = 0 (5.16)
Desta forma, para se encontrar os pontos de equilíbrio, são necessárias as equações dostorques, as quais são mostradas abaixo e foram obtidas no capítulo 3:
τx
τy
τz
=
−Tx sgn (xe6)− bxxe6 + [Tz sgn (xe2) + bxxe2 ] sin(xe3)
MgH cos(xe3)− Ty sgn (xe4)− byxe4 + τey (xe3)
−[Tz sgn (xe2) + bzxe2 ] cos(xe3)
(5.17)
com
τey (xe3) =
−k(xe3 − θi) : xe3 < θi
0 : θi ≤ xe3 ≤ θf
−k(xe3 − θf ) : xe3 > θf
(5.18)
Usando os resultados das eqs. (5.12), (5.14), (5.16) e (5.10) na eq. (5.17), tem-se:
MgH cos(xe3) + τey (xe3) = 0 (5.19)
onde observa-se que xe1 e xe5 não tem restrições, portanto:
xe1 , xe5 ∈ < (5.20)
Para encontrar xe3 , deve-se resolver a eq. (5.19). Os valores deM e g são constantes, jáHpode ser modificado conforme se mudam as posições dos contrapesos. O valor de xe3 dependede H , ou seja, xe3 = xe3(H). A seguir são obtidas as soluções de xe3 para os seguintes casos:
1. posição do centro de massa H < 0: neste caso τey (xe3) > 0 e xe3 < θi e a eq. (5.19) fica
MgH cos(xe3)− k(xe3 − θi) = 0 (5.21)
Para encontrar xe3 quando H < 0, basta resolver a eq. não-linear (5.21).
2. posição do centro de massa H > 0: para este caso τey (xe3) < 0 e xe3 > θf . Partindo daeq. (5.19) chega-se na eq. (5.22)
MgH cos(xe3)− k(xe3 − θf ) = 0 (5.22)
Neste caso xe3 é encontrado resolvendo a eq. não-linear (5.22).
66
3. posição do centro de massa H = 0: nessa situação τey (xe3) = 0 e θi ≤ xe3 ≤ θf . Daeq. (5.19) conclui-se que não há restrições no valor de cos(xe3), logo xe3 é dado por
θi ≤ xe3 ≤ θf (5.23)
Como pode ser observado dos cálculos anteriores, há infinitos pontos de equilíbrio para o sis-tema giroscópico.
A análise de estabilidade dos pontos de equilíbrio é feita pela determinação dos autovalo-res γ obtidos da solução do determinante a seguir:
|J(Xe)− γI| = 0 (5.24)
sendo J(Xe) a matriz Jacobiana, obtida a partir de F (X) da equação diferencial e calculadano ponto de equilíbrio, e I a matriz identidade.
AplicandoXe na matriz Jacobiana obtida no apêndice B, obtém-se:
J(Xe) =
0 1 0 0 0 0
0 − bzIz
0 0 0 0
0 0 0 1 0 0
0 0 −MgH sin(xe3)+deIy
− byIy
0 0
0 0 0 0 0 1
0(bx− bzIx
Iz) sin(xe3 )
Ixx0 0 0 − bx
Ixx
(5.25)
sendo
de =
−k : x3 < θi
0 : θi ≤ x3 ≤ θf
−k : x3 > θf
(5.26)
O determinante, lado esquerdo da eq. (5.24), tem a forma:
|J(Xe)−γI| =
∣∣∣∣∣∣∣∣∣∣∣∣∣∣∣
−γ 1 0 0 0 0
0 − bzIz− γ 0 0 0 0
0 0 −γ 1 0 0
0 0−MgH sin(xe3 )+de
Iy− byIy− γ 0 0
0 0 0 0 −γ 1
0(bx− bzIx
Iz) sin(xe3 )
Ixx0 0 0 − bx
Ixx− γ
∣∣∣∣∣∣∣∣∣∣∣∣∣∣∣(5.27)
Resolvendo a eq. (5.24) com a ajuda da eq. (5.27) tem-se:
γ2(bzIz
+ γ
)(bxIxx
+ γ
)[γ
(byIy
+ γ
)+MgH sin(xe3)− de
Iy
]= 0 (5.28)
67
Assim, a solução da eq. (5.28) é:γ1 = γ2 = 0 (5.29)
γ3 = −bzIz
(5.30)
γ4 = − bxIxx
(5.31)
γ5 e γ6 são obtidos da solução da eq. (5.32), que é parte da eq. (5.28),
γ2 +byIyγ +
MgH sin(xe3)− deIy
= 0, (5.32)
resolvendo-se para os casos em que H = 0, H < 0 e H > 0.
No caso da posição do centro de massaH = 0, θi ≤ xe3 ≤ θf , logoMgH sin(xe3)−de =
0 e a eq. (5.32) torna-se:
γ2 +byIyγ = 0 (5.33)
com as soluções dadas por:γ5 = 0 (5.34)
γ6 = −byIy
(5.35)
Algumas observações podem ser feitas para os casos:
• para H < 0, xe3 < θi e sin(xe3) < 0 logo MgH sin(xe3)− de = MgH sin(xe3) + k > 0;
• para H > 0, xe3 > θf e sin(xe3) > 0 logo MgH sin(xe3)−de = MgH sin(xe3)+k > 0;
Consequentemente para H < 0 ou H > 0, a eq. (5.32) se torna:
γ2 +byIyγ +
MgH sin(xe3) + k
Iy= 0 (5.36)
com a solução na forma:
γ =1
2
byIy
(−1±
√δ)
(5.37)
sendo
δ = 1− 4
(Iyby
)2MgH sin(xe3) + k
Iy(5.38)
Existe três casos a se considerar δ = 0, δ < 0 e δ > 0.
1. δ = 0: Tem-se IyMgH sin(xe3 )+k
b2y= 1
4e as soluções são números reais negativos e iguais,
68
ou seja, a solução tem multiplicidade 2
γ5 = γ6 = −1
2
byIy
(5.39)
2. δ < 0: Tem-se IyMgH sin(xe3 )+k
b2y> 1
4e as soluções γ5 e γ6 são números complexos
conjugados, com parte real negativa
γ5 =1
2
byIy
(−1 + i
√−δ)
(5.40)
γ6 =1
2
byIy
(−1− i
√−δ)
(5.41)
3. δ > 0: Não é difícil ver que 0 < δ < 1 e IyMgH sin(xe3 )+k
b2y< 1
4com as soluções γ5 e γ6
sendo números reais distintos e negativos
γ5 =1
2
byIy
(−1 +
√δ)< 0 (5.42)
γ6 =1
2
byIy
(−1−
√δ)< 0 (5.43)
Estes resultados, a respeito dos sinais dos autovalores do sistema linearizado, estão natabela 5.1.
Tabela 5.1: Resumo dos sinais dos autovalores dos pontos de equilíbrio do modelo do giroscó-pio.
Autovalor H < 0 ou H > 0 H = 0δ < 0 δ = 0 δ > 0
γ1 0 0 0 0γ2 0 0 0 0
γ3 − bzIz
− bzIz
− bzIz
− bzIz
γ4 − bxIxx
− bxIxx
− bxIxx
− bxIxx
γ512
byIy
(−1 + i
√−δ)−1
2
byIy
12
byIy
(−1 +
√δ)
0
γ612
byIy
(−1− i
√−δ)−1
2
byIy
12
byIy
(−1−
√δ)
− byIy
Como dois dos autovalores sempre são nulos, os pontos de equilíbrio são elípticos e por-tanto a análise do sistema não-linear através do estudo dos autovalores do sistema linearizadonão é conclusiva. Para saber a estabilidade dos ponto de equilíbrio do sistema não-linear podese usar o teorema da variedade central ou o método direto de Lyapunov, porem estes fogem doescopo desta dissertação.
69
5.2 Expoentes de Lyapunov
Com os expoentes de Lyapunov, como visto no capítulo 4, é possível caracterizar um sis-tema não-linear, isto é, saber se o mesmo tem comportamento regular ou caótico, ou em outraspalavras verificar a sensibilidade do sistema às condições iniciais. No sistema giroscópico emestudo nesta dissertação tem-se um total de seis expoentes de Lyapunov, sendo um para cadauma das direções canônicas do espaço de fases. Para um sistema ser caótico, basta que umdestes expoentes seja positivo, por isso, e por simplificação, nos gráficos onde os expoentesaparecem, somente contém os três maiores expoentes.
O comportamento do giroscópio pode ser modificado com as mudanças dos valores dosparâmetros associados ao sistema. Isso pode ser verificado através da mudança dos valores dosexpoentes de Lyapunov. A fim de ter um estudo mais completo, esta dissertação adotou umafaixa de valores para alguns desses parâmetros.
Dentre os parâmetros do modelo do giroscópio têm-se os associados ao amortecimentoviscoso (bx, by e bz), ao torque do atrito cinético (Tx, Ty e Tz), torque elástico (k) e os re-lacionados as componentes do tensor de inércia, entre outros. No caso do tensor de inércia,suas componentes podem ser modificadas facilmente com as alterações das posições dos con-trapesos, que neste caso seriam o contrapeso maior com massa m1 e posição h1 e contrapesomenor com massa m2 e posição h2. A mudança das posições dos contrapesos pode alterar tam-bém o amortecimento viscoso do sistema, mas por simplificação admite-se que este parâmetropermanece constante em relação às mudanças das posições dos contrapesos.
Como o modelo do giroscópio contém muitos parâmetros 1, foram escolhidos apenas doispara estudo, o amortecimento viscoso e as posições dos contrapesos que, de certa forma, podemser modificados experimentalmente. O amortecimento viscoso, pode ser modificado na prática,por exemplo, com a mudança da densidade do ar. Mas como seria muito difícil controlar esseparâmetro experimentalmente, só é mostrada a análise dos expoentes de Lyapunov relacionadosao modelo.
Para o amortecimento viscoso foi escolhida uma faixa de valores que no caso, está emtorno dos valores estimados experimentalmente. A faixa de valores do amortecimento viscoso(em N ·m · s/rad) está entre 1,0 × 10−5 e 1,0 × 10−3 com passo de 1,0 × 10−5. Para queo tempo de processamento na obtenção dos expoentes de Lyapunov, para este caso, não fiquemuito alto, os parâmetros relacionados ao torque cinético não foram considerados, ou seja,assumiu-se Tx = Ty = Tz = 0,0.
Quanto a posição dos contrapesos, apenas a do contrapeso maior será variada para a aná-lise tando do modelo quanto do experimento. O contrapeso m2 é fixo o mais próximo possíveldo centro, ou seja, h2 = −0,0827 m. A posição do contrapeso maior h1, varia de −0,2579 m,
1Com valores destes parâmetros obtidos no apêndice A.
70
que corresponde a distância de 15,0 cm entre as faces mais próximas dos contrapesos m1 e m2;até−0,1079 m, neste caso os contrapesos estavam encostados. Esta variação ocorreu com passode 1,0 mm quando a análise foi do modelo e 1,0 cm, quando a análise foi experimental.
O tempo total escolhido para obtenção dos expoentes de Lyapunov do modelo foi de600,0 s. Este tempo foi escolhido de forma que fosse suficiente para que a solução do sistematenha uma grande quantidade de pontos no espaço de fases. O tempo é em média dez vezesmaior que o tempo adotado na parte experimental. Nos experimentos não foi possível ir muitoalém de 1,0 min, pois o cabo que liga o sensor de nutação acabava influenciando no movimento,não sendo possível prolongar muito o experimento. Caso os resultados simulados sejam com-parados com os resultados experimentais, o tempo total sera de 60,0 s. A condição inicial X0,tanto da simulação quanto do experimento, pode ser vista na tabela 5.2.
Tabela 5.2: Condições iniciais para estimativa dos expoentes de Lyapunov.x1 (rad) x2 (rad/s) x3 (rad) x4 (rad/s) x5 (rad) x6 (rad/s)
X0 0,00 0,00 0,00 0,00 0,00 68,49
X′
0 0,01 0,01 0,01 0,01 0,01 68,50
O erro padrão adotado para o método de Runge-Kutta de oitava ordem de passo variável,que foi utilizado nos programas usados na seções 5.2.1 e 5.2.2, possui o valor de 1,0 × 10−12,tanto para o erro relativo quanto para o erro absoluto, que é próximo do limite de casas decimaispara números de precisão dupla; qualquer exceção é explicitada no texto.
A fim de comparar duas soluções próximas no espaço de fases, adotou-se uma segundacondição inicialX
′
0, com diferença de 0,01 rad (ou rad/s), que é o valor aproximado da incer-teza do ângulo, tabela 5.2.
5.2.1 Expoentes de Lyapunov do modelo
A seguir é feita uma descrição dos programas e algoritmos usados para estimar os ex-poentes de Lyapunov. Todos estes programas foram testados com o sistema de Lorentz, paraposterior utilização no sistema giroscópico. Na sequência são obtidos os expoentes de Lyapu-nov para cada método. Comparações entre os métodos assim como visualização do espaço defases para alguns casos, também são feitos.
Para o método padrão foi necessário além da equação diferencial, obtida no capítulo 3, amatriz Jacobiana do sistema, que é apresentada no apêndice B. Para o método padrão foi criadoum programa baseado no algoritmo descrito em (Wolf et al., 1985). Neste programa, para obteras soluções numéricas das equações diferenciais foi usado o método de Runge-Kutta de oitavaordem de passo variável, poderia ser usado o de quarta ordem, mas o importante é o fato dopasso ser variável ou seja, o passo de integração é modificado pelo algoritmo do método para
71
manter a solução numérica dentro de um erro estipulado. Outros parâmetros relacionados aoalgoritmo são:
• Condições iniciais da equação diferencial;
• Tempo total de integração;
• Tempo entre ortonormalização (passo). Este passo foi escolhido como 0,25 s.
Para o método de Eckmann-Ruelle, foi adotado o programa desenvolvido por Paul HenryBryant, cujo código se encontra em (Bryant, 2009). Neste programa foram feitas algumasmodificações no código original, como a possibilidade de gravar os expoentes em arquivo, usode logaritmo na base dois para compatibilidade com o método padrão e utilização do métodode Runge-Kutta de oitava ordem com passo variável. Para esse método, assim como no métodopadrão, é necessária a equação diferencial e a matriz Jacobiana do sistema 2. Os parâmetrosusados pelo programa, como a condição inicial, são os mesmos adotadas pelo método padrãopara que se possa comparar os métodos.
No método de Wolf que utiliza séries temporais, é necessário, além da série temporal, opasso, o tempo de amostragem, o tempo de atraso e a dimensão de imersão. Com isso pode-se obter o atrator reconstruído, e assim, poder aplicar o método. Para a obtenção do tempode atraso foi usado o método de informação mútua. Já para obter a dimensão de imersão, foiusado o método dos falsos vizinhos próximos; em ambos os casos, foi usado os programasdesenvolvidos por Matjaž Perc, que se encontra em (Perc, 2007). Aplicações dos programaspodem ser vistos em (Perc, 2005a; Perc, 2005b; Kodba et al., 2005; Perc, 2006). No programaresponsável por encontrar o tempo de atraso, denominado de mutual, é necessário o número depontos da série. Como existe um limite de 250000 pontos no programa mutual, foi necessárioutilizar um tempo de amostragem de 1,0 × 10−2 para não exceder esse limite; para o tempode atraso máximo foi usado o valor 500(5,0 s). Os outros parâmetros do programa não foramalterados. Agora no programa responsável por encontrar a dimensão de imersão, denominadode fnn, foi necessário indicar a dimensão máxima e mínima, as quais foram ajustadas como 1 e10, respectivamente. Outros parâmetros como número de pontos da série e o tempo de atrasosão obtidos do programa mutual; os parâmetros restantes não foram alterados.
Através do atrator reconstruído, pode-se obter o maior expoente de Lyapunov. Para istofoi usado o programa desenvolvido por Matjaž Perc nomeado de lyapmax e também foi criadoum programa chamado de wolf. Os programas lyapmax e wolf tem como base o algoritmodescrito em (Wolf et al., 1985). Ambos os programas usam os dados obtidos dos programasmutual e fnn, além dos seguintes parâmetros: passo 25 (0,25 s), distância mínima inicial 0,01,distância máxima inicial 2,0, com os demais parâmetros deixados nos valores padrões.
2No algoritmo existe a possibilidade de se usar a matriz Jacobiana numérica obtida, nesse caso, da própriaequação diferencial.
72
A convergência dos dois maiores expoentes de Lyapunov para o método padrão, para ocaso em que a face do contrapeso m1 se encontra na posição 9,0 cm da face de m2, pode serobservada na (figura 5.1(a)). O primeiro expoente tem um convergência assintótica (linha (1) dafigura 5.1(a)), já o segundo expoente (linha (2) da figura 5.1(a)) é assintótica com oscilações. Nocaso do segundo expoente (figura 5.1(b), linha (2)), devido às oscilações, dependendo do tempofinal adotado o valor do expoente pode ser muito diferente. Por exemplo se o tempo final for de596,0 s, o expoente valerá−0,95×10−3 1/s, já se o tempo final for de 597,0 s o expoente será de3,05× 10−3 1/s, o que constitui uma diferença relativa muito grande no expoente de Lyapunovpor uma diferença de apenas 1,0 s no tempo final. Para contornar esse problema, o expoentefoi calculado como a média aritmética simples do último trecho do gráfico da convergência dosexpoentes e foi adotado o trecho final dos últimos 5,0 s que equivale às últimas 20 amostras dográfico (já que uma aproximação dos expoentes ocorre a cada 0,25 s, que é o valor do passoadotado). Neste casso o valor médio do segundo expoente de Lyapunov para um tempo totalde 600,0 s é m = 1,5813 × 10−3 1/s (figura 5.1(b), linha (m)) com desvio padrão da médiad = 3,8116 × 10−4 1/s e erro relativo e = 24,1 %. Adotando-se um trecho maior esse errorelativo pode ficar menor, mas a média não vai mudar muito.
0 100 200 300 400 500 600
0
0.02
0.04
0.06
0.08
0.1
0.12
0.14
0.16
0.18
0.2
tempo [s]
Exp
oen
tes
de
Lya
pu
no
v [1
/ s]
(1)(2)
(1)(2)
(a) Convergência de λ1 e λ2.
594 595 596 597 598 599 600−4
−3
−2
−1
0
1
2
3
4x 10
−3
tempo [s]
Exp
oent
es d
e Ly
apun
ov [1
/ s]
(2)(m)
(m)
(2)
(b) Trecho final da convergência de λ2.
Figura 5.1: Exemplo de convergência para os dois maiores expoentes de Lyapunov usando ométodo padrão.
Um exemplo de gráfico de convergência para o método de Eckmann-Ruelle e para ométodo de Wolf para séries temporais simuladas, nas mesmas condições que foram aplicadasno método padrão, pode ser visto nas figuras 5.2(a) e 5.2(b), respectivamente.
O tempo de processamento para dois dos três exemplos de convergência dos expoentesfoi: método padrão 68,0 s, método de Eckmann-Ruelle 173,0 s; já era esperado que o método deEckmann-Ruelle demandasse mais tempo de processamento, como visto em (Ramasubramanianand Sriram, 2000).
Os expoentes de Lyapunov mostrados nas figuras desta dissertação, são os três maiores, os
73
300 350 400 450 500 550 600−0.05
0
0.05
0.1
0.15
0.2
0.25
0.3
tempo [s]
Exp
oent
es d
e Ly
apun
ov [1
/ s]
(1)(2)
(1)
(2)
(a) Método de Eckmann-Ruelle.
0 100 200 300 400 500 600
−0.2
−0.1
0
0.1
0.2
0.3
0.4
tempo [s]
Exp
oent
es d
e Ly
apun
ov [1
/ s]
(1)(2)
(1)
(2)
(b) Método de Wolf. Série da velocidade de nuta-ção (x4). Linha (1) programa lyapmax e na linha (2)programa wolf.
Figura 5.2: Exemplo de convergência para os dois maiores expoentes de Lyapunov usando osmétodos de Eckmann-Ruelle e Wolf.
restantes em geral, são negativos. No caso isso é feito para que as figuras não fiquem carregadas,já que há um total de seis expoentes.
Iniciando a análise dos expoentes de Lyapunov do modelo, com a variação dos parâmetrosdo amortecimento viscoso, no caso bx, by e bz conforme descrito na seção 5.2. Variou-se umdos parâmetros do amortecimento viscoso de cada vez, mantendo-se os demais com os valorespadrões. Para o controle do erro no método de Runge-Kutta foram usados os seguintes valores:erro relativo 1,0 × 10−12, erro absoluto 1,0 × 10−11; o erro absoluto foi reduzido para que ométodo de Runge-Kutta não apresenta-se mensagens de “erro”, como por exemplo, que o passode integração estivesse muito pequeno, etc.
As figuras 5.3(a) e 5.3(b) mostram os expoentes de Lyapunov calculados com o métodopadrão e com o método de Eckmann-Ruelle, respectivamente, considerando a variação do parâ-metro bx. No método padrão na figura 5.3(a), observou-se que sempre um dos expoentes, o (1)no caso, é sempre zero, os demais expoentes são positivos para a maioria dos parâmetros. Emrelação ao método de Eckmann-Ruelle, os expoentes são mostrados na figura 5.3(b). Observa-seque para toda a faixa de parâmetros os expoentes são positivos. Não há concordância numéricasentre os resultados do método padrão e do método de Eckmann-Ruelle. Por exemplo, enquantoum dos expoentes do método padrão é sempre zero, no método de Eckmann-Ruelle isso nãoocorre.
Na figura 5.4(a) podem ser vistos os expoentes estimados pelo método padrão, em todaa faixa de by. Há pelo menos um dos expoentes positivos e, de forma geral, esses expoentesdecaem com o aumento de by. Já no método de Eckmann-Ruelle, figuras 5.4(b), dois dosexpoentes ficam praticamente constantes, o restante decai com o aumento de by.
Quanto a influência da variação de bz, pode ser constatado analisando a figura 5.5(a)
74
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
x 10−3
0
0.005
0.01
0.015
0.02
0.025
0.03
0.035
Parâmetro bx [N ⋅ m ⋅ s / rad]
Exp
oent
es d
e Ly
apun
ov [1
/ s]
(1)(2)(3)
(1)
(2)
(3)
(a) Métodos padrão.
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
x 10−3
0
0.005
0.01
0.015
0.02
0.025
0.03
0.035
Parâmetro bx [N ⋅ m ⋅ s / rad]
Exp
oent
es d
e Ly
apun
ov [1
/ s]
(1)(2)(3)
(1)
(2)(3)
(b) Métodos Eckmann-Ruelle.
Figura 5.3: Expoentes de Lyapunov pelos métodos padrão e de Eckmann-Ruelle. A variação debx foi de 1,0× 10−5 até 1,0× 10−3.
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
x 10−3
−2
0
2
4
6
8
10
12
14
16
18
x 10−3
Parâmetro by [N ⋅ m ⋅ s / rad]
Exp
oent
es d
e Ly
apun
ov [1
/ s]
(1)(2)(3)
(1)
(3)
(2)
(a) Métodos padrão.
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
x 10−3
−2
0
2
4
6
8
10
12
14
16
18
x 10−3
Parâmetro by [N ⋅ m ⋅ s / rad]
Exp
oent
es d
e Ly
apun
ov [1
/ s]
(1)(2)(3)
(1)
(3)
(2)
(b) Métodos Eckmann-Ruelle.
Figura 5.4: Expoentes de Lyapunov pelos métodos padrão e de Eckmann-Ruelle. A variação deby foi de 1,0× 10−5 até 1,0× 10−3.
que, usando o método padrão o expoentes (3) apresentou um valor máximo e um mínimo, masainda assim positivo, assim como os outros expoentes. Já para o método de Eckmann-Ruellefiguras 5.5(b), houve muita oscilação nos expoentes, mas sempre com pelo menos um positivo,em toda a faixa de valores de bz.
Para finalizar a análise do amortecimento viscoso, nas figuras 5.6(a), 5.6(b) e 5.6(c) foiefetuada a soma dos expoentes de Lyapunov para os parâmetros bx, by e bz, respectivamente.Observa-se nas figuras 5.6(a) e 5.6(b) que, independente do método aplicado, o coeficiente dasretas foi sempre o mesmo e, conforme aumenta-se o parâmetro de amortecimento viscoso, pro-porcionalmente a soma dos expoentes diminui. Neste caso, a soma dos expoentes é negativa,pois o sistema é dissipativo. Já para bz, conforme mostra a figura 5.6(c), diferente dos casos
75
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
x 10−3
−0.005
0
0.005
0.01
0.015
0.02
0.025
0.03
0.035
0.04
Parâmetro bz [N ⋅ m ⋅ s / rad]
Exp
oent
es d
e Ly
apun
ov [1
/ s]
(1)(2)(3)
(1)
(2)
(3)
(a) Métodos padrão.
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
x 10−3
−0.005
0
0.005
0.01
0.015
0.02
0.025
0.03
0.035
0.04
Parâmetro bz [N ⋅ m ⋅ s / rad]
Exp
oent
es d
e Ly
apun
ov [1
/ s]
(1)(2)(3)
(1)
(2)
(3)
(b) Métodos Eckmann-Ruelle.
Figura 5.5: Expoentes de Lyapunov pelos métodos padrão e de Eckmann-Ruelle. A variação debz foi de 1,0× 10−5 até 1,0× 10−3.
anteriores, com o aumento do parâmetro de amortecimento correspondente, no caso do parâme-tro bz, a soma dos expoentes fica oscilando, principalmente no método de Eckmann-Ruelle. Namédia a soma fica mais negativa com aumento de bz.
1 2 3 4 5 6 7 8 9 10
x 10−4
−0.18
−0.17
−0.16
−0.15
−0.14
−0.13
−0.12
−0.11
−0.1
−0.09
−0.08
Parâmetro bx [N ⋅ m ⋅ s / rad]
Som
a do
s E
xpoe
ntes
de
Lyap
unov
[1 /
s]
(1)(2)
(1)
(2)
(a) Variação de bx.
1 2 3 4 5 6 7 8 9 10
x 10−4
−0.1
−0.098
−0.096
−0.094
−0.092
−0.09
−0.088
−0.086
−0.084
−0.082
−0.08
Parâmetro by [N ⋅ m ⋅ s / rad]
Som
a do
s E
xpoe
ntes
de
Lyap
unov
[1 /
s]
(1)(2)
(2)
(1)
(b) Variação de by .
1 2 3 4 5 6 7 8 9 10
x 10−4
−0.022
−0.02
−0.018
−0.016
−0.014
−0.012
−0.01
−0.008
−0.006
Parâmetro bz [N ⋅ m ⋅ s / rad]
Som
a do
s E
xpoe
ntes
de
Lyap
unov
[1 /
s]
(1)(2)
(2)
(1)
(c) Variação de bz .
Figura 5.6: Soma dos expoentes de Lyapunov. Métodos: padrão (1), Eckmann-Ruelle (2). Comparâmetros variando de 1,0× 10−5 até 1,0× 10−3.
Apesar dos expoentes de Lyapunov para os métodos padrão e de Eckmann-Ruelle nãoconcordarem na análise da influência do amortecimento viscoso, a soma dos expoentes diferepouco, indicando que as orientações dos expoentes é que diferem. Como a soma dos expoentestem relação com a expansão ou contração do volume do atrator no espaço de fases, uma somanegativa indica contração do volume do atrator, isso já era esperado pois o sistema é dissipativoe quanto mais dissipativo mais o volume se contraí, como foi observado na diminuição da somados expoentes conforme o amortecimento viscoso aumentava.
Como os expoentes de Lyapunov obtidos pelos métodos padrão e de Eckmann-Ruelle nãoconcordaram, poder-se-ia imaginar que os programas que implementam tais métodos estives-
76
sem com algum erro. Porém tais programas foram testados como o sistema de Lorentz e osresultados foram plenamente satisfatórios como pode ser visto no apêndice C (página 103).
No método padrão 3, as aproximações dos expoentes de Lyapunov são obtidas pela evolu-ção no tempo de n condições iniciais que ocorrem conforme a solução das equações diferenciaislinearizadas do sistema. Em determinados instantes de tempo é necessário fazer uma ortonor-malização na base do sistema linearizado, usando o processo de Gram-Schmidt. O primeirovetor, conforme evolui no tempo, vai tendendo à uma direção que tende a maximizar o expo-ente de Lyapunov correspondente. O processo de Gram-Schmidt normaliza o primeiro vetormas não altera a sua orientação, já os demais vetores mudam sua orientação de forma a ficaremortonormais. Logo, no método padrão, a direção no espaço de fases em que há a maior expan-são 4 ou menor contração ocorre na direção do primeiro vetor. Contudo o que se observou éque o expoente de Lyapunov associado a essa direção é nulo, ou seja o vetor começa unitárioem uma certa direção e com o passar do tempo não muda, fica constante. Devido ao primeirovetor não ser orientado na direção de maior expansão ocorrem diferenças nos resultados dos ex-poentes de Lyapunov calculados pelo método padrão e pelo método de Eckmann-Ruelle. Pararesolver esse problema é necessário alterar o programa de forma que a orientação do primeirovetor seja diferente.
O último parâmetro estudado foi relacionado a posição do contrapeso m1, veja descriçãomais detalhada no inicio desta seção 5.2. Nas figuras 5.7(a) e 5.7(b) apresenta-se os expoentesde Lyapunov estimados usando o método padrão e o método de Eckmann-Ruelle, respectiva-mente. Observou-se que um dos expoentes decresce de valor, conforme o giroscópio fica maisdesequilibrado. Em nenhum momento os expoentes ficaram negativos. A faixa de parâmetrosde h1 não continha nenhum valor que fizesse com que o centro de massa ficasse zero, ou seja,H = 0,0 m, mas o extremo h1 = −0,2579 m esta próximo do equilíbrio. Logo, por simetria 5,o valor máximo do expoente de Lyapunov, linha (2), deve estar na posição h′
1 que faz H = 0 oudeve haver dois máximos em torno de h′
1.
Da análise feita nos parâmetros relacionados ao amortecimento viscoso e a posição docontrapeso maior pode-se ver que a ordem dos expoentes em relação a suas amplitudes sealterna, não somente o primeiro como já foi discutido para o método padrão. Isso acontece paraos dois métodos, mas no sistema de Lorentz isso não ocorreu. Também é observado que osexpoentes encontrados para o sistema giroscópico, embora sejam positivos, tem um amplitudebaixa com um valor máximo da ordem de 0,018 e 0,035 para os métodos padrão e de Eckmann-Ruelle respectivamente. Com um valor baixo para os expoentes de Lyapunov mais o fato de aordem dos expoentes ficar, em muitos casos, alternando, conclui-se que esses expoentes podemestar convergindo para zero. Caso os expoentes de Lyapunov sejam realmente positivos, tem-senessa situação caos fraco já que os expoentes tem um valor pequeno mas positivo, por outro
3Ver seções 4.2.1 e 4.2.2, páginas 52 e 55 respectivamente.4Maior divergência entre duas trajetórias com condições iniciais próximas.5Essa seria uma hipótese a confirmar.
77
−0.26 −0.24 −0.22 −0.2 −0.18 −0.16 −0.14 −0.12 −0.1−2
0
2
4
6
8
10
12
14
16
18
x 10−3
Parâmetro h1 [m]
Exp
oent
es d
e Ly
apun
ov [1
/ s]
(1)(2)(3)
(1)
(2)
(3)
(a) Método padrão.
−0.26 −0.24 −0.22 −0.2 −0.18 −0.16 −0.14 −0.12 −0.1−2
0
2
4
6
8
10
12
14
16
18
x 10−3
Parâmetro h1 [m]
Exp
oent
es d
e Ly
apun
ov [1
/ s]
(1)(2)(3)
(1)
(2)
(3)
(b) Método de Eckmann-Ruelle.
Figura 5.7: Expoentes de Lyapunov com h1 variando de −0,2579 m até −0,1079 m.
lado se os expoentes estiverem convergindo para zero, tem-se um comportamento regular.
Para confirmar a suspeita de que os expoente de Lyapunov positivos estão convergindopara zero, pretende-se usar o teste 0 − 1 que possibilita ter uma ideia se é caos fraco ou se ocomportamento é regular. Para o teste são escolhidos dois expoentes, em relação a variação doparâmetro bx que são obtidos pelo método padrão e que possuam as maiores amplitudes. Nafigura 5.3(a), pode-se ver que o expoente (3) tem um valor máximo em bx = 2,2 × 10−4, já oexpoente (2), apresenta um pico em bx = 7,3× 10−4.
Com os parâmetros bx = 2,2× 10−4 e bx = 7,3× 10−4 obtêm-se duas séries temporais daposição angular de nutação x3. Tal escolha de variável de estado é justificada na próxima seção,que é obtida pela solução numérica da EDO do sistema. A condição inicial é a mesma adotadano cálculo dos expoentes de Lyapunov e o passo de integração foi δt = 1,0× 10−2 s e o tempototal foi t = 800,0 s. A partir das séries temporais os tempos de atraso são obtidos e o teste 0−1
pode ser aplicado aos dados da série que foram re-amostrados utilizando o tempo de atraso 6.Para cada parâmetro dois valores de K são obtidos, um que foi obtido usando a série temporalna íntegra (Kt) e outro que exclui da série os primeiros 10 % dos dados (Kr, série reduzida),o que elimina os movimentos iniciais de maior amplitude. Na tabela 5.3 pode-se observar quetanto Kt quanto Kr, para os dois parâmetros usados, apresenta um valor próximo de zero. Paraconfirmar se o comportamento é regular ou se há caos fraco são utilizados a seguir os gráficosp− q.
Para o parâmetros bx = 2,2 × 10−4, usando todos os dados da série, obtém-se a fi-gura 5.8(a). Já a figura 5.8(b) foi obtida retirando-se 10 % dos dados iniciais da série. Dasfiguras 5.8(a) e 5.8(b) conclui-se então que o comportamento para bx = 2,2 × 10−4 é regular,conforme discussão no final da seção 4.3.
6Detalhes de como o teste 0− 1 funciona pode ser visto na seção 4.3, página 57.
78
Tabela 5.3: Resumo do teste 0− 1 aplicado a duas séries temporais simulada de x3. Sendo τt eKt o tempo de atraso e o resultado do teste 0− 1 da série completa respectivamente. Enquantoτr e Kr estão relacionados a série original em que foi eliminado os primeiros 10 % dos dados.
bx (N ·m · s/rad) τt τr Kt Kr
2,2× 10−4 9 7 2,99×10−1 1,36×10−1
7,3× 10−4 10 7 2,31×10−1 −6,22×10−4
−0.5 0 0.5−0.5
−0.4
−0.3
−0.2
−0.1
0
0.1
0.2
0.3
0.4
0.5
p
q
(a) Série completa.
−0.8 −0.6 −0.4 −0.2 0 0.2−0.2
−0.1
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
p
q
(b) Série sem os primeiros 10% dos dados.
Figura 5.8: Gráfico p− q do teste 0− 1 com bx = 2,2× 10−4 N ·m · s/rad.
Para o parâmetros bx = 7,3 × 10−4 com a série temporal completa observa-se na fi-gura 5.9(a) que existe alguns pontos um pouco desalinhados. Isso acaba gerando um pouco dedúvidas em relação se é caos fraco ou o comportamento regular. Agora com a série reduzida,isto é, sem os movimentos iniciais com amplitudes maiores, pode-se ver da figura 5.9(b) que ocomportamento é regular.
−0.6 −0.4 −0.2 0 0.2 0.4−0.4
−0.3
−0.2
−0.1
0
0.1
0.2
0.3
0.4
0.5
0.6
p
q
(a) Série completa.
−1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4−0.2
0
0.2
0.4
0.6
0.8
1
1.2
p
q
(b) Série sem os primeiros 10% dos dados.
Figura 5.9: Gráfico p− q do teste 0− 1 com bx = 7,3× 10−4 N ·m · s/rad.
Como o comportamento é regular para os dois casos estudados, pelo menos desconside-rando os dados iniciais da série que formaria provavelmente a parte de transiente, conclui-seque os expoentes positivos encontrados estão convergindo para zero.
79
Para finalizar a análise não-linear do modelo são apresentados os gráficos no espaço defases para os casos bx = 2,2× 10−4 e bx = 7,3× 10−4. Como a dimensão do espaço de fases é6 surge a necessidade de projetar a solução da EDO em um espaço de dimensão menor. Nestecaso adota-se a opção de ignorar 3 dimensões, no caso os estados mantidos são aqueles cujas asamplitudes são limitadas, ou seja, a velocidade angular de precessão x2, a posição angular denutação x3 e a velocidade de nutação x4.
A solução da EDO para bx = 2,2 × 10−4 dos primeiros 30,0 s pode ser vista na fi-gura 5.10(a) e os ultimos 30,0 s de um tempo total de 300,0 s é mostrado na figura 5.10(b). Jápara bx = 7,3 × 10−4 a trajetória dos primeiros e dos últimos 30,0 s podem ser vistos nas figu-ras 5.11(a) e 5.11(b) respectivamente de um tempo total de 300,0 s. Na figura 5.11(b) observa-seque, aparentemente, a trajetória está indo para um ciclo limite.
00.5
11.5
0
0.2
0.4−1
−0.5
0
0.5
1
x2 [rad/s]x
3 [rad]
x 4 [rad
/s]
(a) Primeiros 30,0 s.
−3−2
−10
12
34
x 10−30.61
0.62
0.63
0.64
0.65
−0.4
−0.2
0
0.2
0.4
X [rad/s]
Y [rad]
Z [r
ad/s
]
(b) Últimos 30,0 s de um total de 300,0 s.
Figura 5.10: Espaço de fases das variáveis x2, x3 e x4 com bx = 2,2× 10−4 N ·m · s/rad.
01
23
0
0.5
1−1
−0.5
0
0.5
1
x2 [rad/s]x
3 [rad]
x 4 [rad
/s]
(a) Primeiros 30,0 s.
00.2
0.40.6
0.81
x 10−5
0.55
0.6
0.65
0.7−1
−0.5
0
0.5
1
X [rad/s]Y [rad]
Z [r
ad/s
]
(b) Últimos 30,0 s de um total de 300,0 s.
Figura 5.11: Espaço de fases das variáveis x2, x3 e x4 com bx = 7,3× 10−4 N ·m · s/rad.
80
5.2.2 Expoentes de Lyapunov experimental
Para a obtenção dos expoentes de Lyapunov, baseado em dados experimentais, é neces-sário ter uma série temporal de uma das grandezas físicas do giroscópio. Como visto anterior-mente, as grandezas do movimento do giroscópio que puderam ser medidas foram as posiçõesangulares de precessão e nutação. Para a análise dos expoentes de Lyapunov, através de sériestemporais, o parâmetro usado foi a posição do contrapeso m1, conforme descrição da seção 5.2.
Para a série temporal da posição angular de precessão, com o contrapeso m1 nas posições−0,1079 m e −0,1179 m, foram obtidos os seguintes dados: tempo de atraso 248 e 163; dimen-são de imersão 1 e 2, respectivamente. Já para a série temporal da posição angular de nutação,os dados foram: tempo de atraso 168 e 155; dimensão de imersão 5 e 6, respectivamente. Comoa dimensão de imersão relacionada a posição de precessão, encontrada para as duas posiçõesde m1, foram muito baixas em comparação com o valor encontrado para a posição de nutação,conclui-se que a série temporal da posição angular de precessão não representa bem o sistema,em comparação com a série temporal da posição angular de nutação. Por isso a análise seráfeita em relação a esta última série.
Na tabela 5.4 pode-se ver um resumo das informações obtidas, da série temporal da posi-ção angular de nutação. O tempo de amostragem foi de ts = 0,001 s e o passo de 250.
Tabela 5.4: Resumo da série temporal experimental da posição angular de nutação x3. Natabela τ representa tempo de atraso, m representa dimensão de imersão, λl e λw representamos expoente de Lyapunov, usando os programas lyapmax e wolf respectivamente.
Posição pontos τ m λl λwde h1 (m) (1/s)
−0,1079 49793 168 5 −7,70×10−2 8,33×10−3
−0,1179 33655 155 6 7,19×10−2 1,17×10−1
−0,1279 34141 136 10 1,06×10−1 −2,40×10−2
−0,1379 31282 139 10 1,12×10−1 3,09×10−2
−0,1479 30611 159 7 6,84×10−2 1,27×10−1
−0,1579 31674 139 5 9,79×10−2 −1,65×10−1
−0,1679 31136 158 7 8,69×10−2 1,01×10−1
−0,1779 27546 137 5 6,62×10−2 1,13×10−1
−0,1879 32423 133 5 1,33×10−1 1,15×10−1
−0,1979 34329 163 6 2,07×10−1 4,72×10−2
−0,2079 54243 245 10 1,97×10−2 −2,27×10−1
−0,2179 65664 327 5 1,11×10−1 −8,87×10−2
−0,2279 62608 201 10 −2,40×10−3 4,42×10−2
−0,2379 71556 286 6 9,79×10−2 8,90×10−2
−0,2479 65126 304 6 −2,38×10−2 7,10×10−2
−0,2579 117889 433 10 −2,19×10−2 4,31×10−3
Para transformar número de pontos n no tempo total tt, basta usar a expressão: tt =
81
(n − 1)ts, logo o tempo máximo obtido, foi de ts = 117,888 s. Observou-se que a dimensãode imersão varia de 5 a 10. Exemplo de convergência do expoente de Lyapunov, se encontranas figuras 5.12(a) e 5.12(b), para o caso em que o contrapeso m1 esteja na posição h1 =
−0,1979 m.
0 5 10 15 20 25 30 350
0.2
0.4
0.6
0.8
1
1.2
1.4
1.6
1.8
2
Tempo [s]
Con
verg
ênci
a do
Mai
or E
xpoe
nte
de L
yapu
nov
[1 /
s]
(a) Usando o programa lyapmax.
0 5 10 15 20 25 30 35−0.8
−0.6
−0.4
−0.2
0
0.2
0.4
0.6
0.8
Tempo [s]
Con
verg
ênci
a do
Mai
or E
xpoe
nte
de L
yapu
nov
[1 /
s]
(b) Usando o programa wolf.
Figura 5.12: Exemplo de convergência de λ da série temporal experimental de x3, com h1 =−0,1979 m.
Com base na tabela 5.4, obtém-se a figura 5.13. Os resultados não concordam numerica-mente, apesar dos programas lyapmax e wolf serem, a principio, baseados no mesmo algoritmo;em relação aos sinais dos expoentes, os dois programas concordam para certos parâmetros,sendo o sinal que apresenta concordância positivo.
−0.26 −0.24 −0.22 −0.2 −0.18 −0.16 −0.14 −0.12 −0.1−0.25
−0.2
−0.15
−0.1
−0.05
0
0.05
0.1
0.15
0.2
0.25
Parâmetro h1 [m]
Exp
oent
es d
e Ly
apun
ov [1
/ s]
(lyapmax)(wolf)
(wolf)
(lyapmax)
Figura 5.13: Expoentes de Lyapunov da série temporal experimental de x3.
A fim de verificar se o modelo e os resultados experimentais do giroscópio, apresentam omesmo comportamento não-linear, foi feita uma comparação. Para isso simulou-se uma sérietemporal resolvendo a equação diferencial do modelo nas mesmas condições do experimento e
82
usando a posição angular de nutação como uma série temporal. Na tabela 5.5 apresenta-se umresumo das propriedades da série simulada.
Tabela 5.5: Resumo da série temporal simulada da posição angular de nutação x3. Na tabela τrepresenta tempo de atraso, m representa dimensão de imersão, λl e λw representam os expo-ente de Lyapunov, usando os programas lyapmax e wolf respectivamente.
Posição pontos τ m λl λwde h1 (m) (1/s)
−0,1079
60001
294 2 1,01×10−2 −1,19×10−4
−0,1179 289 3 −6,95×10−3 −1,10×10−3
−0,1279 277 3 −1,53×10−2 −1,09×10−3
−0,1379 267 3 −1,08×10−1 8,74×10−4
−0,1479 273 3 −1,69×10−1 −7,18×10−4
−0,1579 257 3 −1,73×10−1 −6,13×10−2
−0,1679 244 8 −2,14×10−2 −9,30×10−2
−0,1779 255 5 −1,17×10−1 −1,00×10−1
−0,1879 222 3 −1,84×10−1 −3,38×10−2
−0,1979 231 3 −3,06×10−1 −3,29×10−2
−0,2079 228 3 −2,62×10−1 −4,22×10−2
−0,2179 215 3 −1,84×10−1 −5,49×10−2
−0,2279 203 3 −1,80×10−1 −7,47×10−2
−0,2379 196 3 −1,83×10−1 −4,46×10−2
−0,2479 195 3 −1,93×10−1 −4,39×10−2
−0,2579 186 3 −1,95×10−1 −4,59×10−2
A dimensão de imersão da série temporal simulada apresenta, na média, um valor maisbaixo comparado com a serie temporal experimental. Nas figuras 5.14(a) e 5.14(b), pode servisto a convergência do expoente de Lyapunov da série simulada, para as mesmas condições dasérie experimental.
0 10 20 30 40 50 60−0.4
−0.3
−0.2
−0.1
0
0.1
0.2
0.3
Tempo [s]
Con
verg
ênci
a do
Mai
or E
xpoe
nte
de L
yapu
nov
[1 /
s]
(a) Usando o programa lyapmax.
0 10 20 30 40 50 60−0.1
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
Tempo [s]
Con
verg
ênci
a do
Mai
or E
xpoe
nte
de L
yapu
nov
[1 /
s]
(b) Usando o programa wolf.
Figura 5.14: Exemplo de convergência de λ da série temporal simulada de x3, com h1 =−0,1979 m.
83
Para ver de forma mais clara os resultados dos expoentes de Lyapunov da série tempo-ral simulada da tabela 5.5, foi criada a figura 5.15. Pode-se verificar nesta figura que não háconcordância numérica entre os programas usados. Porém, em relação ao sinal do expoente, háuma concordância em quase todos os parâmetros, tendo o sinal negativo, ao contrário do queocorreu para a serie temporal experimental, que para a maioria dos parâmetros apresentou sinalpositivo.
−0.26 −0.24 −0.22 −0.2 −0.18 −0.16 −0.14 −0.12 −0.1−0.35
−0.3
−0.25
−0.2
−0.15
−0.1
−0.05
0
0.05
Parâmetro h1 [m]
Exp
oent
es d
e Ly
apun
ov [1
/ s]
(lyapmax)(wolf)
(lyapmax)
(wolf)
Figura 5.15: Expoentes de Lyapunov da série temporal simulada de x3.
Nos resultados dos expoentes de Lyapunov baseados nos programas lyapmax e wolfhouve diferenças. No caso os testes feitos no apêndice C tiveram melhores resultados como programa lyapmax apesar de ambos os programas serem baseados no mesmo algoritmo. Oprograma lyapmax pode ter alguma melhoria em relação ao algoritmo original. Isso poderia serconfirmado conseguindo o código fonte do programa lyapmax com os autores do mesmo.
Pode-se observar que os expoentes de Lyapunov experimentais oscilam muito com a va-riação da posição h1. Como os expoentes de Lyapunov do modelo (figuras 5.7(a) e 5.7(b))nas mesmas condições as variações do maior expoente não foram bruscas e as amplitudes dosexpoente menores, provavelmente devido a ser usado um tempo bem maior.
Os experimentos não puderam ter um tempo mais longo isso acaba afetando e muito osresultados. Por isso os resultados não foram muito conclusivos em relação aos expoentes deLyapunov aplicados a séries temporais experimentais.
5.3 Teste 0− 1
No teste 0− 1 utilizou-se os dados de séries temporais, assim como o tempo de atraso. Avantagem do teste 0− 1 é que não é necessária a reconstrução do atrator do sistema. Neste mé-todo o tempo de atraso é utilizado para se evitar o problema de sobreamostragem. O programa
84
utilizado encontra-se no endereço eletrônico (Matthews, 2009).
Na tabela 5.6 apresentam-se os dados relevantes a respeito da série temporal experimentalda posição angular de nutação, assim como os valores de K. O tempo de amostragem originalda série foi de ts = 0,001 s e a série foi re-amostrada conforme o tempo de atraso empregado.Para a série temporal simulada os dados se encontram na tabela 5.7.
Tabela 5.6: Resumo da série temporal experimental da posição angular de nutação x3. Na tabelaτ representa tempo de atraso, K é o resultado do teste 0− 1.
Posição de h1 (m) pontos τ K
−0,1079 49793 168 −0,118−0,1179 33655 155 −0,059−0,1279 34141 136 −0,023−0,1379 31282 139 −0,049−0,1479 30611 159 −0,087−0,1579 31674 139 −0,141−0,1679 31136 158 −0,074−0,1779 27546 137 −0,054−0,1879 32423 133 −0,105−0,1979 34329 163 −0,058−0,2079 54243 245 −0,183−0,2179 65664 327 −0,166−0,2279 62608 201 −0,186−0,2379 71556 286 −0,135−0,2479 65126 304 −0,142−0,2579 117889 433 −0,134
Conforme as tabelas 5.6 e 5.7 pode-se ver que os valore de K estão mais próximos dezero, portanto pelo teste 0 − 1, em ambos os casos o comportamento do sistema dinâmico éregular. Porém as séries temporais são curtas, podendo os valores de K não terem convergido.Para tirar essa dúvida, obtém-se o gráfico p− q para uma das séries simuladas e uma das sériesexperimentais. Adota-se o parâmetro h1 = −0,1979 m que corresponde ao meio termo, ou seja,não se encontra nas extremidades da haste. Na figura 5.16(a) mostra-se o caso simulado, deacordo com o que foi visto no final da seção 4.3 o sistema apresenta dinâmica regular. Para asérie experimental tem-se a figura 5.16(b). Não está claro, pela figura, se a dinâmica é regularou apresenta caos fraco, pois há poucos pontos na mesma.
5.4 Conclusão
No sistema giroscópico em estudo foi verificado que existem infinitos pontos de equilíbrioe da análise da estabilidade dos mesmos constatou-se que eles eram elípticos, pois a matrizJacobiana do sistema apresentou autovalores nulos e por isso a análise do sistema linear através
85
Tabela 5.7: Resumo da série temporal simulada da posição angular de nutação x3. Na tabela τrepresenta tempo de atraso, K é o resultado do teste 0− 1.
Posição de h1 (m) pontos t K
−0,1079
60001
294 −0,142−0,1179 289 −0,151−0,1279 277 −0,099−0,1379 267 −0,115−0,1479 273 −0,109−0,1579 257 −0,117−0,1679 244 −0,122−0,1779 255 −0,101−0,1879 222 −0,083−0,1979 231 −0,063−0,2079 228 −0,073−0,2179 215 −0,065−0,2279 203 −0,050−0,2379 196 −0,053−0,2479 195 −0,041−0,2579 186 −0,022
−1 −0.5 0 0.5 1−1
−0.8
−0.6
−0.4
−0.2
0
0.2
0.4
0.6
0.8
1
p
q
(a) Caso simulado.
−0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6−0.6
−0.4
−0.2
0
0.2
0.4
0.6
0.8
p
q
(b) Caso experimental.
Figura 5.16: Gráfico p− q do teste 0− 1 com h1 = −0,1979 m.
dos autovalores não foi conclusiva.
Para o estudo dos expoentes de Lyapunov em relação ao amortecimento viscoso, verificou-se que as somas dos expoentes de Lyapunov na média, independente do método aplicado, fica-vam mais negativas conforme aumentavam os parâmetros de amortecimento, pois o sistema édissipativo.
Apesar dos expoentes de Lyapunov para os métodos padrão e de Eckmann-Ruelle nãoconcordarem na análise do amortecimento viscoso a soma dos expoentes difere pouco, indi-cando que as orientações dos expoentes é que não coincidem.
Da análise feita nos parâmetros relacionados ao amortecimento viscoso e a posição do
86
contrapeso maior utilizando os métodos padrão e de Eckmann-Ruelle, observou-se que a ordemdos expoentes em relação a suas amplitudes se alterna. Também é observado que os expoentesencontrados para o sistema giroscópico, embora sejam positivos, tem um amplitude baixa comum valor máximo da ordem de 0,035 para o método de Eckmann-Ruelle. Com um valor baixopara os expoentes de Lyapunov mais o fato de a ordem dos expoentes ficar alternando conclui-se que esses expoentes estão convergindo para zero, isso foi confirmado pelo teste 0− 1, sendoassim o comportamento do sistema é regular.
Na análise experimental do giroscópio utilizou-se para a série temporal, dados da posiçãoangular de nutação. Para os expoentes de Lyapunov foi utilizado como parâmetro a posição docontrapeso maior m1.
Para as séries temporais simuladas e experimentais não houve concordância entre os ex-poentes de Lyapunov obtidos com os programas wolf e lyapmax. Com os dados simuladosa maior parte dos expoentes foram negativos, já com os dados experimentais a maioria dosexpoentes teve sinal positivo.
Observou-se que os expoentes de Lyapunov tanto simulados quanto experimentais osci-laram muito com a variação da posição h1. Como os expoentes de Lyapunov do modelo nasmesmas condições as variações do maior expoente não foram brusca e as amplitudes dos expo-ente foram menores, provavelmente devido ao tempo utilizado ser bem maior.
Os experimentos não puderam ter um tempo mais longo, por isso os resultados não forammuito conclusivos em relação aos expoentes de Lyapunov aplicados as séries temporais tantosimuladas quanto experimentais.
O teste 0− 1 foi aplicado nas mesmas séries temporais utilizadas para obter os expoentesde Lyapunov experimentais, no teste os valores de K ficaram próximos de zero, portanto peloteste 0 − 1 o comportamento do sistema dinâmico deve ser regular, porém as séries temporaissão curtas, podendo os valores de K não terem convergido. Visualizando os gráficos p − q foiconcluído que no caso simulado o sistema tem comportamento regular, já o caso experimentalnão foi possível dizer se a dinâmica é regular, pois não havia dados o suficiente.
Capítulo 6
Conclusões Finais e Trabalhos Futuros
Neste trabalho foi feito a modelagem do giroscópio da Pascor e foi verificado a validadedo mesmo com dados experimentais. Com o modelo matemático foram obtidos os pontos deequilíbrio e também foi verificada a estabilidade do sistema giroscópico nestes pontos. Porfim, os expoentes de Lyapunov foram obtidos tanto do modelo matemático quando dos dadosexperimentais.
O giroscópio da Pascor foi modelado utilizando conceitos da mecânica de Newton-Euler,sendo levado em consideração efeitos elásticos e dissipativos. As equações diferenciais obtidasdescrevem o movimento de precessão, nutação e spin do giroscópio. As equações diferenciaisforam resolvidas numericamente com o método de Runge-Kutta de oitava ordem com passovariável e as soluções foram comparadas com os dados experimentais das posições angularesde precessão e de nutação.
Em relação a obtenção dos dados experimentais observou-se que o cabo do sensor daposição angular de nutação induz um amortecimento que dificulta uma aquisição de dados maislonga e confiável. Experimentalmente só foi possível obter os dados das posições angulares deprecessão e de nutação e a velocidade de spin foi estimada para a condição inicial com a ajudade um tacômetro.
Os resultados mostraram que o modelo proposto apresenta uma boa concordância com omovimento real do giroscópio e pode ser utilizado para as finalidades de análise não-linear.
Em relação ao estudo da estabilidade dos pontos de equilíbrio não foi possível verificar,através da análise do sistema linearizado, se os mesmos são estáveis ou instáveis pois os pontosde equilíbrio são elípticos.
Já no estudo dos expoentes de Lyapunov, verificou-se que as somas dos expoentes sãoinversamente proporcionais aos parâmetros do amortecimento viscoso. Também foi verificadoque os expoentes de Lyapunov, em relação ao amortecimento viscoso, obtidos pelos métodospadrão e de Eckmann-Ruelle não apresentaram concordância, mas a soma dos expoentes dife-riram pouco, indicando que as orientações dos expoentes não eram as mesmas. Na análise dosexpoentes de Lyapunov do modelo, observou-se que a ordem dos expoentes em relação a suas
87
88
amplitudes se alternavam e que embora os expoentes sejam positivos, estes tem uma amplitudebaixa com um valor máximo da ordem de 0,035 para o método de Eckmann-Ruelle. Com umvalor baixo para os expoentes de Lyapunov mais o fato de a ordem dos expoentes ficar alter-nando conclui-se que esses expoentes estão convergindo para zero, isso foi confirmado peloteste 0− 1, sendo assim o comportamento do sistema é regular.
Na análise experimental do giroscópio utilizou-se, para a série temporal, dados da posiçãoangular de nutação. Para os expoentes de Lyapunov foi utilizado como parâmetro a posição docontrapeso maior m1.
Para as séries temporais simuladas e experimentais não houve concordância entre os ex-poentes de Lyapunov obtidos com os programas wolf e lyapmax. Com os dados simuladosa maior parte dos expoentes foram negativos, já com os dados experimentais a maioria dosexpoentes teve sinal positivo.
Observou-se que os expoentes de Lyapunov tanto simulados quanto experimentais osci-laram muito com a variação da posição h1 do contrapeso.
Os experimentos não puderam ter um tempo mais longo, por isso os resultados não forammuito conclusivos em relação aos expoentes de Lyapunov aplicados as séries temporais tantosimuladas quanto experimentais.
O teste 0− 1 foi aplicado nas mesmas séries temporais utilizadas para obter os expoentesde Lyapunov experimentais e simulados. No teste os valores de K ficaram próximos de zero,portanto pelo teste 0−1 o comportamento do sistema dinâmico deve ser regular, porém as sériestemporais são curtas, podendo os valores de K não terem convergido. Visualizando os gráficosp − q foi concluído que no caso simulado o sistema tem comportamento regular, já o casoexperimental não foi possível dizer se a dinâmica é regular, pois não haviam dados suficientes.
Com base no que foi visto neste trabalho e seguindo a mesma linha de pesquisa, pode-sefazer sugestões para trabalhos futuros em dois aspectos:
Aspecto experimental: Utilizar um giroscópio mais preciso, com mais sensores, melhor placade aquisição e métodos para obter condições iniciais quaisquer;
Aspecto de análise não-linear: Incluir o estudo de bifurcações, mapas de Poincaré, outros ti-pos de expoentes de Lyapunov.
Referências Bibliográficas
Aballay, E. and Avilés, E. (2002). Estudio del Movimiento Giroscópico.
Acar, C. and Shkel, A. (2008). MEMS Vibratory Gyroscopes: Structural Approaches to ImproveRobustness, Mems Reference Shelf, Springer.
Addison, P. S. (1997). Fractals and Chaos: An Illustrated Course, Institute of Physics Pub.
Angular (2012). Foto: precessão, nutação e spin. Último acesso em 06/08/2012.Disponível em: http://en.wikipedia.org/wiki/File:Praezession.svg
Arnold, V. I. (1987). Métodos Matemáticos da Mecânica Clássica, Editora Mir Moscovo.
Bryant, P. H. (2009). Método de Eckmann-Ruelle. Código de programa: cálculo de expoentede Lyapunov. Último acesso em 06/08/2012.Disponível em: http://inls.ucsd.edu/~pbryant/
Campanharo, A. S. L. O., Macau, E. E. N. and Ramos, F. M. (2005). Detectando a Presençade Caos em uma Série Temporal, XXVIII CNMAC - Congresso Nacional de MatemáticaAplicada e Computacional, São Paulo.
Carrera, D. H. Z. and Weber, H. I. (2009). Dynamics of Bodies in Space Rotating across Stabi-lity Borders, 8th Brazilian Conference on Dynamics, Control and Applications, DINCON.Baurú. Proceedings of DINCON 2009. Rio Claro : SBMAC.
Carrera, D. H. Z., Weber, H. I. and Morrot, R. (2010). Alternatives of Long Term Behavior of aGyroscope Comsidering Damping, PACAM XI- 11th Pan-American Congress of AppliedMechanics. Foz do Iguaçu. Proceedings of PACAM 2010. São Carlos : USP.
Chen, H.-K. and Ge, Z.-M. (2005). Bifurcations and chaos of a two-degree-of-freedom dissi-pative gyroscope, Chaos, Solitons & Fractals 24(1): 125 – 136.
Crabtree, H. (1909). An elementary treatment of the theory of spinning tops and gyroscopicmotion, Longmans, Green and Co. in London, New York . Editor: Universidade de Mi-chigan.
Eckmann, J. P. and Ruelle, D. (1985). Ergodic theory of chaos and strange attractors, Rev. Mod.Phys. 57: 617–656.
Falconer, I., Gottwald, G. A., Melbourne, I. and Wormnes, K. (2007). Application of the 0-1 Test for Chaos to Experimental Data, SIAM Journal on Applied Dynamical Systems6(2): 395–402.
89
90
Fazanaro, F. I., Soriano, D. C., Madrid, M. K., Suyama, R., Attux, R. and de Oliveira, J. R.(2010). Cálculo do Espectro de Lyapunov via Dinâmicas Clonadas e sua Aplicação noCircuito de Chua, XVIII Congresso Brasileiro de Automática (CBA), pp. 2511–2518. Bo-nito, MS, Brasil.
Fiedler-Ferrara, N. and do Prado, C. P. C. (1994). Caos - Uma Introdução, 1 edn, EditoraEdgard Blücher LTDA.
Fraser, A. M. and Swinney, H. L. (1986). Independent coordinates for strange attractors frommutual information, Physical Review A 33: 1134–1140.
Ge, Z. M., Chen, H. K. and Chen, H. H. (1996). The Regular and Chaotic Motions of aSymmetric Heavy Gyroscope with Harmonic Excitation, Journal of Sound and Vibration198(2): 131 – 147.
Goldstein, H. (1981). Classical Mechanics, 2 edn, Addison-Wesley Publishing Company.
Gottwald, G. A. and Melbourne, I. (2005). Testing for chaos in deterministic systems withnoise, Physica D: Nonlinear Phenomena 212(1-2): 100–110.
Gottwald, G. A. and Melbourne, I. (2009a). On the Implementation of the 0-1 Test for Chaos,SIAM Journal on Applied Dynamical Systems 8(1): 129.
Gottwald, G. A. and Melbourne, I. (2009b). On the Validity of the 0-1 Test for Chaos, Nonline-arity 22(6): 1–21.
Gottwald, G. and Melbourne, I. (2004). A new test for chaos in deterministic systems, Procee-dings of the Royal Society of London A 460: 603–611.
Hairer, E., Nørsett, S. P. and Wanner, G. (1993). Solving Ordinary Differential Equations I:Nonstiff Problems, 2 edn, Springer-Verlag.
Hairer, E. and Wanner, G. (1991). Solving Ordinary Differential Equations II: Stiff andDifferential-Algebraic Problems, Springer-Verlag.
Kantz, H. and Schreiber, T. (2004). Nonlinear Time Series Analysis, Cambridge nonlinearscience series, Cambridge University Press.
Ke-Hui, S., Xuan, L. and Cong-Xu, Z. (2010). The 0-1 test algorithm for chaos and its applica-tions, Chinese Physics B 19(11): 110513.
Kodba, S., Perc, M. and Marhl, M. (2005). Detecting chaos from a time series, EuropeanJournal of Physics 26(1): 205–215.
Kostov, S. and Hammer, D. (2010). “It Has to Go Down A Little, In Order to Go Around” -Following Feynman on the Gyroscope.
Lawrence, A. (1998). Modern Inertial Technology: Navigation, Guidance, and Control, Me-chanical Engineering Series, Springer.
Lemos, N. A. (2007). Mecânica Analítica, 2 edn, Editora Livraria da Física.
Lorenz, E. N. (1963). Deterministic Nonperiodic Flow, Journal of the Atmospheric Sciences20(2): 130–141.
91
Marion, J. B. and Thornton, S. T. (1995). Classical Dynamics of Particles and Systems, 4 edn,Harcourt College Publishers.
Matthews, P. (2009). Teste 0-1. Código de programa: teste 0-1. Último acesso em 06/08/2012.Disponível em: http://www.mathworks.com/matlabcentral/fileexchange/25050-0-1-test-for-chaos/all_files
McGlynn, E. (2007). Introducing Gyroscopes Quantitatively without Puting Students into aSpin, European Journal of Physics 28(3): 479–486.
Monteiro, L. H. A. (2006). Sistemas Dinâmicos, 2 edn, Editora Livraria da Física.
Mörters, P., Peres, Y., Schramm, O. and Werner, W. (2010). Brownian Motion, CambridgeSeries in Statistical and Probabilistic Mathematics, Cambridge University Press.
Neto, J. B. (2004). Mecânica Newtoniana, Lagrangeana e Hamiltoniana, 1 edn, Editora Livra-ria da Física.
PASCO (1996). Rotary Motion Sensor for ULI – Model CI-6625, PASCO Scientific.
PASCO (1997). DataStudio from PASCO Scientific Model, CI-6859C, PASCO Scientific.
Perc, M. (2005a). Nonlinear time series analysis of the human electrocardiogram, EuropeanJournal of Physics 26(5): 757–768.
Perc, M. (2005b). The dynamics of human gait, European Journal of Physics 26(3): 525–534.
Perc, M. (2006). Introducing Nonlinear Time Series Analysis in Undergraduate Courses, FI-ZIKA A 15(2): 91–112.
Perc, M. (2007). Reconstrução do atrator, expoente de Lyapunov. Programas: maior expo-ente de Lyapunov pelo método de Wolf, tempo de atraso pelo método de informação mú-tua, dimensão de imersão pelo método dos falsos vizinhos próximos. Último acesso em06/08/2012.Disponível em: http://www.matjazperc.com/ejp/time.html
Ramasubramanian, K. and Sriram, M. S. (2000). A comparative study of computation of Lya-punov spectra with different algorithms, Physica D: Nonlinear Phenomena 139(1-2): 72– 86.
Santos, I. F. (2001). Dinâmica de Sistemas Mecânicos - Modelagem, Simulação, Visualização,Verificação, 1 edn, Makron Books LTDA.
Savi, M. A. (2006). Dinâmica Não-linear e Caos, 1 edn, Editora E-papers.
Small, M. (2005). Applied Nonlinear Time Series Analysis: Applications in Physics, Physiologyand Finance, World Scientific Series on Nonlinear Science: Monographs and Treatises,World Scientific.
Strogatz, S. H. (1994). Nonlinear Dynamics And Chaos: With Applications To Physics, Biology,Chemistry, And Engineering, Studies in Nonlinearity, Westview Press.
Symon, K. R. (1996). Mecânica, 3 edn, Editora Campus LTDA.
92
Takwale, R. G. and Puranik, P. S. (1979). Introduction to Classical Mechanics, TATAMacGraw-Hill Publishing Company Ltd.
Titterton, D. H. and Weston, J. L. (2004). Strapdown Inertial Navigation Technology, Iee RadarSeries, 2 edn, Institution of Electrical Engineers.
Vuolo, J. H. (1996). Fundamentos da Teoria de Erros, 2 edn, Editora Edgard Blücher LTDA.
Wiggins, S. (2003). Introduction to Applied Nonlinear Dynamical Systems and Chaos, Texts inApplied Mathematics, Springer.
Wolf, A., Swift, J. B., Swinney, H. L. and Vastano, J. A. (1985). Determining Lyapunov expo-nents from a time series, Physica D: Nonlinear Phenomena 16(3): 285 – 317.
Apêndice A
Obtenção dos Parâmetros Indiretos doModelo
Na equação diferencial do giroscópio obtida no capítulo 3 pode-se observar a existênciade muitos parâmetros, muitos destes podem ser medidos diretamente como as posições doscentros de massa e as massas do contra peso grande (h1,m1), contra peso pequeno (h2,m2),disco (h,m) e haste (hh,mh). O valor de mh é obtido indiretamente, outras massas como ado sensor de posição angular de nutação são adicionadas a mh isto é feito por simplificação, ocálculo de mh consiste em por o giroscópio em equilíbrio (H = 0) e, através da expressão docentro de massa tem-se:
mh =−m1h1 −m2h2 −mdhd
hh. (A.1)
Na equação diferencial do giroscópio existem vários outros parâmetros que não podemser medidos diretamente como a aceleração da gravidade g, o parâmetro do torque elástico k,os parâmetros do amortecimento viscoso (bx, by e bz) e os parâmetros do torque cinético (Tx, Tye Tz). A seguir será descrito como foi obtido cada um destes parâmetros.
Uma maneira de obter os parâmetros de uma EDO é através do método dos mínimosquadrados, mas para esse método é necessário conhecer a forma geral das soluções da EDO.Como o sistema é não-linear não há possibilidade de obter as soluções analíticas e, portanto,não dá pra aplicar o método dos mínimos quadrados.
Apesar de não haver possibilidade de aplicar o método dos mínimos quadrados, é possívelusar um conceito do qual ele se baseia, tal conceito chama-se conceito de máxima verossimi-lhança (Vuolo, 1996). O conceito de máxima verossimilhança consiste em minimizar as dis-tâncias entre os resultados experimentais e os resultados da EDO, mais precisamente encontrarum conjunto de parâmetros µ que minimize S2, tem-se:
S2 =N−1∑i=0
[xi − x(ti,µ)]2 (A.2)
93
94
eµ =
µ1 µ2 . . . µk
T(A.3)
sendo xi os dados experimentais e x(ti,µ) a solução da EDO que no caso será numérica.
A variável S2 será minimizada utilizando força bruta. De forma resumida, a partir doprimeiro parâmetro µ1 escolhe-se uma faixa para o mesmo (µ11, . . ., µ1i, . . .) e para cada µ1i
resolve-se a EDO com um conjunto de parâmetros iniciais e verifica-se para qual dos parâmetrosµ1i o valor de S2 é minimizado, esse parâmetro é guardado. Repete-se o processo para os outrosparâmetros, aos poucos a variável S2 vai sendo minimizada. O tempo de processamento é umpouco alto, mas a implementação é relativamente simples. Existem, é claro, outros métodosmais sofisticados como o algoritmo genético, mas optou-se por um mais simples que pudesseser aplicado de forma imediata.
Além dos parâmetros indiretos já citados, as próprias condições iniciais são consideradascomo parâmetros, já que há incertezas nas mesmas. Os dados experimentais usados foram daposição angular de precessão e de nutação com o contrapeso maior na posição h1 = −0,1979 m
e o menor na posição h2 = −0,0827 m, as condições iniciais foram as descritas na seção 3.3 eo tempo total do experimento foi de t = 34,32 s.
Para resolver a EDO, o tempo de amostragem da solução foi t = 0,01 s, este foi escolhidode forma que o tempo de processamento não fosse muito elevado. O erro adotado para a soluçãonumérica foi 1,0 × 10−11 para o erro relativo e 1,0 × 10−7 para o erro absoluto de forma queo método de Runge-Kutta não apresentasse problemas, os parâmetros diretos que serão usadospara resolver a EDO são os da seção 3.3.
Nas figuras A.1(a) e A.1(b) podem ser vistas as comparações dos dados experimentais daprecessão com os obtidos pela solução numérica da EDO usando as condições iniciais e os pa-râmetros ajustados pelo método de força bruta. Da figura A.1(b) percebe-se que as frequênciasdos dados experimentais e simulados diferem. Poderia-se tentar implementar algum mecanismoque leve as frequências dos sinais em consideração com o intuito de melhorar a otimização.
No movimento de nutação, as comparações dos dados experimentais com os simuladospodem ser vistas nas figuras A.2(a) e A.2(b). Existe uma oscilação extra nos dados experi-mentais que modula a oscilação principal, tal oscilação extra não foi prevista pelo modelo, ainterferência não parece ser causada pelo fio do sensor de nutação, pois este deveria interferirno movimento de precessão, é mais provável que alguma hipótese não tenha sido consideradadurante a modelagem como por exemplo considerar a base onde se localiza o giroscópio comonão inercial. Também foi verificado a necessidade de aplicar algum mecanismo que faça comque as frequências dos movimentos experimental e simulados fiquem mais próximas.
O valor obtido para a aceleração da gravidade g pelo ajuste de força bruta foi g =
14,45 m/s2, um valor acima do esperado. Erros nos valores das posições e das massas pro-pagam erros no valor do centro de massa H que pode ser uma das causas do valor de g acima
95
0 5 10 15 20 25 30 350
5
10
15
20
25
30
35
40
tempo [s]
x 1 [ra
d]
SimuladoExperimental
(a) Comparação usando o tempo total.
0 1 2 3 4 5 6 7 80
1
2
3
4
5
6
7
8
9
tempo [s]
x 1 [ra
d]
SimuladoExperimental
(b) Comparação usando um quarto do tempo total.
Figura A.1: Comparação dos dados experimentais e simulados da precessão após ajuste.
0 5 10 15 20 25 30 35−0.1
0
0.1
0.2
0.3
0.4
0.5
0.6
tempo [s]
x 3 [rad
]
SimuladoExperimental
(a) Comparação usando o tempo total.
0 1 2 3 4 5 6 7 8−0.05
0
0.05
0.1
0.15
0.2
0.25
0.3
tempo [s]
x 3 [rad
]
SimuladoExperimental
(b) Comparação usando um quarto do tempo total.
Figura A.2: Comparação dos dados experimentais e simulados da nutação após ajuste.
do esperado. Daqui por diante será usado o valor adotado no Laboratório de Física da Unioestecujo valor é g = 9,78 m/s2 que foi obtido em experimentos com pendulo simples.
Os parâmetros indiretos obtidos nesse apêndice encontra-se na tabela A.1 e as condiçõesiniciais encontradas são:
X0 =
0,000 0,000 0,000 0,519 0,000 60,777T
(A.4)
Uma estimativa para o “erro” da simulação em relação aos dados experimentais é dado
96
Tabela A.1: Parâmetros indiretos do giroscópio da Pascor.T ( N ·m) b( N ·m · s/rad) g( m/s2) k( N/rad)
x 2,933× 10−03 5,604× 10−05
y 1,366× 10−03 4,652× 10−14
z 5,360× 10−03 4,645× 10−03
9,78 16,81
por:
S =
√∑N−1i=0
[xpi − xp(ti,µ)]2 + [xni − xn(ti,µ)]2
N
(A.5)
sendo xpi e xni os dados experimentais das posições angulares de precessão e nutação respec-tivamente. A estimativa do “erro” pela eq. (A.5) foi S = 0,9080, isto da uma ideia de quãopróximos estão os dados experimentais dos simulados.
Apêndice B
Jacobiana da Equação Diferencial doGiroscópio
A matriz jacobiana J é utilizada para linearizar uma equação diferencial não-linear emtorno de um ponto. Nesta dissertação, a matriz jacobiana é usada para determinar os expoen-tes de Lyapunov do modelo do giroscópio assim como verificar a estabilidade dos pontos deequilíbrio.
Para um sistema dinâmico não-linear cuja equação diferencial tem a forma:
X = F (X); X ∈ <n (B.1)
a matriz jacobiana é calculada como segue:
J = [aij]n×n; aij =∂fi∂xj
(B.2)
A equação diferencial do giroscópio foi obtida no capítulo 3, para facilitar os cálculos amesma é reescrita abaixo:
x1
x2
x3
x4
x5
x6
=
f1(x1, x2, x3, x4, x5, x6)
f2(x1, x2, x3, x4, x5, x6)
f3(x1, x2, x3, x4, x5, x6)
f4(x1, x2, x3, x4, x5, x6)
f5(x1, x2, x3, x4, x5, x6)
f6(x1, x2, x3, x4, x5, x6)
(B.3)
sendo as equações diferenciais no espaço de estado:
x1 = x2 (B.4)
97
98
x2 =τz + x4[Ixxx6 + (−Ix + Iy + Iz)x2 sin(x3)]
Iz cos(x3)(B.5)
x3 = x4 (B.6)
x4 =τy − Ixxx2x6 cos(x3) + (Ix − Iz)x22 cos(x3) sin(x3)
Iy(B.7)
x5 = x6 (B.8)
x6 = 1Ixxτx + (Ix + Iy − Iz)x2x4 cos(x3) + Ix
Izτz + x4[Ixxx6+
(−Ix + Iy + Iz)x2 sin(x3)] tan(x3)(B.9)
e as componentes do torque resultante são:τx
τy
τz
=
−Tx sgn (x6) + [Tz sgn (x2) + bxx2] sin(x3)− bxx6MgH cos(x3)− Ty sgn (x4)− byx4 + τey (x3)
−[Tz sgn (x2) + bzx2] cos(x3)
(B.10)
a componente na direção 2 do torque elástico que está na eq. (B.10) é:
τey (x3) =
−k(x3 − θi) : x3 < θi
0 : θi ≤ x3 ≤ θf
−k(x3 − θf ) : x3 > θf
(B.11)
e o valor dos ângulos que limitam o movimento de nutação são:
θi = −50π/180 rad; θf = 35π/180 rad (B.12)
Para encontrar as componentes da matriz jacobiana aij são necessárias as derivadas par-ciais das funções fi da equação diferencial, ou seja, ∂fi
∂xie as derivadas parciais dos torques.
Começando então pelas últimas:
Derivada do torque τx:
τx = τx(x2, x3, x6) = −Tx sgn (x6) + [Tz sgn (x2) + bxx2] sin(x3)− bxx6 (B.13)
como τx não depende de x1, x4 e x5, como observado da eq. (B.13), obtém-se de forma diretaos resultados:
∂τx∂x1
=∂τx∂x4
=∂τx∂x5
= 0 (B.14)
99
para os outros casos tem-se:∂τx∂x2
= bx sin(x3) (B.15)
∂τx∂x3
= [Tz sgn (x2) + bxx2] cos(x3) (B.16)
∂τx∂x6
= −bx (B.17)
A derivada parcial do torque elástico é:
de =∂τe∂x3
=
−k : x3 < θi
0 : θi ≤ x3 ≤ θf
−k : x3 > θf
(B.18)
como τe só depende de x3:
∂τe∂x1
=∂τe∂x2
=∂τe∂x4
=∂τe∂x5
=∂τe∂x6
= 0 (B.19)
e a derivada parcial do torque τy é:
τy = τy(x3, x4) = MgH cos(x3)− Ty sgn (x4)− byx4 + τey (x3) (B.20)
∂τy∂x1
=∂τy∂x2
=∂τy∂x5
=∂τy∂x6
= 0 (B.21)
usando a eq. (B.18):∂τy∂x3
= −MgH sin(x3) + de (B.22)
∂τy∂x4
= −by (B.23)
Já derivando em relação a τz tem-se:
τz = τz(x2, x3) = −[Tz sgn (x2) + bzx2] cos(x3) (B.24)
∂τz∂x1
=∂τz∂x4
=∂τz∂x5
=∂τz∂x6
= 0 (B.25)
∂τz∂x2
= −bz cos(x3) (B.26)
∂τz∂x3
= [Tz sgn (x2) + bzx2] sin(x3) (B.27)
100
Com os resultados das derivadas acima, pode-se calcular as derivadas parciais ∂fi∂xj
, come-çando pelas funções fi mais simples tem-se:
f1 = x2; f3 = x4; f5 = x6 (B.28)
∂f1∂x1
=∂f1∂x3
=∂f1∂x4
=∂f1∂x5
=∂f1∂x6
= 0 (B.29)
∂f3∂x1
=∂f3∂x2
=∂f3∂x3
=∂f3∂x5
=∂f3∂x6
= 0 (B.30)
∂f5∂x1
=∂f5∂x2
=∂f5∂x3
=∂f5∂x4
=∂f5∂x5
= 0 (B.31)
∂f1∂x2
=∂f3∂x4
=∂f5∂x6
= 1 (B.32)
Calculando ∂f2∂xj
utilizando as eqs. (B.25), (B.26) e (B.27), partindo da função f2 = τz +
x4[Ixxx6 + (−Ix + Iy + Iz)x2 sin(x3)]/[Iz cos(x3)], ou seja, f2 = f2(x2, x3, x4, x6) tem-se:
a21 =∂f2∂x1
= 0 (B.33)
a22 =∂f2∂x2
=−bz + (−Ix + Iy + Iz)x4 tan(x3)
Iz(B.34)
a23 =∂f2∂x3
=[Tz sgn (x2) + bzx2] tan(x3) + (−Ix + Iy + Iz)x2x4
Iz+ f2 tan(x3) (B.35)
a24 =∂f2∂x4
=Ixxx6 + (−Ix + Iy + Iz)x2 sin(x3)
Iz cos(x3)(B.36)
a25 =∂f2∂x5
= 0 (B.37)
a26 =∂f2∂x6
=Ixxx4
Iz cos(x3)(B.38)
Para encontrar a derivada parcial de f4 usa-se as eqs. (B.21), (B.22) e (B.23), com a funçãof4 = f4(x2, x3, x6) = [τy − Ixxx2x6 cos(x3) + (Ix − Iz)x22 cos(x3) sin(x3)]/Iy obtém-se:
a41 =∂f4∂x1
= 0 (B.39)
101
a42 =∂f4∂x2
=−Ixxx6 cos(x3) + 2x2(Ix − Iz) sin(x3) cos(x3)
Iy(B.40)
a43 =∂f4∂x3
=(Ixxx2x6 −MgH) sin(x3) + de
Iy+
(Ix− Iz)(2 cos2 x3 − 1)x22Iy
(B.41)
a44 =∂f4∂x4
=−byIy
(B.42)
a45 =∂f4∂x5
= 0 (B.43)
a46 =∂f4∂x6
=−Ixxx2 cos(x3)
Iy(B.44)
Por último encontra-se ∂f6∂xj
utilizando as eqs. (B.14), (B.15), (B.16) e (B.17), com f6 =
f6(x2, x3, x4) = 1Ixx
[τx + (Ix + Iy − Iz)x2x4 cos(x3) + Ixf2 sin(x3)] resultado:
a61 =∂f6∂x1
= 0 (B.45)
a62 =∂f6∂x2
=(bx + Ixa22) sin(x3) + (Ix + Iy − Iz)x4 cos(x3)
Ixx(B.46)
a63 =∂f6∂x3
=[Tz sgn (x2) + bxx2 + Ixf2] cos(x3) + [Ixa23 − (Ix + Iy − Iz)x2x4] sin(x3)
Ixx(B.47)
a64 =∂f6∂x4
=(Ix + Iy − Iz)x2 cos(x3) + Ixa24 sin(x3)
Ixx(B.48)
a65 =∂f6∂x5
= 0 (B.49)
a66 =∂f6∂x6
=−bx + Ixa26 sin(x3)
Ixx(B.50)
Finalmente a matriz jacobiana é encontrada:
J =
0 1 0 0 0 0
0 a22 a23 a24 0 a26
0 0 0 1 0 0
0 a42 a43 a44 0 a46
0 0 0 0 0 1
0 a62 a63 a64 0 a66
, (B.51)
102
os valores de a22, a23, etc, não foram substituídos explicitamente na equação (B.51), pois amesma não caberia na página.
Apêndice C
Teste dos Programas de AnáliseNão-Linear
Para verificar o bom funcionamento dos métodos empregados no capítulo 5, é usado parateste um sistema dinâmico não-linear bem conhecido, no caso o sistema adotado é o de Lo-rentz (Lorenz, 1963).
As equações diferenciais desse sistema dinâmico não-linear são dadas a seguir:
X =
x
y
z
=
σ(y − x)
x(ρ− z)− yxy − βz
(C.1)
e sua jacobiana
J =
−σ σ 0
(ρ− z) −1 −xy x −β
. (C.2)
Com base nas equações diferenciais, constata-se que o sistema de Lorentz possui 3 pontos
de equilíbrioXe =xe ye ze
T:
Xe1 =
0
0
0
, Xe2 =
√β(ρ− 1)√β(ρ− 1)
ρ− 1
, Xe3 =
−√β(ρ− 1)
−√β(ρ− 1)
ρ− 1
, (C.3)
com ρ > 1.
Para verificar a estabilidade desses pontos de equilíbrio, basta fazer o estudo do sistemalinearizado η = Jη em torno dos pontos de equilíbrio, isto é feito através da verificação dosautovalores do sistema linearizado det (J(Xe)− γI) = 0, mais abaixo a estabilidade dospontos de equilíbrio será verificada para dois casos particulares.
Para o sistema de Lorentz será adotado a condição inicial X0 =
0,0 1,0 0,0T
com
103
104
tempo total de evolução (integração da trajetória fiducial) de 200,0 s e os parâmetros σ = 10,0
e β = 8/3, os métodos que serão empregados a seguir levam em conta que ρ pode variar de 1,0
até 100,0, os valores de ρ são variados com passo de uma unidade (1,0).
Para os expoentes de Lyapunov utilizando o método padrão e o método de Eckmann-Ruelle, o maior expoente pode ser visto na figura C.1, da figura pode-se ver que as diferençasentre os métodos é pequena, tal diferença deve ficar menor utilizando-se um tempo de evoluçãomaior. Com base no gráfico da figura C.1, constrói-se os atratores e verifica-se a estabilidadedos pontos de equilíbrio para dois valores particulares de ρ.
0 20 40 60 80 100−2
−1.5
−1
−0.5
0
0.5
1
1.5
2
2.5
3
Parâmetro ρ
Exp
oent
es d
e Ly
apun
ov [1
/ s]
Método padrãoMétodo de Eckmann−Ruelle
Figura C.1: Maior expoente de Lyapunov para o sistema de Lorenz, usando os métodos padrãoe de Eckmann-Ruelle.
Para ρ = 13,0 tem-se um comportamento regular (λ < 0), o atrator pode ser visto nafigura C.2. Para a estabilidade dos pontos de equilíbrio tem-se:
Xe1 =
0 0 0T
: γ1 = −17,7577, γ2 = 6,7577, γ3 = −2,6667, (C.4)
este ponto de equilíbrio é instável pois a parte real de um dos autovalores é positiva;
105
Xe2 =
5,6569 5,6569 12,0000T
:
γ1 = −12,7848, γ2 = (−0,4409 + 7,0615i), γ3 = (−0,4409− 7,0615i), (C.5)
Xe3 =−5,6569 −5,6569 12,0000
T:
γ1 = −12,7848, γ2 = (−0,4409 + 7,0615i), γ3 = (−0,4409− 7,0615i), (C.6)
para os pontos Xe2 e Xe3 , que possuem os mesmos autovalores, são estáveis pois a parte realde todos os autovalores é negativa e há dois autovalores complexos conjugados.
05
1015 0
510
150
5
10
15
20
YX
Z
Figura C.2: Atrator do sistema de Lorenz com comportamento regular e com condição inicialX0 =
0,0 1,0 0,0
T e parâmetros σ = 10,0, β = 8/3 e ρ = 13,0.
Já para ρ = 28,0 tem-se um comportamento caótico (λ > 0), o atrator correspondenteestá na figura C.3. A estabilidade dos pontos de equilíbrio é visto abaixo:
Xe1 =
0 0 0T
:
γ1 = −22,8277, γ2 = 11,8277, γ3 = −2,6667, (C.7)
para este ponto o equilíbrio é instável pois a parte real de um dos autovalores é positiva;
Xe2 =
8,4853 8,4853 27,0000T
:
γ1 = −13,8546, γ2 = (0,0940 + 10,1945i), γ3 = (0,0940− 10,1945i); (C.8)
106
Xe3 =−8,4853 −8,4853 27,0000
T:
γ1 = −13,8546, γ2 = (0,0940 + 10,1945i), γ3 = (0,0940− 10,1945i), (C.9)
para os pontos Xe2 e Xe3 , que possuem os mesmos autovalores, são instáveis pois a parte realdos autovalores complexos conjugados é positiva.
−20−10
010
20 −40−20
020
400
10
20
30
40
50
YX
Z
Figura C.3: Atrator do sistema de Lorenz com comportamento caótico e com condição inicialX0 =
0,0 1,0 0,0
T e parâmetros σ = 10,0, β = 8/3 e ρ = 28,0.
Para os parâmetros ρ = 13,0 e ρ = 28,0 serão criadas duas séries temporais da variável x,através da solução numérica da equação diferencial com tempo de amostragem de δt = 0,001 s.As séries temporais são utilizadas na obtenção do maior expoente de Lyapunov usando o métodode Wolf com os programas lyapmax e wolf (veja seção 5.2.1 na página 70) e também será usadano teste 0 − 1. Os resultados desses métodos juntamente com o resultado do método padrãocomo referência, está resumido na tabela C.1.
Tabela C.1: Resumo da série temporal simuladas da variável x do sistema de Lorentz. Na tabelaτ representa tempo de atraso, m representa dimensão de imersão, λl expoente de Lyapunovobtida do programa lyapmax, λw expoente de Lyapunov obtida do programa wolf, λp expoentede Lyapunov pelo método padrão e K é o resultado do teste 0− 1.
ρ τ m λl λw λp K Comportamento13,0 156 2 −0,118 −0,622 −0,617 0,274 regular28,0 159 3 1,15 0,0497 1,19 0,999 caótico
107
De acordo com a tabela C.1, pode se verificar que o programa wolf apresenta resultadomais próximo do valor de referência do método padrão com ρ = 13,0, já para ρ = 28,0 oprograma lyapmax apresenta um resultado mais próximo do respectivo valor de referência, nogeral os resultados do programa lyapmax estão melhores. Para o teste 0 − 1 o resultado para ocaso caótico foi um valor próximo de 1,0, já para o caso regular o valor não é tão próximo dezero assim, com uma série temporal com evolução temporal maior o valor de K deve tender azero mais rapidamente.
