Analise Dinamica Pelo Metodo Dos Elementos Finitos

Análise Dinâmica pelo Método de ElementosFinitos

Prof. Paulo de Tarso R. Mendonça, Ph.D.Universidade Federal de Santa Catarina - UFSCGrante - Depto. Engenharia Mecânica, CP 476,

Florianópolis, SC. CEP 88040-900,[email protected]

23 de Agosto de 2006

Resumo

O estudo que relaciona as forças que atuam sobre um corpo com o movimento, tanto do corpocomo um todo quanto de suas partes relativamente umas às outras, é denominado dinâmica.As equações que representam este movimento em velocidades não relativisticas são as leis domovimento de Newton. Um tipo particular de comportamento dinâmico é o movimento vi-bratório ou simplesmente a vibração, onde o sistema oscila em torno de uma certa posiçãode equilíbrio. Este texto lida com a simulação numérica por elementos Þnitos de vibrações emcorpos sólidos. Nas primeiras seções desse capítulo e no próximo, apresentamos uma revisãodos fundamentos de vibrações mas deve-se ter claro que este não é um substitutivo a um cursoformal no assunto. O objetivo destes tópicos aqui consiste apenas em homogeneizar o materiale a notação de forma a permitir o tratamento do problema pelo método dos elementos Þnitos.

Conteúdo1 Equação do Movimento em um Grau de Liberdade 3

2 Vibrações Livres de Sistemas não-Amortecidos 4

3 Vibrações Livres de Sistemas Amortecidos 5

4 Carregamento Harmônico 9

5 Resposta a Carregamentos não-Periódicos 125.0.1 Exemplo 1 Sistema Amortecido sob Carregamento Exponencial . . . . . . . 16

6 Sistemas com mais de um Grau de Liberdade 19

7 Elementos Finitos em Dinâmica 217.1 Princípio de DAlembert . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 217.2 Princípio dos Trabalhos Virtuais em Barras . . . . . . . . . . . . . . . . . . . . . . . 22

1

CONTEÚDO 2

7.3 Matrizes de Elementos Finitos de Barras . . . . . . . . . . . . . . . . . . . . . . . . . 237.4 Equações do Movimento de Lagrange . . . . . . . . . . . . . . . . . . . . . . . . . . . 26

7.4.1 Exemplo 2 Equações de Movimento com E.F. de Barra . . . . . . . . . . . . 27

8 Análise Modal 288.1 Vibrações Livres Não-Amortecidas . . . . . . . . . . . . . . . . . . . . . . . . . . . . 308.2 Propriedades dos Autovetores e Autovalores . . . . . . . . . . . . . . . . . . . . . . . 34

8.2.1 Ortogonalidade . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 348.2.2 Normalização e Ortonormalidade . . . . . . . . . . . . . . . . . . . . . . . . . 358.2.3 Exemplo 3 Freqüências Naturais . . . . . . . . . . . . . . . . . . . . . . . . . 378.2.4 Exemplo 4 Modos Naturais . . . . . . . . . . . . . . . . . . . . . . . . . . . 388.2.5 Exemplo 5 Normalização de Autovetores . . . . . . . . . . . . . . . . . . . . 398.2.6 Exemplo 6 Solução Analítica de Vibrações Livres em Barra . . . . . . . . . . 408.2.7 Autovetores Linearmente Independentes . . . . . . . . . . . . . . . . . . . . . 42

8.3 Análise Modal para Excitação Inicial - Sistema não-Amortecido . . . . . . . . . . . . 438.3.1 Exemplo 7 Resposta para Deslocamento Inicial pelo MEF . . . . . . . . . . 478.3.2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49

8.4 Análise Modal Geral . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 518.4.1 Exemplo 9 Solução pelo MEF de barra sob Carga Variável no Tempo . . . . 52

8.5 Resumo do Método de Análise Modal . . . . . . . . . . . . . . . . . . . . . . . . . . . 53

9 Determinação do Amortecimento 569.1 Um Grau de Liberdade . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 569.2 n-Graus de Liberdade Elementos Finitos . . . . . . . . . . . . . . . . . . . . . . . . 579.3 Métodos Experimentais . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58

9.3.1 Exemplo 10 Determinação Experimental de Amortecimento . . . . . . . . . 599.4 Método Analítico 1 Rayleigh . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 609.5 Método Analítico 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63

9.5.1 Exemplo 11 Determinação Experimental da Matriz de Amortecimento . . . . 639.5.2 Exemplo 12 Vibração Amortecida de Barra sob Deslocamento Inicial . . . . 649.5.3 Exemplo 13 Vibração Forçada Amortecida pelo MEF . . . . . . . . . . . . . 65

10 Lista de Exercícios 6610.0.4 Exemplo 14 Problema não-linear - Pêndulo duplo . . . . . . . . . . . . . . . 74

11 Princípio de Hamilton 7511.1 Cálculo variacional . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75

11.1.1 Operador variação . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7911.2 Energia cinética de um sistema e coordenadas generalizadas . . . . . . . . . . . . . . 8011.3 Princípio de Hamilton . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80

11.3.1 Equação de Lagrange para sistemas concervativos . . . . . . . . . . . . . . . . 8111.4 Sistemas contínuos - vibrações livres em vigas . . . . . . . . . . . . . . . . . . . . . . 8211.5 Vibrações livres de vigas pelo MEF e princípio de Hamilton . . . . . . . . . . . . . . 84

12 Lista de Exercícios 85

Paulo de Tarso R Mendonça - Grante - EMC - UFSC, 2004 3

1 Equação do Movimento em um Grau de LiberdadeConsideremos o sistema mecânico ilustrado na Figura 1, denominado sistema massa-amortecedor-mola, ou K-C-M . Observemos que nosso interesse Þnal consiste em analisar o comportamento decorpos sólidos, estruturas. O sistema mostrado na Figura em geral é tomado apenas como uma repre-sentação, uma idealização em vez de um sistema real. Como será visto, a formulação e a compreensãodo comportamento deste sistema é a peça básica sobre a qual são construídas as formulações dossistemas contínuos de corpos tridimensionais.

k

c

F(t)m

x(t)

m

(a) (b)

Fk(t)

dF (t)

F(t)

Figura 1: (a) Sistema idealizado k-c-m exitado, (b) diagrama de corpo livre.

No sistema da Figura 1, a massa m é considerada rígida, a mola de rigidez linear k é consideradasem massa e o amortecedor linear da constante c é considerado sem massa ou rigidez. Estes sãoevidentemente idealizações. A rigidez da mola signiÞca que Fk = k δk onde δk é o deslocamento entreas extremidades da mola provocado pela força Fk. O amortecedor é tal que, Fd = c úδd, onde úδd é avelocidade de afastamento entre as extremidades do amortecedor, como ilustrado na Figura 2.

c

c

dk δ

(a) (b)

δdcdF =

k

kδk=kF

δ k

Figura 2: (a) Força numa mola proporcional ao deslocamento; (b) força num amortecedor propor-cional à velocidade.

Observando o diagrama de corpo livre na Þgura e usando a segunda lei de Newton, obtemosa equação de movimento do sistema como:

F (t)− Fk(t)− Fd(t) = mx(t), (1)

onde x(t) é o deslocamento da massa, medido a partir da posição de equilíbrio. Essa posiçãocorresponde à situação onde a mola está descarregada. Substituindo as expressões para as forçasobtemos a equação do movimento na forma:

mx(t) + c úx(t) + k x(t) = F (t). (2)


Essa é uma equação diferencial linear, ordinária de coeÞcientes constantes m, c e k, que deÞnemas características do sistema físico sendo simulado. O carregamento aplicado sobre o sistema érepresentado pela força F (t), função do tempo t.

2 Vibrações Livres de Sistemas não-Amortecidos

O chamado problema de vibrações livres é aquele em que o sistema se move em ausência deforças de excitação, isto é, quando na eq. (2) se tem F (t) = 0 para todo t > 0. Nesse caso, a eq. (2)é dita estar em sua forma homogênea. Fisicamente, um sistema pode permanecer em movimentodurante algum tempo após a aplicação e subsequente remoção de força. Também é possível colocá-lo em movimento aplicando um deslocamento ou velocidade de curta duração. Por outro lado, asolução deste problema fornece subsídeos para a solução de problemas excitados, o que constitui aoutra razão pela qual ele é sempre estudado.É costumeiro reescrever a equação de movimento (2) em sua forma homogênea não-amortecida

como

..x(t) + ω2n x(t) = 0, ω2n =

k

m. (3)

A solução deste problema é conhecida e tem a forma geral

x(t) = A1 cosωnt+A2 senωnt (4)

onde A1 e A2 são constantes de integração a serem determinados a partir dos valores dados dodeslocamento e velocidade iniciais x(0) e

.x(0). Usando relações trigonométricas, (cos(a − b) =

cos a cos b+sena sen b), a solução também pode ser posta na forma equivalente

x(t) = A cos (ωnt− φ) (5)

com

A =A1cosφ

=A2senφ

, tanφ =A2A1

(6)

As novas constantes A e φ tem signiÞcado físico mais evidentes que A1 e A2, são a amplitudee ângulo de fase do movimento. O sistema realiza uma oscilação harmônica simples com fre-qüência natural ωn (em rad/segundo), isto é, a massa move-se para frente e para trás sempre coma mesma amplitude A e com freqüência de ωn/2π ciclos por segundo. O tempo gasto em cada ciclo,o período, é T = 2π/ωn segundos.DeÞnindo as condições iniciais °°°° x(0) = xo,.

x(0) = υo,(7)

aplicando-as em (4), calculam-se as constantes A1 e A2. A solução (4) então torna-se

x(t) = xo cosωnt+υoωnsenωnt. (8)

A amplitude e o ângulo de fase são


A =

sx2o +

µυoωn

¶2, tanφ =

υoxoωn

. (9)

Além das formas (4) e (5) existe ainda uma terceira forma para a solução do problema (3), dadapor

x(t) = A est, (10)

onde A e s são constantes a serem determinadas. Derivando (10) é possível ver que ela realmentesatisfaz a equação diferencial (3). Fazendo a substituição obtemos

s2A est + ω2nA est = 0. (11)

Como A e est são não nulos para uma solução não trivial, podemos dividir toda a equação por A est,obtendo a chamada equação característica do problema:

s2 + ω2n = 0. (12)

Essa equação tem duas soluções, dadas por

s = ± iωn, (13)

onde i =√−1 é a unidade complexa A solução da equação do movimento é então uma combinação

linear das duas formas resultantes da substituição de (13) em (10):

x(t) = A1eiωn +A2e

−iωn. (14)

3 Vibrações Livres de Sistemas AmortecidosQuando o amortecimento do sistema não é nulo, a equação de movimento (2) é reescrita para umaforma análoga a (3):

x(t) + 2ζωn úx(t) + ω2n x(t) = 0 (15)

onde

ζ =c

2mωn. (16)

é o chamado quociente de amortecimento viscoso. ζ tem signiÞcado físico deÞnido e será vistona seção 9. A solução do problema é aquela mostrada na eq. (10). Substituindo a solução (10) em(15) e simpliÞcando obtemos a equação característica

s2 + 2ζωn s+ ω2n = 0. (17)

Dois valores de s satisfazem esta equação:

s1s2

¾=

µ−ζ ±

qζ2 − 1

¶ωn. (18)


Cada raiz produz uma solução em (10). Da teoria de equações diferenciais lineares, temos que asolução do problema é uma combinação linear das soluções:

x(t) =hA1e

√ζ2−1ωnt +A2e−

√ζ2−1ωnt

ie−ζωn t. (19)

Os coeÞcientes A1 e A2 podem ser complexos. Claramente as raízes s1 e s2 podem ser reais oucomplexas dependendo do valor do amortecimento ζ. Por exemplo, para 0 < ζ < 1 a solução (19) de(15) pode ser posta na forma

x(t) =£A1e

iωdt +A2e−iωdt¤ e−ζωnt, (1a. forma), (20)

onde ωd = ωnp1− ζ2 é a chamado freqüência da vibração livre amortecida.

Como A1 e A2 podem ser complexos, pode-se coloca-los naforma polar A1 = |A1| eiα e A2 =|A2| eiβ, onde |A1| e |A2| são os módulos e α e β os argumentos. Então (20) Þca

x(t) =£|A1| ei(ωdt+α) + |A2| ei(−ωdt+β)¤ e−ζωnt. (21)

Usando as formas e±ic = cos c± i sen c, temos

x(t)eζωnt = |A1| [cos(ωdt+ α) + i sen (ωdt+ α)] + |A2| [cos(−ωdt+ β) + i sen (−ωdt+ β)] . (22)

Expandindo asfunçõestrigonométricas e separando as partes real e imaginária tem-se:

x(t)eζωnt = [(|A1| cosα+ |A2| cosβ) cosωdt+ (− |A1| sen α+ |A2| sen β) sen ωdt]+i [(|A1| cosα− |A2| cosβ) sen ωdt+ (|A1| sen α+ |A2| sen β) cosωdt] . (23)

Toma-se a solução do problema físico de vibrações livres como aparte real ou a parteimaginária de x(t). De (23) nota-se que amas são deÞnidas pela mesma base de funções, sen ωdt ecosωdt. (Que essas funções são solução pode ser provado simplesmente substituindo-as na equaçãodiferencial.) Então,

x(t) = [C1 cosωdt+ C2 sen ωdt] e−ζωnt, (2a. forma)úx(t) = −ζωn [C1 sen ωdt+ C2 cosωdt] e−ζωnt

+ωd [C1 cosωdt− C2 sen ωdt] e−ζωnt(24)

onde C1 e C2 são constantes reais, a serem determinadas pelas condições iniciais.

Outra forma útil para a solução pode ser obtida observando que o termo entre chaves na eq.(20)pode ser visto como uma soma vetorial no plano complexo, cujo resultado é também um vetor,comcomprimento A e ângulo φ em relação ao vetor eiωdt. Assim, em vez da forma (20), a solução podeser posta como

x(t)eζωnt = A ei(ωdt−φ) (25)

onde A é um coeÞciente real. Colocando em forma polar e separando as partes real e imagináriatem-se:


x(t)eζωnt = A cos(ωdt− φ) + iA sen (ωdt− φ),

=

⎡⎣A cosφ| z c1

cosωdt+A sen φ| z c2

sen ωdt

⎤⎦+i

⎡⎣A cosφ| z d1

sen ωdt−A sen φ| z d2

cosωdt

⎤⎦ . (26)

Se a solução do problema físico for tomada como a parte real de x(t), tem-se

x(t) = e−ζωnt [c1 cosωdt+ c2 sen ωdt] . (27)

Mas essa é exatamente a forma (24). A relação entre as constantes c1 e c2 com A e φ é dada em (26)como ½

c1 = A cosφc2 = A sen φ

=⇒(

tanφ =c2c1

A2 = c21 + c22.

(28)

Com isso, a parte real de (25) resulta em

x(t) = A e−ζωnt cos(ωdt− φ), (3a. forma)úx(t) = −A e−ζωnt [ζωn cos(ωdt− φ) + ωd sen (ωdt− φ)] (29)

Observe que em todas estas manipulações estamos apenas recombinando e mudando o signiÞcadodas constantes, mas elas permanecem sempre duas, ou A1 e A2, ou A e φ, ou c1 e c2. Cada uma destasdiferentes formas para a solução é mais adequada a diferentes tipos de interpretações e manipulaçõesalgébricas, algumas das quais faremos uso ao longo do texto. Esta última forma é obtida usando asdeÞnições de c1 e c2 em (28), que se torna:Caso 1 Para 0 < ζ < 1, velocidade inicial prescrita ( úx(0) = vo) e deslocamento inicial nulo

(x(0) = xo = 0). De (29) temos que cosφ = 0, o que dá φ = π/2. Também, a amplitude é A = vo/ωde a solução Þca

x(t) =voωde−ζωnt senωdt, ωd = ωn

p1− ζ2. (30)

Essa é a resposta do sistema à velocidade inicial vo, denominada solução transiente doproblema. O termo transiente se refere ao fato de que ela consiste em uma função periódica, oseno, que por si tem amplitude constante igual a 1 para todo t > 0. Porém o amortecimento no termoexponencial faz com que o fator multiplicando seno decresça ao longo do tempo. Assim, as oscilaçõesvão diminuindo de amplitude, como ilustrado na Figura 3. (Ali mostramos o caso de deslocamentoinicial, em vez de velocidade inicial, para melhor facilidade de visualização.)

Caso 2 Consideramos o problema (15) onde o sistema está submetido a um deslocamento inicialx(0) = xo, e velocidade inicial nula, úx(0) = 0. Da solução geral (29) temos:

x(0) = xo = A cosφ,úx(0) = vo = −Aζωn cosφ+Aωd senφ = 0. (31)


0.00 1.00 2.00 3.00 4.00 5.00

-4.00

-2.00

0.00

2.00

4.00

6.00

e−ξ ω nt

x(0)

t2πω

d

4πωd

x(t)

Figura 3: Esboço de solução de vibração transiente amortecida, eq.(35). (Usados os valores ωn =5 rad/s, ζ = 0, 2, xo = 5 mm, ωd = 4, 899, A = 5, 1, φ = 0, 201)

Este é um sistema de duas equações e duas incógnitas, A e φ. Elimina-se A cosφ na primeira equaçãocom uso da segunda, o que resulta em

tanφ =ζωnωd

e A =xocosφ

. (32)

Tomando ωd = ωnp1− ζ2 de (30), tem-se

tanφ =ζp1− ζ2

. (33)

De trigonometria temos que, da relação acima, cosφ =p1− ζ2. (considere um triângulo retângulo

com catetos ζ ep1− ζ2). Então, da segunda das eqs. (32),

A =xop1− ζ2

(34)

Então a solução do problema de vibração livre amortecida com x(0) = xo e úx(0) = 0 é obtidalevando (33) e (34) à solução (29) (ver Figura 3):

x(t) = xo√1−ζ2

e−ζωnt cos (ωdt− φ). (35)

Caso 3 Por Þm, consideramos o caso mais geral deÞnido pelo problema de valor inicial:°°°°°°x(t) + 2ζωn úx(t) + ω

2n x(t) = 0,

x(0) = uo,úx(0) = vo.

(36)


Uma vez que o problema é linear podemos simplesmente sobrepor a solução (30) obtida parax(0) = 0 e úx(0) = vo, com a solução (35) obtida para x(0) = uo e úx(0) = 0. Então a soluçãocompleta é

x(t) =

∙uo√1−ζ2

cos (ωdt− φ) + voωdsenωdt

¸e−ζωnt (37)

onde ωd e φ são deÞnidos em (30) e (33). Esta solução pode ainda ser posta nas formas alternativas:

x(t) =

∙vo + uoζωn

ωdsenωdt+ xo cos ωdt

¸e−ζωnt, ou ainda (38)

x(t) =

sµvo + uoζωn

ωd

¶2+ x2o cos (ωdt− φ) e−ζωnt, com tanφ =

vo + uoζωnxoωd

.

4 Carregamento HarmônicoConsideremos agora a solução particular do problema (2). O caso mais simples de carregamento é ochamado carregamento harmônico, isto é, um que varia segundo um seno ou cosseno ao longo dotempo:

F (t) = kf(t) = kA cosωt. (39)

A força aumenta e diminui ao longo do tempo com amplitude kA constante e freqüência constanteω. A eq. (2) pode ser dividida por k gerando uma equação similar a (15):

x(t) + 2ζωn úx(t) + ω2n x(t) = ω

2nA cosωt. (40)

A solução desse problema tem a seguinte forma

x(t) = X cos (ωt− φ) , (41)

onde X e φ são a amplitude e o ângulo de fase do movimento. Substituindo na eq. (40) obtemosa equação característica do problema

X£¡ω2n − ω2

¢cos (ωt− φ)− 2ζωnω sen (ωt− φ)

¤= ω2nA cosωt. (42)

Usam-se em seguida as seguintes relações trigonométricas½cos (ωt− φ) = cosωt cosφ+ senωt senφ,sen (ωt− φ) = senωt cosφ− cosωt senφ. (43)

Substiuindo em (42) pode-se igualar os coeÞcientes de cosωt de ambos os lados da igualdade e fazero mesmo com os coeÞcientes de senωt. Obtem-se então duas equações:½

X[(ω2n − ω2) cosφ+ 2ζωωn senφ] = ω2nA,X[(ω2n − ω2) senφ− 2ζωωn cosφ] = 0.


0.0 2.0

cos ψ

t

Figura 4: Visualização de um carregamento dado pela parte real de (45): f(t) = cos(ωt + ψ), comω = 5s−1 e ψ = 1, 2.

Essas equações podem ser resolvidas para as incógnitas do problema, X e φ:

X = A

⎡⎣Ã1−µ ωωn

¶2!2+

µ2ζω

ωn

¶2⎤⎦−1/2 e tanφ =2ζω/ωn

1−µω

ωn

¶2 . (44)

Levando (44) a (41) temos que o sistema responde com a mesma freqüência ω do carregamento,com uma amplitude X dependente da amplitude A do carregamento.

Uma outra forma de colocar a solução, consiste em usar uma forma complexa em lugar de(41). Representemos o carregamento por

f(t) = A ei(ωt+ψ). (45)

A é real e ψ o ângulo defase do carregamento. Usando a representação polar, o carregamentoconsiderado é da forma

f(t) = A cos(ω t+ ψ) + i A sen (ω t+ ψ). (46)

Obtemos a solução desse problema com o conhecimento prévio de que o carregamentoÞsicamente aplicado consiste ou na parte real de f(t), isto é, Re(f(t)) = A cos(ω t+ ψ), ouna parte imaginária, Im(f(t)) = sen (ω t+ ψ). Da solução complexa tomaremos também a partereal ou imaginária como representação física do movimento efetivamente realizado.Nesse aspecto, a dedução vista aqui é a mesma daquela para o carregamento cossenoidal da

eq.(39). Também, ali poderiamos ter feito a dedução para um carregamento em seno, em vez decosseno. Na presente dedução, pretende-se aumentar a generalidade, permitindo que o carregamentotenha uma fase, o que permite que o seno (ou cosseno) se inicie num valor diferente de sua posiçãonatural. Isso pode ser apreciado na Figura 4, que ilustra um carregamento dado pela parte real de(45). Note que os máximos da curva não ocorrem nas mesmas posições em que ocorreria no cosseno.A equação do movimento Þca

x(t) + 2ζωn úx(t) + ω2n x(t) = ω

2nA e

i(ωt+ψ). (47)


Toma-se a solução e suas derivadas na forma geral:

x(t) = X eiωt,úx(t) = iωX eiωt = iωx(t),x(t) = −ω2X eiωt = −ω2x(t).

(48)

Deve-se notar que, no caso geral, X é complexo. Isso signiÞca que, se colocamos X na formaX = X e−iφ, onde X é o módulo e φ o argumento, a solução se torna x(t) = X ei(ωt−φ). Então, φ éo ângulo de fase da resposta em relação à excitação.Substituindo (48) em (47), com X = X e−iφ, obtém-se

X£ω2n − ω2 + i(2ζωnω)

¤= ω2nAe

i(ψ+φ) (49)

Separando as partes real e imaginária temos duas equações reais:¯X (ω2n − ω2) = ω2nA cos(ψ + φ),X (ζωnω) = ω

2nA sen (ψ + φ),

que podem ser resolvidas para a amplitude X e para a fase φ:

tan(ψ + φ) =2ζω/ωn

1− (ω/ωn)2e X =

⎡⎣µ2ζωωn

¶2+

Ã1−

µω

ωn

¶2!2⎤⎦−1/2 . (50)

Note que essas são as mesmas expressões (44), porém agora válidas para um carregamento har-mônico com ângulo de fase ψ. Então, a parte real da solução correspondente a f(t) = A cos(ω t+ψ)é Re [x(t)], isto é:

x(t) = X Re£ei(ωt−φ)

¤= X cos(ω t− φ) (51)

Outra forma para a solução pode ser obtida substituindo (48) em (47), quando se determina asolução no tempo, correspondente a F (t) = kA cos(ωt+ ψ):

Re (x(t)) = Re

µω2nA e

i(ωt+ψ)

ω2n − ω2 + i2ζωωn

¶= Re

⎛⎜⎝ A ei(ωt+ψ)

1−³ωωn

´2+ i2ζ

³ωωn

´⎞⎟⎠

| z x(t)

.(52)

Observamos que a amplitude complexa X é uma função da freqüência do movimento, isto é, X =X(ω). Podemos multiplicar o numerador e o denominador do termo entre parênteses pelo conjugadocomplexo do denominador, o que resulta em:

x(t) =Aei(ωt+ψ)

(1− β2) + i(2ζβ) =A

(1− β2)2 + (2ζβ)2

⎧⎪⎨⎪⎩£(1− β2) cos(ωt+ ψ) + 2ζβ sen (ωt+ ψ)¤| z Re(x(t))

+ i£(1− β2) sen (ωt+ ψ)− 2ζβ cos(ωt+ ψ)¤ª (53)


onde β = ω/ωn. O termo real, o primeiro dentro do colchete, corresponde à solução transienteamortecida do carregamento F (t) = kA cos(ωt + ψ), e o termo imaginario corresponde ao carrega-mento F (t) = kA sen (ωt+ ψ).Então, a parte real da solução é:

x(t) =A£(1− β2) cos(ωt+ ψ) + 2ζβ sen (ωt+ ψ)¤

(1− β2)2 + (2ζβ)2 (Solução particular). (54)

A solução geral do sistema não-amortecido sob carregamento F (t) = kA cosωt (isto é, ψ = 0)é obtida sobrepondo a solução homogênea (24) e a solução particular (54):

x(t) = C1 cosωnt+ C2 senωnt| z Solução homogênea

+A cosωt

(1− β2)| z Solução particular

. (55)

Aplicando as condições iniciais x(0) = úx(0) = 0, obtém-se C2 = 0 e C1 = −A(1− β2), o que resultana solução:

x(t) =A

(1− β2) [cosωt− cosωnt]. (56)

5 Resposta a Carregamentos não-PeriódicosEstamos interessados em obter a resposta do sistema a um carregamento qualquer como o ilustradona Figura 5. Vários métodos existem para estimar a solução deste problema, e nos concentraremosaqui no método baseado na integral de convolução, também chamada integral de Duhamel.

F(τ)

0

F(t)

tt

Δττ

Figura 5: Carregamento temporal genérico.

Primeiramente introduzimos o conceito de função Delta de Dirac. Observe a função ilustradana Figura 6 é uma função ga(t) deÞnida em ∀t ∈ R, tal que¯

¯ ga(t) = 0, ∀ t < a, e ∀t > a+ E,ga(t) =

1

E, ∀ t ∈ [a; a+ E] . (57)


t

g (t)a

a ε

1ε_

Figura 6: Função com integral unitária.

É visível que

In =

Z ∞

−∞ga(y) dy = 1, ∀E ∈ R. (58)

Uma vez que a integral será sempre unitária para qualquer valor de E pode-se deÞnir uma pseudofunção (ou função generalizada), denominada função Delta de Dirac δ (t− a) , como

δ (t− a) = lim7→0

ga(t). (59)

Esta função é então

δ (t− a) = 0 para ∀t 6= a, (60)

e indeÞnida em t = a. Sua integral é tal que

Propriedade −→ R∞−∞ δ (t− a) dt = 1. (61)

Esta função é melhor compreendida como um operador com a seguinte propriedade (decorrentede (61)):

Propriedade −→ J =R∞−∞G(t)δ(t− a) dt = G(a), (62)

isto é δ(t− a) é tal que a quando multiplicado por qualquer função contínua e limitada tem integraligual ao valor desta função no ponto a. Esta característica pode ser entendida com a ajuda da Figura7. Note que o resultado do produto G(t) g(t− a) é não-nulo apenas no intervalo [a; a+ E]. Então

J =

Z a+7

a

G(t)lim7→0

ga(t) dt = lim1

E

Z a+7

a

G(t) dt = lim1

EG(a)E = G(a),

que resulta (62). Lembremos que as operações acima são apenas formais, e que resta provar algumasdelas, como a passagem do limite para fora da integral. Aquela relação constitui a principal utilidadeda função Delta de Dirac.

Consideremos agora um outro conceito, o impulso, deÞnido como

Impulso −→ I =R∞t=0F (t) dt. (63)


G(t)

1/ε

G(a)

t

a ε

G(t)

Figura 7: Funções G(t) e ga(t).

No momento nos interessa uma força de curta duração, uma força impulsiva como a função gailustrada na Figura 6 quando E → 0. Apesar da duração da força tender a zero, desejamos que seuimpulso, sua integral no tempo seja constante, igual a um valor dado F.Então a força impulsiva aplicada no instante a, de impulso F , é:

Força impulsiva −→ F (t) = Fδ (t− a). (64)

É visível que, substituindo esta equação em (63) temos I = F se a > 0. F tem unidade Newton·segundo.A resposta impulsiva h(t) é deÞnida como a resposta do sistema a uma força impulsiva unitária

aplicada no instante t = 0, isto é, h(t) = Fδ(t), com F = 1 Ns, sob condições iniciais nulas.Calculemos a resposta impulsiva à força F (t) = Fδ(t). Tomamos a eq. (2):°°°°°°

mx(t) + c úx(t) + kx(t) = Fδ(t), t > 0,x(0) = 0,úx(0) = 0.

(65)

Integra-se cada termo no intervalo t ∈ [0; E] e faz-se o limite E→ 0:

lim7→0

R 70Fδ(t) dt = F,

lim7→0+

R 70mx dt = lim

7→0+(m úx)|70 = lim

7→0+m ( úx (E)− úx (0)) = m úx(0+),

lim7→0+

R 70c úx dt = lim

7→0+(cx)|70 = lim

7→0+c (x (E)− x (0)) = 0,

lim7→0+

R 70kx dt = 0.

(66)

úx (0+) é a velocidade logo após o instante inicial. Note que em t = 0 a velocidade é nula, mas épossível aplicar uma variação de velocidade em um intervalo 4t bastante curto. Por outro lado, nãoé possível aplicar uma variação de deslocamento num intervalo de tempo curto. Assim, não apenasx (0) = 0, mas também x (0+) = 0.O resultado das integrações em (66) é

úx¡0+¢=F

m. (67)


Isto mostra que a aplicação de uma força impulsiva F produz uma mudança instantânea develocidade. Isto é interessante porque estamos interessados em obter a resposta do sistema a F . Masdevido a (67) podemos usar a solução que já dispomos para a resposta do sistema a uma velocidadeinicial. Usamos então a eq. (37) com vo = F/m e uo = 0:

x (t) =F

mωde−ζωnt sinωdt, ωd = ωn

¡1− ζ2¢1/2 , t > 0. (68)

Fazendo F = 1 temos a chamada resposta impulsiva

h(t) =1

mωde−ζωnt senωdt, t > 0,

h(t) = 0, t ≤ 0.(69)

Como o problema de vibrações sendo considerado é linear, a resposta do sistema a um impulsto Ié x(t) = I h(t). Então, a solução no instante t devido a um impulso I = F (0)4τ aplicado no instanteinicial τ = 0 é:

x(t) = F (0)4τ h(t). (70)

Essa é a resposta no instante t devido a um carregamento de curta duração, aplicado no instantet = 0. (Considere que a faixa hachurada na Figura 5 estivesse localizada na origem.)

Imagine-se agora o histórico de carga representado pela curva F (t) na Figura 5, substituido poruma seqüência de retângulos como aquele hachurado, cada um iniciando num instante τ , de duração4τ , e altura F (τ).Pergunta: qual a resposta no instante t devido a um impulso aplicado no instante τ? Observando

a eq.(70), notamos que a solução num dado instante t depende apenas do tempo decorrido entre oinstante do impulso e o instante t. Isto porque, obviamente, o corpo não sofre nenhum efeito doimpulso antes dele ter sido aplicado, isto é em t < τ . Mais ainda, a resposta do sistema no instantet depende apenas do lapso de tempo desde a aplicação do impulso, isto é, da extensão de tempodecorrido t − τ . Então, se o impulso foi aplicado no instante τ , a solução em t pode ser obtidasimplesmente usando a solução (70) substituindo t por (t− τ), isto é

4x(t) = F (τ) 4τ h (t− τ) = F (τ)

mωde−ζωn(t−τ) senωd (t− τ)4τ . (71)

Mas observe que no instante t o valor do deslocamento não é apenas devido a este intervalo decarregamento aplicado entre τ e τ + 4τ . Existem também parcelas provenientes dos impulsos deduração 4τ aplicados desde o instante 0 até t que compõem a curva F (t). Então a resposta total emt é

x(t) =X

F (τ) h (t− τ)4τ .. (72)

Fazendo 4τ → 0 o somatório tende à integral e temos

Integral de Convolução −→ x(t) =

Z t

τ=o

F (τ) h (t− τ) dτ . (73)


Esta integral aparece em diversas áreas das ciências físicas e é logicamente objeto de estudomatemático em busca de suas propriedades. É a chamada integral de convolução. Uma daspropriedades mais úteis desta integral, que apresentamos sem demonstração, é queZ t

τ=o

F (τ) h (t− τ) dτ =Z t

τ=o

F (t− τ) h(τ) dτ. (74)

Substituindo a deÞnição de h(t) em (73) temos a solução do sistema a um carregamento genérico:

x(t) =1

mωd

Z t

o

F (τ)e−ζωn(t−τ) senωd(t− τ) dτ. (75)

No estudo de vibrações esta é a chamada integral de Duhamel. Esta solução é chamadasolução particular ou solução de regime permanente. Lembremos que esta é apenas parte dasolução geral, válida no caso em que x(0) = úx(0) = 0.O problema geral °°°°°°

mx+ c úx+ κx = F (t), t > 0,x(0) = uo,úx(0) = vo,

(76)

tem solução obtida sobrepondo a solução de regime permanente (75) com a solução de transiente(37) obtendo a solução geral:

x(t) =uo e

−ζωntp1− ζ2

cos (ωdt− φ) + voωde−ζωnt senωdt+

1

mωd

Z t

0

F (τ)e−ζωn(t−τ) senωd(t− τ) dτ,

ωd = ωnp1− ζ2, ω2n =

k

m, tanφ = ζ√

1−ζ2.

(77)

5.0.1 Exemplo 1 Sistema Amortecido sob Carregamento Exponencial

Considere um sistema de um grau de liberdade, como na Figura 1a, amortecido, submetido a umcarregamento do tipo °°°° F (t) = mω2de−ζωnt para t ≥ 0,

F (t) = 0 para t < 0.(78)

Calcule a resposta do sistema para condições iniciais nulas, isto é, x(0) = úx(0) = 0.

Solução:Toma-se a solução geral do sistema, eq. (77). Para uo = vo = 0 Þca-se apenas com a integral de

convolução, que substituindo (78) Þca

x(t) =1

mωd

Z t

τ=0

F (τ) e−ζωn(t−τ) senωd (t− τ) dτ

= ωd e−ζωnt

Z t

o

senωd (t− τ) dτ.


0 2 4 6t

0

5

10

15

20

F(t

)

Figura 8: Carregamento exponencial do Exemplo 1 dado pela eq. (78), usando os valores:m =1 kg, ξ = 0, 2, ωn = 5 s−1, ωd = 4, 9 s−1.

Integrando e aplicando os limites temos a resposta.

x(t) = e−ζωnt [1− cosωdt] para t > 0. (79)

A Figura 9 ilustra a curva de resposta ao longo do tempo. É interessante notar que o movimentoda massa não é oscilante em torno do ponto de equilíbrio x = 0, mas sofre um movimento oscilanteonde a posição mínima é x = 0. A massa atinge essa posição periodicamente com período

tn =2nπ

ωd. (80)

0.0 1.0 2.0 3.0 4.0 5.0t

0.0

0.2

0.4

0.6

0.8

1.0

1.2

x(t)

Figura 9: Resposta do sistema do Exemplo 1 para o carregamento exponencial da eq.(78), usandoos valores: m =1 kg, ξ = 0, 2, ωn = 5 s−1, ωd = 4, 9 s−1.

Observemos algumas propriedades das soluções obtidas pela integral de Duhamel. Observea deÞnição em (77). Consideremos o caso em que o carregamento seja dado como uma combinaçãode dois outros carregamentos:


F (τ) = aF1(τ) + b F2(τ), (81)

como por exemplo as funções ilustradas nas Figuras 10a e b, onde F1(t) = senωt e F2(t) = 0 parat < t1 e F1(t) = senωt e F2(t) = senω(t − t1) para t ≥ t1. Podemos deÞnir uma função F (t) comoF (t) = F1(t)−F2(t), como ilustrado na Figura 10c. Da integral de convolução temos a resposta parao carregamento (81):

Z t

τ=0

F (τ) e−ζωn(t−τ) senωd (t− τ) dτ = a

Z t

τ=0

F1(τ) e−ζωn(t−τ) senωd (t− τ) dτ +

b

Z t

τ=0

F2(τ) e−ζωn(t−τ) senωd (t− τ) dτ. (82)

1

1

t t

t

t

t

t

F

(a)

(b)

(c)

1

F

1

x(t)

x (t)

F 2

t1

x 2(t)

1t

1

Figura 10: Solicitações e respostas com shift e sobreposição.

Uma vez que, freqüentemente, é impossível realizar analiticamente a integral de convolução deuma função arbitrária, em geral ela é realizada numericamente. Entretanto, algumas propriedadesda integral de convolução permitem alargar um pouco o leque de situações passíveis de ter soluçãoanalítica. Por exemplo, suponha que se tenha conseguido obter a solução para um carregamentoF1(t) como mostrado na Figura 10(a) e (b). Se transladarmos F1(t) em t1 e deÞnirmos assim afunção F2(t), a solução x2(t) é a solução x1(t) transladada, isto é, x2(t) = x1(t − t1) para t > t1 ex2(t) = 0 para t < t1. Assim, a solução associada a F = aF1 + b F2 (Figura 10c)é a combinação dasduas soluções:

x(t) = a x1(t) + b x2(t). (83)


Esta possibilidade de combinação é resultante da linearidade da equação do movimento usada. Estasoperações também necessitam que as condições iniciais uo e vo sejam combinadas da mesma formaatravés das mesmas constantes a e b.

6 Sistemas com mais de um Grau de Liberdade

Poucos são os sistemas físicos na engenharia que consistem realmente de um grau de liberdade comodescrito nas seções acima, composto um bloco rígido de massa m, ligado a uma base por uma molae um amortecedor. Estamos interessados principalmente em determinar o comportamento dinâmicodos sistemas contínuos, isto é, corpos e estruturas sólidas, tridimensionais, com sua forma própria,sua massa e sua capacidade de amortecimento interno. Entretanto, a teoria uniaxial vista nas seçõesacima é de fato usada como parte de vários métodos de análise de corpos tridimensionais, como serávisto nas seções que se seguem.Considere o corpo com forma genérica ilustrado na Figura 11a submetido a um conjunto de forças

variantes ao longo do tempo. Caso sua forma, apoios e carregamento sejam simples, regulares, é pos-sível uma modelagem analítica que resulta na solução exata de sua resposta. Alguns problemas ondeo corpo tenha forma de barra, vigas, placas circulares ou retangulares, dependendo do carregamento,podem ser tratadas desta forma. Uma série de livros clássicos tratam destas soluções, como por ex-emplo Langhaar [15], Meirovitch [18], Clough [4] e outros. Frequentemente porém, os componentese sistemas usados em engenharia são de forma e carregamento complexos e não podem ser trata-dos pelas fórmulas analíticas simples disponíveis. Da mesma forma que em problemas estáticos, amaneira hoje padrão de se tratar destes problemas consiste em abrir mão do desejo de obter umasolução exata e buscar uma solução aproximada do problema.Para tratar então do problema contínuo como o do corpo tridimensional da Figura 11a, cri-

aremos um modelos discretizado como o ilustrado na Figura 11b, onde o corpo é simulado por umacoleção de massas discretas unidas por molas e amortecidedores entre si. A forma de realizar esteprocesso de discretização não é óbvio, e existem diversos métodos, dentre os quais o próprio métodoque estamos tratando, o de elementos Þnitos. No momento consideramos que, de alguma forma,temos já realizado esta discretização e temos disponível um modelo como o da Figura 11b, com nmassas discretas. Cada uma dessas massas pode ser considerada um corpo rígido, de forma quepodemos aplicar a ela a segunda lei de Newton. A Figura 11c representa um diagrama de corpo livrede uma massa genéricami. Sobre ela atuam uma força externa Fi(t) e as forças internas provenientesdos deslocamentos relativos às outras massas. Estas forças internas são as forças elásticas fe, rela-cionadas à rigidez das molas Ki e Ki+1, e as forças de amortecimento fa relacionadas às constanteCi e Ci+1 dos amortecedores. Pela segunda lei de Newton, a resultante de todas estas forças deve serigual à força de inércia mixi. Então a equação do movimento para uma massa mi interna qualqueré a seguinte:

Fi + Ci+1 ( úui+1 − úui) +Ki+1 (ui+1 − ui)− Ci ( úui − úui−1)−Ki (ui − ui−1) = miui. (84)

Podemos rearranjar os termos colocando a parte conhecida, a força externa Fi(t), do lado direito:

miui − Ci+1 úui+1 + (Ci + Ci+1) úui − Ci úui−1 −Ki+1ui+1 + (Ki +Ki+1)ui −Kiui−1 = Fi. (85)

Podemos expandir estas equações e colocar o sistema em forma matricial:


(a)

(b

(c)

F (t) (t)FnF (t)1

F2(t)

1F (t)

mnmi

iF (t)

m1

i+1(t)

i-1F

u (t)i-1u (t)1u (t)i+1u (t)

i nu (t)

i+1mmi-1

k 1 ik

i+1c

nk

1c ic nc

k i+1

miik u (t)i i-1u (t)

i+1c u (t)i+1 u (t)i

(t)iF

u (t)i

m

-( )

( )-

u (t)( )u (t)-k i+1i+1 i

u (t)( )ic i u (t)- i-1

u i

Figura 11: a) Corpo sólido tri-dimensional qualquer; b) modelo discretizado; c) diagrama de corpolivre da massa mi indicando força externa aplicada, força de inércia, forças elásticas de mola e forçasde amortecimento.

⎡⎢⎢⎢⎢⎢⎣m1

m2

m3

. . .mn

⎤⎥⎥⎥⎥⎥⎦

⎡⎢⎢⎢⎢⎢⎢⎢⎣

u1u2u3u4...un

⎤⎥⎥⎥⎥⎥⎥⎥⎦+

⎡⎢⎢⎢⎢⎢⎣(C1 + C2) −C2−C2 (C2 + C3) −C3

−C3 (C3 + C4) −C4. . .

⎤⎥⎥⎥⎥⎥⎦

⎡⎢⎢⎢⎢⎢⎣úu1úu2úu3...úu1

⎤⎥⎥⎥⎥⎥⎦

+

⎡⎢⎢⎢⎢⎢⎢⎢⎣

(K1 +K2) −K2

−K2 (K2 +K3) −K3

−K3 (K3 +K4) −K4

−K4 (K4 +K5) −K5

. . .

⎤⎥⎥⎥⎥⎥⎥⎥⎦

⎡⎢⎢⎢⎢⎢⎢⎢⎣

u1u2u3u4...un

⎤⎥⎥⎥⎥⎥⎥⎥⎦=

⎡⎢⎢⎢⎢⎢⎢⎢⎣

F1(t)F2(t)F3(t)F4(t)...Fn(t)

⎤⎥⎥⎥⎥⎥⎥⎥⎦(86)


Este sistema pode então ser escrito de forma compacta como

M u(t) +C úu(t) +Ku(t) = F(t), (87)

que é a equação algébrica do movimento do sistema discreto da Figura 11b. As matrizes são aschamadas matriz massa ou de inérciaM, matriz de amortecimento C e a matriz de rigidez K, todassimétricas. Embora neste exemplo aM seja diagonal, de forma geral isto não é assim.

7 Elementos Finitos em DinâmicaNas próximas seções trataremos de métodos para determinar ou estimar a solução do problemaalgébrico (87). Antes disso daremos uma amostra do processo geral de como aquelas matrizes sãodeterminadas pelo método de elementos Þnitos na discretização de um corpo sólido, isto é, um sistemacontínuo.Basicamente, o processo de determinação da equação matricial de movimento num caso dinâmico

pelo método de elementos Þnitos é o mesmo procedimento usado nos capítulos anteriores na deter-minação da equação matricial de equilíbrio. Em ambos os casos usaremos o princípio do trabalhovirtual (PTV), onde no caso dinâmico fazemos uma alteração em seu enunciado pelo uso do chamadoprincípio de DAlembert, descrito a seguir. Um outro procedimento a ser apresentado, além doPTV, é a obtenção das equações de matriciais de movimento pelo uso das equações de movimentode Lagrange. Estas equações, porém, são apenas uma forma derivada do mesmo PTV aplicado àdinâmica.

7.1 Princípio de DAlembert

Julgando-se apenas pelo seu enunciado, este princípio é de uma simplicidade enorme. Sua utilidadeentretanto é também enorme na engenharia. Considere a equação do movimento de uma partículade massa m, dada pela segunda lei de Newton:

nXi=1

Fi +mb = ma, (88)

isto é, a resultante de todas as n forças externas aplicadas Fi, incluindo a força do corpomb, deve serigual a força da inércia, dada pela massa vezes a aceleração a desenvolvida pela massa m; b é umaforça de corpo por unidade de massa. Quando as forças são tais que a aceleração é nula, as forçasestão em equilíbrio e esta equação é chamada equação de equilíbrio. Obviamente, o tratamentode problemas de equilíbrio é mais simples que o tratamento de problemas dinâmicos. DAlembert,de certa forma, realizou uma operação bastante simples. Ele transferiu a força de inércia do ladodireito de (88) e passou-a para o lado esquerdo obtendo

nXi=1

Fi +m(b− a) = 0. (89)

Agora, a forma da equação é exatamente a mesma de uma equação de equilíbrio estático, e tudo oque existe desenvolvido para os problemas estáticos pode ser adaptado para os problemas dinâmicos.O princípio de DAlembert então aÞrma que as forças de inércia podem ser incorporadas às forças decorpo e o problema pode ser tratado como uma equação estática, embora, agora, as forças de corpoassociadas à inécia sejam desconhecidas, por incorporarem as acelerações.


7.2 Princípio dos Trabalhos Virtuais em Barras

Lembre que nos vários capítulos anteriores, o P.T.V. foi desenvolvido e aplicado aos diversos tiposde componentes para o comportamento estático. Com o uso do princípio de DAlembert as mesmasexpressões do P.T.V. podem ser expandidas ao problema dinâmico de forma bastante simples.

dxx

A

Figura 12: Elemento diferencial de volume de barra.

Tomemos por exemplo o P.T.V. para o problema estático de barras:

AE

Z L

o

du

dx

du

dxdx−A

Z L

o

bu dx−Af u(L) = 0, ∀u ∈ V ar, (90)

onde b é a força de corpo por unidade de volume, e f a força concentrada aplicada na extremindadex = L. Considere um elemento diferencial de volume numa barra, como ilustrado na Figura 12. Amassa deste elemento é dm = ρAdx onde ρ é a densidade do material, em kg/m3, por exemplo.

x u(x,t)

X(x,t)

F(0)

F(t)

Figura 13: Posição inicial P e posição Þnal p num dado instante t e deslocamento u(x, t) de umelemento diferencial numa barra sob solicitação dinâmica.

Observe na Figura 13 o comportamento dinâmico de uma barra sob carga axial. O elementodiferencial inicialmente encontra-se a uma distância x da origem. Num outro instante t a posiçãodaquela porção de material encontra-se a uma posição X da origem. A posição atual é função daposição inicial e varia instante a instante. Então podemos escrever que

X(x, t) = x+ u(x, t), (91)

isto é, a posição atual X do ponto é igual à posição inicial x mais o valor u(x, t) do deslocamentosofrido. Como a posição inicial não se altera, diferenciando (91) temos a velocidade e a aceleração:


⎧⎪⎪⎨⎪⎪⎩úX(x, t) =

∂u

∂t(x, t) = úu(x, t), (velocidade)

X(x, t) =∂2u

∂t2(x, t) = u(x, t) (aceleração).

(92)

Isto signiÞca que a taxa de variação da posição é a mesma do deslocamento. O elemento diferencialde massa sofre uma aceleração u(x, t) e sua força de inércia é

Aρu(x, t)| z binércia

dx, (força de inércia) (93)

que pode ser colocado na forma Abinérciadx, onde binércia = ρ u(x, t) é uma psudo-força de corpoassiciada à inércia. Assim, a força de inércia pode ser incluída no P.T.V. da equação (90) substituindoa força de corpo estática b(x) por (b(x, t)−binércia), isto é, por (b(x, t)− ρu(x, t)), levando à expressãodo PTV dinâmico

AE

Z L

o

∂u(x, t)

∂x

du(x)

dxdx−A

Z L

o

(b(x, t)− ρu(x, t)) u(x) dx−Af(t) u(L) = 0, ∀u ∈ V ar. (94)

7.3 Matrizes de Elementos Finitos de Barras

Como no caso estático, consideramos o problema de uma barra sujeita a uma força f na extremidadee forças de corpo b distribuídas ao longo de sua extensão, como ilustrado na Figura 12, mas agora asforças podem ser função do tempo. A solução do problema consiste na função u(x, t) que satisfaz àexpressão do PTV, eq.(94). A cada instante t a aceleração possui um valor, u(x, t) e os carregamentostem valores deÞnidos b(x, t) e f(t). Tem-se então o PTV estático neste instante, onde se deve buscara solução também estática, u(x, t). O tratamento por elementos Þnitos consiste ainda em discretizar ocorpo em elementos e aproximar os campos por funções de interpolação. Considere pois um elementoÞnito linear de dois nós, e suas funções de interpolação implicitas,

ϕe1(x) =Le − xLe

e ϕe2(x) =x

Le, (funcões de interpolação lineares) (95)

associados aos nós intrinsecos 1 e 2. (Le é o comprimento do elemento.) O campo de deslocamentoaxial no elemento e é então interpolado por

u(x, t) =2Xi=1

ui(t)ϕei (x). (96)

Observe que os deslocamentos nodais em (96), os u0is, agora variam com o tempo. De fato, essaexpanssão representa uma separação de variáveis, em tempo t e espaço x.A função peso é interpolada pela mesma base de funções de interpolação ϕei (x):

u(x) =2Xi=1

uiϕei (x). (97)

Entretanto, os valores nodais ui da função peso não variam no tempo. De (96), temos claramente aaproximação da aceleração no elemento:


u(x, t) =2Xi=1

ui(t)ϕei (x). (98)

Levamos essas expressões ao P.T.V. da eq. (94):

AE

Ze

Ã2Xi

ui(t)dϕeidx(x)

!| z

u(x,t)

Ã2Xi

uidϕeidx(x)

!| z

u(x)

dx−

A

Ze

Ãb(x, t)− ρ

2Xi

ui(t)ϕei (x)

!dx−Af(t) u2 = 0, ∀u ∈ V ar. (99)

As integrais são feitas ao longo do comprimento do elemento e. Faz-se o produto entre os doissomatórios na primeira integral e transfere-se o termo independente de x para fora da integral (oscoeÞcientes ui(t) e ui). Primeiramente obtém-se a forma expandida dos somatórios

AERe

µu1(t)

dϕe1dx

+ u2(x)dϕe2dx

¶| z

u(x,t)

µu1dϕe1dx

+ u2dϕe2dx

¶| z

u(x)

dx−A Reb(x, t) (u1ϕ

e1 + u2ϕ

e2) dx

−ρA Re(u1(t)ϕ

e1 + u2(t)ϕ

e2) (u1ϕ

e1 + u2ϕ

e2) dx−Af(t)u2 = 0 ∀u1, u2.

(100)

Como esta igualdade deve ser satisfeita para quaisquer valores de u1 e u2, podemos fazer u1 = 1e u2 = 0 e obter uma equação algébrica. Em seguida pode-se fazer u1 = 0 e u2 = 1, obtendo umasegunda equação:

=⇒∙AE

Ze

dϕe1dx

dϕe1dx

dx

¸| z

Ke11

u1(t) +

∙AE

Ze

dϕe1dx

dϕe2dx

dx

¸| z

Ke12

u2(t)+

∙ρA

Ze

ϕe1ϕe1 dx

¸| z

Me11

u1(t) +

∙ρA

Ze

ϕe1ϕe2 dx

¸| z

M212

u2(t) = AReb(x, t)ϕe1 dx,

=⇒∙AE

Ze

dϕe2dx

dϕe1dx

dx

¸| z

Ke21

u1(t) +

∙AE

Ze

dϕe2dx

dϕe2dx

dx

¸| z

Ke22

u2(t)+

∙ρA

Ze

ϕe2ϕe1 dx

¸| z

Me21

u1(t) +

∙ρA

Ze

ϕe2ϕe1 dx

¸| z

Me22

u2(t) = AReb(x, t)ϕe2 dx+Af(t).

(101)

Pode-se reconhecer as primeiras duas integrais em cada equação como os termos da matriz derigidez do elemento Þnito de análise estática de barra usado em análise estática, visto na eq. (??).Os termos no terceiro e quarto colchete em cada equação formam a chamada matriz massa oumatriz de inércia do elemento. Os termos à direita da igualdade formam o vetor de carregamento,o mesmo mostrado em (??). Então, as equações acima podem ser postas na forma


(Ke11u1(t) +K

e12u2(t) + Me

11u1(t) +Me12u2(t) = F

e1 (t),

Ke21u1(t) +K

e22u2(t) + Me

21u1(t) +Me22u2(t) = F

e2 (t),

(102)

ou, em notação matricial,

Meue(t) +Keue(t) = Fe(t). (103)

Comparando com o caso estático vemos que agora existe uma força de inércia, Meue(t), e queagora o carregamento pode variar no tempo. Mas mesmo que o carregamento seja estático, a presençado termo de inércia permitirá que a solução seja variante com o tempo. A deÞnição de cada termo éa seguinte ¯

¯¯¯

Keij = AE

Ze

dϕeidx

dϕejdx

dx,

Meij = ρA

Ze

ϕeiϕej dx,

F ei (t) = A

Ze

b(x, t)ϕi(x) dx+ f(t)ϕi(Le).

(104)

Para outros tipos mais complexos de elementos, como os de placa ou sólidos, as integrações acimapodem ser inviáveis de serem realizadas analiticamente e são realizadas de forma numérica. Aqui,entretanto, as funções de interpolação são simples, lineares como vistas na eq. (95). A integraçãoanalítica da matriz de inércia é a seguinte¯

¯¯¯

Me11 =

ρA

L2e

Z xe2

xe1

(Le − x)2 dx = ρALe3,

Me21 =M

e12 =

ρA

L2e

Z xe2

xe1

(Le − x)x dx = ρALe6,

Me22 =

ρA

L2e

Z xe2

xe1

x2 dx =ρALe3.

(105)

A matriz de rigidez e o vetor força do elemento são os mesmos já integrados nas equações (??) e(??), isto é,

Ke =AE

Le

∙1 −1

−1 1

¸; Me =

ρALe6

∙2 11 2

¸; Fe = Af

∙01

¸. (106)

Como no caso estático, essas são apenas as matrizes de um elemento, e devem ser rotacionadas esobrepostas nas matrizes globais para gerar as equações discretas de movimento que representamo sistema sendo modelado:

Mu(t) +Ku(t) = F(t). (107)

Observe que este é um sistema de n equações diferenciais, ordinárias, de coeÞcientes constantes,em termos do tempo, não homogêneo (F (t) 6= 0), com n funções incógnitas u1(t), . . . , un(t). Difer-entemente do caso estático, esta não é uma equação algébrica, portanto não pode ser simplesmenteresolvida por inversão de matriz. É um sistema de equações diferenciais e deve ser integrado paradar a resposta do sistema, após a aplicação das condições de contorno e condições iniciais do sistema.


7.4 Equações do Movimento de Lagrange

As equações do movimento (107), obtidas pelo PTV, podem ser também obtidas com o uso dasequações de Lagrange. Considere que podemos expressar a energia de deformação U e a energiacinética T de um corpo elástico ou sistema, em termos de valores nodais de deslocamento ui(t) eúui(t), isto é, se temos as funções

W = U (u1(t), u2(t), . . . , un(t)) ,

T = T (u1(t), u2(t), . . . , un(t), úu1(t), úu2(t), . . . , úun(t) .(108)

Partindo do princípio dos trabalhos virtuais, é possível deduzir as chamadas equações de Lagrange.Não apresentaremos aqui sua dedução, que pode ser achada em textos clássicos de dinâmica ??.Estas equações são as equações do movimento do sistema, em termos dos valores nodais ui(t) e úui(t).Frequentemente é mais fácil obter as equações do movimento através das equações de Lagrange quetentando aplicar diretamente a segunda lei de Newton. As equações de Lagrange são

d

dt

µ∂T

∂ úui

¶+∂T

∂ui+∂W

∂ui= Fi. (109)

Para uma barra, a energia de deformação é:

W =AE

2

Z L

o

µ∂u(x, t)

∂x

¶2dx. (110)

A energia cinética pode ser obtida da seguinte forma. Lembre que a energia cinética de umamassa pontual m é, por deÞnição, Ec = mv2/2. Agora considere o elemento diferencial de barra decomprimento dx das Figuras 12 e 13. Este elemento tem massa diferencial dm = ρAdx e velocidadeaxial úu(x, t). Então sua energia cinética é úu(x, t)2dm/2, isto é, ρA úu(x, t)2dx/2. A energia cinéticada barra completa é então

T =ρA

2

Z L

o

( úu(x, t))2 dx. (111)

Tendo a W e T , cabe agora fazer a discretização de elementos Þnitos. Dividimos a barra emelementos, o que signiÞca simplesmente particionar o intervalo de integração nas deÞnições de We T em uma soma de integrais realizadas em cada elemento. Em cada elemento interpolamos odeslocamento e velocidade usando (96). Então as energias em cada elemento se tornam:

W e =AE

2

Ze

µu1(t)

dϕe1dx

+ u2(t)dϕe2dx

¶2dx,

T e =ρA

2

Ze

( úu1(t)ϕe1 + úu2(t)ϕ

e2)2 dx.

(112)

Se usarmos as funções de interpolação lineares (95) no elemento essas energias Þcam1

W e =AE

2Le

∙u1(t)u2(t)

¸T ∙1 −1

−1 1

¸ ∙u1(t)u2(t)

¸,

T e =ρA

12Le

∙úu1(t)úu2(t)

¸T ∙2 11 2

¸ ∙úu1(t)úu2(t)

¸.

(113)

1O expoente T indica transposto de vetor ou matriz.


Se compararmos estas expressões a (106) vemos que as matrizes acima são proporcionais a rigideze massa do elemento respectivamente. DeÞnindo o vetor de deslocamentos nodais do elementocomo ue(t) = [u1(t);u2(t)]

T , a equação acima pode ser posta na forma½W e = 1

2ueT (t)Ke ue(t),

T e = 12úueT (t)Me úue(t).

(114)

Tendo então as expressões discretizadas para as energias interna e cinemática, podemos passarao uso das equações de Lagrange. Para um dado elemento a equação de Lagrange em (109) reduz-sea duas equações, para i = 1 e i = 2, correspondentes aos dois graus de liberdade do elemento. Asequações Þcam:

Equações de Lagrange

⎧⎪⎪⎨⎪⎪⎩d

dt

µ∂T e

∂ úu1

¶+∂W e

∂u1= F e1 ,

d

dt

µ∂T e

∂ úu2

¶+∂W e

∂u2= F e2 .

(115)

Observe que cada uma das equações (114) é uma forma quadrática, que expande-se em(W e = 1

2[Ke

11u21 +K

e12u1u2 +K

e21u1u2 +K

e22u

22,

T e = 12[Me

11 úu21 +M

e12 úu1 úu2 +M

e21 úu1 úu2 +M

e22 úu2 úu2.

Fazendo as derivações indicadas em (115) e recolocando o resultado em forma matricial obtemos∙Me11 Me

12

Me21 Me

22

¸ ∙u1(t)u2(t)

¸+

∙Ke11 Ke

12

Ke21 Ke

22

¸ ∙u1(t)u2(t)

¸=

∙F e1 (t)F e2 (t)

¸, (116)

que é exatamente a equação do movimento (103) obtida anteriormente usando o PTV e o princípiode DAlembert para um elemento de barra. Para outros tipos de problemas (vigas, placas, cascas,elementos sólidos etc.) o procedimento é análogo para a determinação das equações de movimento.

7.4.1 Exemplo 2 Equações de Movimento com E.F. de Barra

Determine a equação do movimento para uma barra de comprimento L, área de seção transversal A,densidade ρ e módulo de elasticidade E. Obtenha as matrizes para dois e três elementos Þnitos.

Solução:A Figura 14 ilustra os nós e graus de liberdade do modelo de três elementos. Para dois elementos

a equação do movimento é obtida sobrepondo as matrizes em (106):

ρAL

12

⎡⎣ 2 1 01 4 10 1 2

⎤⎦⎡⎣ u1(t)u2(t)u3(t)

⎤⎦+ 2EAL

⎡⎣ 1 −1 0−1 2 −10 −1 1

⎤⎦⎡⎣ u1(t)u2(t)u3(t)

⎤⎦ =⎡⎣ F1(t)F2(t)F3(t)

⎤⎦ ,e, mesma forma, para três elementos,

ρAL

18

⎡⎢⎢⎣2 1 0 01 4 1 00 1 4 10 0 1 2

⎤⎥⎥⎦⎡⎢⎢⎣u1(t)u2(t)u3(t)u4(t)

⎤⎥⎥⎦+ 3EAL⎡⎢⎢⎣1 −1 0 0−1 2 −1 00 −1 2 −10 0 −1 1

⎤⎥⎥⎦⎡⎢⎢⎣u1(t)u2(t)u3(t)u4(t)

⎤⎥⎥⎦ =⎡⎢⎢⎣F1(t)F2(t)F3(t)F4(t)

⎤⎥⎥⎦ .


y

x

z

E, A, ρ

L

1 3 4

2 3

2

1 2 3

2

F1(t)

(t)2F

F (t)4

F (t)3

u2(t) (t)

3u

L2 = L/3

Figura 14: Modelo de barra com três elementos do Exemplo 2.

8 Análise ModalConsidere o caso em que normalmente se considera como estático, onde o carregamento não variacom o tempo. Devemos lembrar que certamente houve um período inicial em que a carga teve queser aplicada, onde ela teve que variar de zero até seu valor Þnal. Quando este período é suÞcien-temente longo, as acelerações desenvolvidas pela estrutura são baixas o suÞciente para poder seremdesprezadas e a análise pode ser feita como estática, sem o primeiro termo de (107), onde força edeslocamento são agora constantes no tempo, constituindo-se no problema estático de obter o deslo-camento Þnal uf deKuf = F. Isto corresponde, por exemplo, a soltar uma carga sobre a carroceria deum caminhão com inÞnito cuidado. A carroceria baixaria suave e lentamente até atingir sua posiçãoÞnal, como na Figura 15a. Na situação oposta a carga seria simplesmente jogada. A carroceria entãooscilaria várias vezes sobre a suspensão. Devido ao amortecimento, essas oscilações gradualmente sereduziriam enquanto o sistema tenderia a sua posição Þnal de repouso como na Figura 15b.Nota-se então que a classiÞcação de problemas como estático ou dinâmicos nem sempre é simples

e direta. Mesmo que o carregamento varie com o tempo não necessariamente se tem um problemadinâmico. Por exemplo, considere um carregamento cíclico com baixa freqüência. Novamente, se afreqüência de carregamento é baixa, as acelerações, que também são cíclicas, serão baixas. Isto podeser visto de (48). Então as acelerações do sistema podem ser desprezadas em (107), resultando numsistema algébrico dado porKu(t) = F(t). Este é o chamado problema quasi-estático porque, em-bora não tenha o termo de inércia, a resposta varia com o tempo como se fosse um problema estático.Para classiÞcar um problema como quasi-estático ou não basta saber se a freqüência de excitaçãoé baixa o suÞciente. Este pequeno é geralmente quantiÞcado de forma um tanto arbitrária. Se afreqüência de excitação for menor que aproximadamente um terço da menor freqüência natural do


u(t)

u

t to

f

to t

u(t)

(a) (b)

Figura 15: (a) Exemplo de aplicação sem efeitos dinâmicos apreciáveis, (tipicamente para to > 100 sou to > 10Tmax, onde Tmax = 2π/ω1 e ω1 é a menor freqüência natural); (b) Aplicação com respostadinâmica, usuamente para tempos de carregamento da ordem de to < 0, 01 s.

sistema, isto é,

ω . ω13, (117)

então o problema pode ser tratado como quase estático com precisão aceitável.A outra situação é quando as freqüências de carregamento são altas e as forças de inércia devem

ser consideradas, o que constitui o problema da dinâmica. Dois grandes tipos de problemas existem,os problemas de propagação da onda e os de dinâmica estrutural. Os problemas de propagaçãode onda ocorrem em situações de impacto ou de explosões ou de acústica entre outros, onde tanto ocarregamento quanto a resposta são de alta freqüência e o período de duração da análise é em geralcurto, da ordem de um período da onda que cruza a estrutura. Por outro lado, quando a freqüênciade carregamento não é alta, no sentido de que é da mesma ordem, ou apenas algumas vezes maiorque a primeira freqüência natural do sistema, o problema é dito de dinâmica estrutural.A Figura 16 mostra um esboço dos diversos tipos de problemas e análises possíveis, embora

na realidade diversas outras situações existam. Os problemas de dinâmica estrutural, por sua vez,podem ser classiÞcados, pelo menos, em três grandes tipos: (a) determinação de freqüência emodos naturais, (b) análise de resposta temporal e (c) análise de freqüências.As freqüências e modos naturais de uma estrutura são determinados por uma série de motivos.

Numa situação de projeto, frequentemente interessa que a freqüência de carregamento Þque abaixoda primeira freqüência natural, ou pelo menos interessa evitar que a freqüência de excitação Þquepróxima a uma das freqüências do sistema.Na análise da resposta temporal buscamos determinar a resposta do sistema, instante a instante

para um dado histórico de carga. Dois grandes métodos existem para realizar esta análise, (a)análise modal e (b) integração direta. O método de análise modal usa as freqüências e modosnaturais, comentados no parágrafo anterior, enquanto o de integração direta faz uma discretizaçãode diferenciais Þnitos no tempo na equação diferencial do movimento, (107), e faz uma integraçãonumérica. O método de análise modal é um método baseado fundamentalmente na linearidade dosistema. Por outro lado, quando o sistema físico é modelado matematicamente, por elementos Þnitospor exemplo, levando em conta efeitos nãolineares, como plasticidade emmetais, grandes deformações


como em processos de conformação, ou grandes deslocamentos, o processo adequado a ser usado é ode integração direta no tempo das equações de movimento, embora existam formas de circunscreveras limitações da análise modal em alguns casos.Os métodos de análise de freqüências não determinam a resposta do sistema a cada instante,

mas sua composição em freqüências, e não será tratado nesse capítulo. O texto a seguir trataráprimeiramente das freqüências e modos naturais da estrutura que em seguida serão usadas no processode análise modal.

Tipos de coportamento

Quase EstáticoEstático Dinâmico

Propagação de Onda

Dinâmica Estrutural

Resposta Temporal Respostas em Frequência

Análise Modal Integração Direta

Freq. e Modos Naturais

Κu = F Ku(t) = F(t)

Figura 16: ClassiÞcação aproximada do comportamento, tipos de análises e métodos em dinâmica.

8.1 Vibrações Livres Não-Amortecidas

Consideremos F(t) = 0, isto é, o sistema de n equações diferenciais (107) descarregado,

Mu(t) +Ku(t) = 0. (118)

A única coordenada nesse sistema é o tempo, uma vez que as coordenadas espaciais xyz foramjá discretizados. Este tipo de equação é bastante conhecida e estudada em matemática, uma vezque toda uma série de fenômenos físicos é modelada por sistemas de equações diferenciais ordináriasdesse tipo. Uma classe de solução tem a seguinte forma

u(t) = φf(t). (119)


Consideremos o signiÞcado desta função. A Figura 17 ilustra o caso de uma viga modelada portrês elementos Þnitos e quatro nós. A eq.(119) representa os deslocamentos de cada nó ao longo dotempo: ¯

¯ u1(t) = φ1f(t),u2(t) = φ2f(t),u3(t) = φ3f(t),u4(t) = φ4f(t).

(120)

Suponha que num dado instante, t = 0 por exemplo, os deslocamentos sejam φ = [0; 0, 2; 0, 7; 1, 1]T

como na Þgura. Se num certo instante t = t1, f(t1) = f1, e num outro instante, t = t2, f(t2) = f2, porexemplo, isto signiÞca, por esta hipótese que os todos os deslocamentos nodais no instante, t = t2são f2/f1 vezes maior que no instante t = t1. Estes valores nodais são sempre proporcionais entresi a qualquer instante.

1 2 3 4

0.2

0.7

1.1φ2

φ 3

φ4

U1 0

==

=

=

φ =01

f(t = 0) = 1

f(t = 1s) = 4

f(t = 2s) = 0.5

4=4.4φ

0.55φ =4

φ 0.8=2

φ 0.102=

0φ =1 U 2.8=3

=3φ 0.35

Figura 17: Exemplo de deslocamentos nodais proporcionais a um fator comum f(t) que varia notempo.

Observe que (119) não é a solução de (118), mas apenas sua forma geral. Antes desconhecíamosos valores nodais da função do tempo u(t). Agora temos u(t) expresso em termos de um perÞlde deslocamentos nodais φ, independente do tempo, e de um fator comum, f(t), ambos tambémdesconhecidos. A diferenca é que antes tínhamos n funções do tempo a determinar, agora as nincógnitas são constantes, as componentes de φ, e apenas uma função incógnita dependente dotempo, f(t). Substituíndo (119) em (118) temos

Mφ f(t) +Kφ f(t) = 0. (121)

Se consideramos f(t) 6= 0 para qualquer valor de t, podemos dividir tudo por f(t) obtendo:


Mφf(t)

f(t)|z−λ

= −Kφ. (122)

Lembremos que estas são n equações diferenciais em f(t). Mφ é uma matriz coluna de n termos,tanto quanto Kφ. Uma vez que tanto Mφ quanto Kφ são independentes de t, também f(t)/f(t)deve sê-lo. Deve então ser igual a uma certa constante, −λ ainda a ser determinada. Com isto seobtém uma nova equação, em termos apenas de f(t). Se f(t)/f(t) = −λ, de (122) temos os doisproblemas:

f(t) + λf(t) = 0,[K− λM]φ = 0. (123)

1. O primeiro problema é uma equação com forma bastante conhecida, cuja solução tem a forma

f(t) = Aest, (124)

onde s é uma constante a ser determinada. Substituindo em (123)1 obtém-se

As2 est + λAest = 0. (125)

Como A e est não podem ser nulos, eles podem ser simpliÞcados resultando a chamada equaçãocaracterística do problema:

s2 + λ = 0. (126)

Raizes reais Uma primeira solução é obtida supondo-se λ < 0, o que daria duas soluções reais,s1 = −s2 =

√−λ = s. Teríamos duas soluções, f1(t) = Aest e f2(t) = Ae−st, isto é, uma solução(119) crescente exponencialmente no tempo e outra solução decrescente. Mas o carregamento é nuloe o sistema é dito conservativo, isto é, não possui dissipação de energia, amortecimento. A únicaforma deste sistema se mover é simplesmente continuar ummovimento iniciado anteriormente atravésde um impulso rápido aplicado no instante inicial. O movimento deve ser tal que a quantidade deenergia total do sistema deve permanecer constante. Isto não permite que a amplitude do movimentocresça ou diminua ao longo do tempo. Como consequência deve-se ter que λ não pode ser negativo.Como o caso λ = 0 nos remeteria a um problema estático, deve-se então ter λ > 0.Se λ = 0, tem-se um problema estático.Raizes complexas Se λ > 0, faz-se λ = ω2 e as soluções de (126) são

s1 = iω e s2 = −iω,com i =

√−1, e a solução de (123) é uma combinação linear das duas soluções:

f(t) = A1eiωt +A2e

−iωt. (127)

As constantes A1 e A2 são complexas. Pode-se igualmente representar f(t) na forma f(t) =Cei(ωt−φ), e tomar apenas a parte real:

f(t) = C cos (ωt− φ) , (128)


onde C é uma constante arbitrária real, φ o chamado ângulo de fase e ω a freqüência do movi-mento, que é harmônico. Essa freqüência é ainda incógnita, a ser determinada pela solução doproblema (123)2. Uma vez determinados φ e ω de (123), e usando f(t) de (128), a solução doproblema de vibrações livres (118) vem de (119) como:

u(t) = Cφ cos (ωt− φ). (129)

Observe que esse é um movimento oscilante, chamado harmônico. Cada ponto i oscila comfreqüência ω uniforme no tempo, em torno de uma certa posição, com amplitude Cφi.

2. Considera-se agora o problema (123)2

£K− ω2M¤| z

A

φ = 0. (130)

Observe que o fator multiplicando φ é uma matriz A ≡ K − ω2M de ordem n × n. A equaçãomatricial é um sistema algébrico de n equações e n incógnitas, os φi, i = 1, n. De álgebra linearsabe-se que se A for uma matriz quadrada real não-singular, a nulidade do lado direito da equação,(F = 0 em Aφ = F), implica que a única solução possível é φ = 0, isto é, φ1 = φ2 = . . . φn = 0. Aúnica maneira de se ter uma solução não-nula é que K− ω2M seja uma matriz singular, isto é:

detA = 0, isto é, det£K− ω2M¤ = 0. (131)

Como as freqüências são ainda desconhecidas, podemos usar a própria condição (131) paradeterminá-las, bastando que procuremos os valores de ω para os quais o determinante de A sejanulo. Note que o determinante de A é uma função, um polinômio em termos de ω2:

det£K− ω2M¤ = p ¡ω2¢ = 0. (132)

Esse é o chamado polinômio característico, (ou equação de freqüência) associado ao chamadodeterminante característico. Se tivermos um sistema pequeno, 2 × 2 por exemplo, as raízes dopolinômio podem ser obtidas analiticamente. Para sistemas da ordem de milhares, como é comum emelementos Þnitos, existem alguns métodos numéricos e dezenas de variações, que estimam algumasou todas as raízes. Para o caso 2× 2, por exemplo, o problema (131) Þca

p (λ) = det

∙K11 − λM11 K12 − λM12

K21 − λM21 K22 − λM22

¸= 0,

p (λ) = (K11 − λM11) (K22 − λM22)− (K21 − λM21) (K12 − λM12) = 0. (133)

Observe que esse é um polinômio do 2o. grau em λ. De forma geral, a algébra mostra que umsistema de n equações reais possui um polinômio característico de grau n e possui n raízes. Oproblema (130) é denominado problema de autovalor ou autoproblema, enquanto as raízes dopolinômio característico, os ω0is, são chamados autovalores do problema. Calculamos então os nautovalores λj = ω2j do problema. A cada autovalor substituído em (130) poderemos resolver e obterum distinto vetor solução φj, isto é£

K− ω2jM¤φj = 0 j = 1, 2, · · · , n. (134)


Cada vetor solução φj é chamado autovetor do problema, e o par¡ωj;φ

j¢é um autopar. φj

é também denominado modo de vibração do sistema. Uma vez que não temos apenas um parde solução do problema de autovalor (130), mas n autopares, a solução do problema dinâmico (118)não é apenas (130), mas uma combinação linear de todos os modos na forma:

u(t) = C1φ1 cos (ω1t− φ1) + C2φ2 cos (ω2t− φ2) + . . .+ Cnφn cos (ωnt− φn) ,

isto é,

u(t) =Pn

j=1Cjφj cos

¡ωjt− φj

¢. (135)

As várias constantes Cj devem ser determinadas de acordo com as condições iniciais do sistema,como sera visto posteriormente.

8.2 Propriedades dos Autovetores e Autovalores

Nos próximos itens exploraremos as características, usos e signiÞcados físicos das freqüências e modosnaturais de um sistema. Antes disso porém, vamos tratá-los simplesmente como entidades matemáti-cas, autovalores e autovetores, e observar suas propriedades.

8.2.1 Ortogonalidade

A primeira propriedade a ser demonstada é a seguinte: considere dois distintos autopares de (130),isto é, (ωr;φ

r) e (ωs;φs), que satisfazem

Kφr = ω2rMφr e Kφs = ω2sMφ

s. (136)

Se multiplicarmos a primeira equação pelo transposto de φs, isto é, φsT , e a segunda por φrT

obtemos2

φsTKφr = ω2rφstMφr e φrTKφs = ω2sφ

rtMφs. (137)

Observe que, enquanto a equação (136)1 consiste de uma igualdade entre dois vetores, isto é, Kφr

é igual a um certo vetor Vr eMφr é igual a um certo vetor Ur. De forma similar para a eq. (136)2.Quando pré-multiplicamos (136)1 por um autovetor φ

sT , isto equivale a um realizar produto escalarφsT Vr, cujo resultado é um escalar. Podemos transpor uma das duas equações (137), a segundapor exemplo, e o sistema Þca

φsTKφr = ω2rφsTMφr e φsTKTφr = ω2sφ

sTMTφr. (138)

Como K eM são matrizes simétricas, os termos se tornam idênticos entre as duas equações. Sesubtrairmos a primeira da segunda equação temos

0 =¡ω2r − ω2s

¢φsTMφr| z

a

. (139)

φsTMφr é um escalar a. Se as freqüências naturais forem distintas, ωr 6= ωs, é então necessário que

φsTMφr = 0 para qualquer r 6= s se ωr 6= ωs. (140)

2O sobre-índice T indica transposto de um vetor ou matriz.


Este resultado é chamado de condição de ortogonalidade dos vetores modais. Caso ωr = ωs,prova-se que o correspondente par de autovetores φr e φs não necessariamente são ortogonais. En-tretanto,esses vetores podem ser ortogonalizados, usando, por exemplo,o método de ortogonalizaçãode Gran-Schmidt.

Assim,deforma geral, considera-se sempre que se tem o conjuntodos n autovetores ortogonais pela massa.

A operação φsTMφr pode ser vista como um tipo diferente de produto escalar entre os vetoresφs e φr. Esta tipo de produto escalar entre dois vetores U e V é deÞnido com o uso de umamatrizpeso, no caso é usado a matriz massa, de forma que U · V ≡ UTMV, enquanto a forma maisconhecida do produto interno é o chamado produto euclidiano, dado por U ·V ≡ UTV. Observeque se dois distintos vetores modais φr e φs são M-ortogonais, isto é, satisfazem (140), eles sãotambém K-ortogonais. A demonstração é feita simplesmente levando (140) para o lado direito de(138)1, o que resulta em:

φsTKφr = 0, (141)

isto é, se dois vetores são ortogonais em relação à massa também o serão em relação à rigidez.

8.2.2 Normalização e Ortonormalidade

Se temos deÞnido um produto escalar, também chamado produto interno entre dois vetores, temosentão uma deÞnição de comprimento, ou norma kφrk de um vetor φr, deÞnida por

kφrk =pφr φr,

=pφrTMφr

(142)

É visível que se o termo dentro do radical puder ser negativo para algum vetor φr, a deÞniçãoperde o signiÞcado, uma vez que não se poderia interpretar como comprimento um valor negativo.Então esta norma só pode ser deÞnida com uma matriz peso que tenha a propriedade de ser positivadeÞnida. Uma matriz A é dita positiva deÞnida se°°°° VTAV ≥ 0, para qualquer V, e

VTAV = 0⇐⇒ V = 0,(143)

isto é, algumas matrizes sempre terão o resultado VTAV positivo, qualquer que seja o vetor con-siderado, exceto no caso dele ser nulo. Agora observe novamente o autoproblema (130). Suponhaque já encontramos um autopar (ωr;φ

r) do problema, que, como tal, satisfaz:

[K− ωrM]φr = 0.Podemos multiplicar estas n equações por uma constante d qualquer e colocar a equação na forma

[K− ωrM] (cφr) = 0.Concluímos que se φr é um autovetor, dφr também o será. Assim, após a determinação de cadaautovetor fazemos sua normalização, isto é, calculamos sua norma por (142) e fazemos

φr =1

kφrkφr. (144)


Visivelmente, agora φr tem norma unitária, isto é, φrTMφr = 1. Como os diversos autovetoressão ortogonais entre si, podemos escrever a seguinte relação

φrTMφs = δrs, r, s = 1, 2, 3, ..., n (145)

onde δrs é um símbolo conhecido como delta de Kronecker. Sua deÞnição é a de que δrs = 1 ser = s, isto é, se tivemros r = s = 1 ou 2, etc. Por outro lado, se r 6= s, (por exemplo r = 1 e s = 2),por deÞnição tem-se que δrs = 0 . Resumindo,

δrs = 1 se r = s,= 0 se r 6= 0. (146)

Um conjunto de vetores que possui a propriedade mostrada em (145) é dito um conjuntoortonormal de vetores, isto é, cada um é normalizado para norma unitária e também é ortogonala todos os demais. A relação (145) está colocada na chamada forma indicial, com os índices r e spodendo assumir valores entre 1 e n. Essa relação pode também ser colocada numa forma matricialcompleta. Para isso deÞne-se a chamada matriz modal como:

Φ =£φ1φ2 . . . φn

¤, (147)

isto é, Φ é a matriz n × n em que cada coluna é composta por um dos autovetores do problema.Assim, a relação de ortonormalidade (145) pode ser colocada na forma matricial

ΦTMΦ = I (148)

onde I é uma matriz indentidade de ordem n × n. Consideremos novamente o autoproblema. Emvez de representar um autopar de solução a cada vez, como em (136)1, aplicamos todos os autoparessimultaneamente. Isto é feito da seguinte forma:

KΦ =MΦΛ2 (149)

onde Λ2 é uma matriz diagonal, denominada matriz espectral, composta pelos autovalores:

Λ2 =

⎡⎢⎢⎢⎣ω21

ω22. . .

ω2n

⎤⎥⎥⎥⎦ . (150)

De forma expandida, (149) pode ser posto como

⎡⎢⎢⎢⎣K11 K12

K21 K22

. . .Knn

⎤⎥⎥⎥⎦⎡⎢⎢⎢⎣φ11 φ21 · · · φn1φ12 φ22...

.... . .

...φ1n φ2n φnn

⎤⎥⎥⎥⎦

=

⎡⎢⎢⎢⎣M11 M12

M21 M22

. . .Mnn

⎤⎥⎥⎥⎦⎡⎢⎢⎢⎣φ11 φ21 · · · φn1φ12 φ22...

.... . .

...φ1n φ2n φnn

⎤⎥⎥⎥⎦⎡⎢⎢⎢⎣ω21

ω22. . .

ω2n

⎤⎥⎥⎥⎦ .


Note que (149) representam as n equações associadas a cada autovetor. Então tem-se de faton× n equações algébricas. Podemos em seguida pré-multiplicar (149) por ΦT obtendo:

ΦTKΦ = ΦtMΦΛ2.

Mas com a ortonormalidade dos autovetores em relação à matriz massa, eq. (148), o lado direitosimpliÞca-se e temos

ΦTKΦ =Λ2. (151)

Já tínhamos visto em (141) a ortogonalidade dos vetores em relação à rigidez, isto é, φrKφs = 0.Agora temos também que φrTKφr = ω2r, isto é, a norma de um autovetor em relação a matriz derigidez é o quadrado do correspondente autovalor. A seção ?? descreve também outras propriedadesdas matrizes e autopares do problema.

8.2.3 Exemplo 3 Freqüências Naturais

Considere a barra do Exemplo 2, engastada na extremidade esquerda, modelada por três elementos.Obtenha a aproximação de elementos Þnitos para sua primeira e segunda freqüência natural. UseE = 200.000MPa, ρ = 7.800kg/m3, A = 1, 0 cm2 e L = 1, 0 m.

Solução:As freqüências naturais são as raizes ω2j do polinômio característico deÞnido em (132) pelo deter-

minante det[K − ω2M] = 0. Da solução do Exemplo 2, o problema de autovalor para um modo jé: ⎧⎪⎪⎨⎪⎪⎩

3EA

L

⎡⎢⎢⎣1 −1 0 0−1 2 −1 00 −1 2 −10 0 −1 1

⎤⎥⎥⎦− αj⎡⎢⎢⎣2 1 0 01 4 1 00 1 4 10 0 1 2

⎤⎥⎥⎦⎫⎪⎪⎬⎪⎪⎭⎡⎢⎢⎣φj1φj2φj3φj4

⎤⎥⎥⎦ =⎡⎢⎢⎣0000

⎤⎥⎥⎦ . (152)

onde αj = ω2jρAL18.

Deve-se primeiramente aplicar as condições de contorno para vincular a barra. Uma vez que elaestá engastada pelo nó 1, qualquer que seja seu movimento vibratório este deve ser tal que u1(t) = 0.Então todos os modos de vibração devem ser tais que φj1 = 0. Levando este valor à equação signiÞcaeliminar a primeira coluna de cada matriz junto com o termo φj1. Em seguida eliminamos a primeiralinha, Þcando então com matrizes 3× 3. Quanto às constantes multiplicativas, dividimos a equaçãopor 3EA/L e deÞnimos

αj = ω2j

ρAL/18

3EA/L= ω2j

ρL2

54E(153)

O polinômio característico então Þca

p(αj) = det

⎧⎨⎩⎡⎣ 2 −1 0−1 2 −10 −1 1

⎤⎦− αj⎡⎣ 4 1 01 4 10 1 2

⎤⎦⎫⎬⎭ = 0, (154)

que pode ser simpliÞcado para:


Tabela 1: Freqüências naturais em modelo de barra em balanço.Freqüência em Hz (Erro %)

Freq. Analítico (Hz) 1 elemento 2 elementos 3 elementosω1

14L

pE/ρ = 1.265, 9 1.395,8 (10,3 %) 1.298,3 (2,56 %) 1.280,4 (1,10 %)

ω234L

pE/ρ = 3.797, 8 4.536,5 (19,5 %) 4.187,6 (10,3 %)

ω354L

pE/ρ = 6.329, 5 7.597,0 (16,7 %)

p(αj) = (2− 4αj)2(1− 2αj)− (1 + αj)2(1− 2αj)− (1 + αj)2(1− 4αj) = 0.As três raizes são

α1 =11− 6√3

13= 0, 0467458m−2s−4 −→ ω1 = 1.280, 4Hz,

α2 =1

2= 0, 0467458m−2s−4 −→ ω2 = 4.187, 6Hz, (155)

α3 =11 + 6

√3

13= 1, 64556m−2s−4 −→ ω3 = 7.597, 0Hz.

Observe que usando dois elementos as freqüências aproximadas são: ω1 = 1.298, 3Hz e ω2 =4.536, 5Hz, enquanto usando um único elemento a primeira freqüência é aproximada por ω1 =1.395, 8Hz. A comparação dos resultados com a solução analítica é vista na Tabela 1.Observe que os modos iniciais convergem primeiro e os modos mais altos sempre

requerem malha mais reÞnada para atingir precisões aceitaveis. Isto é regra geral nasaproximações por elementos Þnitos.

8.2.4 Exemplo 4 Modos Naturais

Considere a barra engastada do Exemplo 2 modelada por três elementos Þnitos. Determine os modosnaturais de vibração da barra. As matrizes do problema são (a deÞnição de αj é dada em (153)):⎧⎨⎩

⎡⎣ 2 −1 0−1 2 −10 −1 1

⎤⎦− αj⎡⎣ 4 1 01 4 10 1 2

⎤⎦⎫⎬⎭⎡⎣ φj2φj3φj4

⎤⎦ =⎡⎣ 000

⎤⎦ . (156)

Solução:O autovetor φj é obtido substituindo o valor de ωj da Tabela 1 (Exemplo 3) em (152) e resolvendo

o sistema para cada modo j. Para omodo 1, usamos ω1 = 1.280,4 Hz = 8.045 s−1, o que correspondea α1 = 0,04674m−2s−4. A eq.(156) para o modo j = 1 Þca⎡⎣ 1, 81302 −1, 04674 0

−1, 04674 1, 81302 −1, 046740 −1, 04674 0, 90651

⎤⎦⎡⎣ φ12φ13φ14

⎤⎦ =⎡⎣ 000

⎤⎦ .Triangularizando a matriz por fatorização gaussiana temos:


⎡⎣ 1, 81302 −1, 04674 00 1, 20868 −1, 046740 0 0

⎤⎦⎡⎣ φ12φ13φ14

⎤⎦ =⎡⎣ 000

⎤⎦ .Podemos fazer φ14 = 1, 0. Neste caso resolvemos φ

12 = 1, 5 e φ

13 = 1, 15471, isto é,

φ1 =

⎡⎣ 0, 50, 8661, 0

⎤⎦ .Para o modo 2, o autovalor é α2 = 0, 0467458m−2s−4. Os autovetores são

φ2 =

⎡⎣ 1, 00

−1, 0

⎤⎦ e φ3 =

⎡⎣ 0, 5−0, 8661, 0

⎤⎦ .

8.2.5 Exemplo 5 Normalização de Autovetores

Considere o autoproblema do Exemplo 4. Normalize os autovetores, forme as matrizes massa erigidez.

Solução:Os autopares obtidos foram:

(α1;φ1) =

⎛⎝11− 6√32

;

⎡⎣ 0, 5√3/21, 0

⎤⎦⎞⎠ , (α2;φ2) =

⎛⎝0, 5;⎡⎣ 10−1

⎤⎦⎞⎠ ,(α3;φ

3) =

⎛⎝11 + 6√32

;

⎡⎣ 0, 5

−√3/21, 0

⎤⎦⎞⎠ .Formamos a matriz modal não-normalizada

Φ =

⎡⎣ 1/2 1 1/2√3/2 0 −√3/21 −1 1

⎤⎦ .O teste de ortononalidade com a massa é feito por

A = ΦTMΦ =ρAL

18ΦT

⎡⎣ 4 1 01 4 10 1 2

⎤⎦Φ = ρAL

18

⎡⎣ 8, 59808 0 00 6, 0000 00 0 3, 40192

⎤⎦ .Essa é uma matriz diagonal, o que conÞrma a ortogonalidade dos autovetores. Os termos da

diagonal são os quadrados das normas dos autovetores, isto é, Ajj =°°φj°°2 = φjTMφj. Então,

pode-se obter a matriz modal normalizada dividindo cada coluna j de A pela norma°°φj°°. Isto

resulta em:


Φ =

r18

ρAL

⎡⎣ 0, 170518 1/√6 0, 271037

0, 295345 0 −0, 4695360, 341035 −1/√6 0, 542173

⎤⎦ . (157)

Observe que

ΦTKΦ =3EA/L

ρAL/18

⎡⎣ 0, 046746 0 00 0, 5 00 0 1, 64556

⎤⎦ = Λ2,que é a matriz com as freqüências naturais da barra, como pode ser visto nas eqs.(153), (155) e (150).

8.2.6 Exemplo 6 Solução Analítica de Vibrações Livres em Barra

Considere o problema de vibrações livres de uma barra engastada numa das extremidades como noExemplo 2. a) Encontre a solução analítica para o problema; b) Plote a solução analítica e a deelementos Þnitos para os primeiros modos.

Solução:a) Para a solução analítica, primeiro usamos a equação diferencial do problema estático, eq.(??):

AEd2u(x)

dx2+Abx(x) = 0 ∀x ∈ (0, L), (158)

e aplicamos o princípio de DAlembert, tomando a força de corpo bx como sendo a força de inércia,−ρu. A equação do movimento para vibrações livres, e condições iniciais e de contorno Þcam então:¯

¯¯AE

d2u(x, t)

dx2− ρA d

2u(x, t)

dt2= 0, ∀x ∈ (0, L), t > 0,

u(0, t) = 0, t > 0,

AEdu

dx(L, t) = Afx = 0, t > 0.

(159)

(fx é a força na extremidade x = L da barra.) Esse problema pode ser resolvido pelo método deseparação de variáveis, cujo ponto de partida consiste em supor que a solução pode ser escrita naforma u(x, t) = X(x)T (t), que substituida em (159)1 produz:

X 00(x)T (t)− ρ

EX(x)T (t) = 0.

Pode-se separar os termos dependentes de x e t da seguinte forma:

ρ

E

X 00(x)X(x)

=T (t)

T (t)= ω2. (160)

Uma vez que lado esquerdo da igualdade é função apenas de x enquanto o lado direito dependeapenas de t, conclui-se que cada lado deve ser constante. Denominemos essa constante de ω2. Istogera duas equações diferenciais ordinárias:


⎧⎨⎩ X 00(x) + ω2ρ

EX(x) = 0,

T (t) + ω2T (t) = 0.

(161)

A primeira dessas equações, quando juntadas às condições de contorno, gera o primeiro prob-lema de valor no contorno. As condições de contorno são de que uma das extremidades éengastada, u(0, t) = 0, e a outra é livre de forças. Isto resulta no problema°°°°°°°°°°°°°°

Ed2X(x)

dx2− ω2ρ|z

Eβ2

X(x) = 0,

X(0) = 0,

dX

dx

¯x=L

= 0.

(162)

Esse é um problema de autovalor contínuo, não discretizado com aqueles vistos anteriormente.Para compactar a notação, deÞne-se β2 = ω2ρ/E. Essa equação diferencial tem solução bastanteconhecida. Se ω for conhecido, a solução é:

X(x) = A senβx+B cosβx. (163)

Resolvendo para a primeira condição de contorno obtemos B = 0. Para a segunda condição,

dX

dx

¯x=L

= βA cosβx = 0.

Como βA 6= 0, é necessário que βL = (2n−1)π/2, para n = 1, 2, 3, ...Estes são então os autovaloresda barra:

βn =(2n− 1)π2L

, n = 1, 2, 3, ...,

ou, da deÞnição de β, as freqüências naturais da barra engastada são:

ωn =(2n− 1)π2L

rE

ρ, n = 1, 2, 3, ..., (164)

Os correspondentes modos naturais de vibrações são obtidos levando as freqüências para (163) comB = 0:

Xn(x) = An senβnx, n = 1, 2, 3, ... (165)

b) A solução aproximada para os modos de vibrações obtida por três elementos é vista na eq.(157),normalizada pela massa. Esses são valores nodais da solução. Já os valores de deslocamentos empontos arbitrários x da barra são obtidos com o uso das funções de interpolação. Uma vez queestas funções são lineares em x ao longo de cada elemento, obtém-se a distribuição de deslocamentosde forma quebrada, como mostrada na Figura 18. Os valores de βn na solução analítica sãoβ1 = 1, 57, β2 = 4, 71 e β3 = 7, 85. A solução analítica normalizada de forma a que o módulo do


0.00 0.25 0.50 0.75 1.00x

-0.50

-0.25

0.00

0.25

0.50

modo 2

modo1

modo 3

φ(x),

L

X(x)

Figura 18: Solução analítica Xn(x) e aproximada φn(x), obtida por três elementos Þnitos de barra,para os três primeiros modos naturais de vibrações.

deslocamento Xj(L) coincida com o valor de φj em x = L, isto é, em (165) a constante Aj foi obtida

por Aj = φj4/senβjL. Então, A1 = 0, 3410, A2 = −0, 4082 e A3 = 0, 5422.

Observe que a precisão obtida decai progressivamente para os modos mais altos, tal como asfreqüências naturais, vists no Exemplo 2. Isto é regra no método de elementos Þnitos. Nota-se quetambém nesse exemplo que os valores nodais foram sempre exatos. Isto ocorre sempre, masapenas em problemas unidimensionais. Problemas bi- ou tridimensionais como os de placa ousólidos apresentam valores nodais apenas aproximativos da solução analítica.Note que a plotagem da Figura 18 representa, na realidade, deslocamentos axiais, e não transver-

sais,como a Þgura pode sugerir.

8.2.7 Autovetores Linearmente Independentes

Monstraremos que o conjunto dos n autovetores associados às matrizes massa e rigidez do sistemamecânico formam um conjunto de vetores linearmente independentes. Um conjunto linearmenteindependente signiÞca que qualquer um de seus elementos não pode ser escrito como uma combinaçãolinear dos demais. Para a prova consideraremos o contrário, isto é, que o conjunto é linearmentedependente, tal que, se tomarmos, por exemplo, o primeiro vetor, poderiamos escreve-lo como umacombinação linear dos demais:

φ1 = a2φ2 + a3φ

3 + . . . anφn, (166)

com constantes a2, . . . , an não-nulas. Passando os termos da direita para a esquerda e mudando osnomes das constantes tem-se

b1φ1 + b2φ

2 + . . .+ bnφn = 0.


Isto pode ser colocado em forma matricial como Φb = 0, onde b = [b1b2 . . . bn]T . Pré-multiplicando

por ΦTM tem-se

ΦTMΦ| z I

b = 0 =⇒ b = 0,

isto é, b1 = b2 = . . . = bn = 0, devido a ortonormalidade dos vetores. Então (166) só podeser satisfeito se todos os ai = 0. Mas b1 = 1, logo (166) não é possível, isto é, φ1 não pode serrepresentado como uma combinação linear dos demais vetores. Repetindo o procedimento para osdemais vetores, conclui-se que o conjunto de autovetores Φ é linearmente independente.Da teoria de álgebra linear tem-se então a seguinte consequência: qualquer vetor φ de ordem n

pode ser escrito como uma combinação linear dos vetores modais φr, isto é, qualquer φ pode serrepresentado por:

φ = d1φ1 + d2φ

2 + . . .+ dnφn =

nXj=1

djφj,

ou ainda, numa forma matricial, o teorema de expansão pode ser representado por:

φ = Φd, (167)

onde d é o vetor coluna composto pelos n coeÞciente dj tal que d = [d1, d2, . . . dn]T .

O conjunto Φ de autovetores é então dito ser uma base do espaço de dimensão Þnita n-dimensional gerado pelas matrizes massa e rigidez. A expressão (167) caracteriza o chamadoteorema de expansão. Ele é fundamental ao método de análise modal visto a seguir.

8.3 Análise Modal para Excitação Inicial - Sistema não-Amortecido

A aplicação do método de análise modal mais simples é no caso em que a única excitação no sistemaé aquela aplicada no instante inicial, e ele permanece em movimento permanentedevido à ausênciade dissipação de energia. Ao longo do tempo o carregamento é nulo, isto é, F(t) = 0 para t > 0, e aequação de movimento discretizada (107) reduz-se à eq.(118):

Mu(t) +Ku(t) = 0. (168)

Primeiramente formamos o problema de autovalor (130) e obtemos sua solução, os n auto-valoresΛ2 e os correspondentes autovetores Φ. Usamos então o teorema da expansão (167), isto é,se qualquer vetor no espaço n-dimensional gerado por K e M pode ser expandido comouma combinação linear dos modos naturais, então a solução u(t) de (168) também devepoder, já que este vetor também pertence ao mesmo espaço. Então pode-se representar asolução u(t) por:

u(t) = Φη(t), (169)

onde η(t) = [η1(t), η2(t) . . . ηn(t)]T são os coeÞcientes, os ds de (167). u(t) são os deslocamentos

nodais físicos do sistema, enquanto η(t) são coordenadas generalizadas, aqui denominadas tam-bém coordenadas modais. Nesse caso, como u(t) é um vetor que muda a cada instante, tambémos coeÞcientes ηj(t) devem ser função do tempo. Aparentemente não se conseguiu nenhum beneÞciocom a transformação vetorial acima. Simplesmente passamos de n funções incógnitas uj(t) paraoutras n funções incógnitas ηj(t). Mas a seguir se poderá visualizar sua utilidade.


Realizamos então as seguintes operações: substituímos (169) em (168) e pré-multiplicamos estaúltima por ΦT , o que resulta

ΦtM Φ η(t)+Φt K Φ η(t) = 0.

Se o conjunto Φ de autovetores tiver sido ortonormalizado, com o uso de (148) e (151) a equaçãoacima Þca

Iη(t) + Λ2η(t) = 0 (170)

Observe que, comoΛ2 é umamatriz diagonal, este é um conjunto de n equações diferenciais ordinárias,homogêneas, com coeÞcientes constantes, desacoplado, isto é, a j-ésima equação tem a forma

ηj(t) + ω2j ηj(t) = 0, (171)

ou ainda, de forma extendida, as n equações são:¯¯¯η1(t) + ω

21η1(t) = 0,

η2(t) + ω22η2(t) = 0,...

ηn(t) + ω2nηn(t) = 0.

Compare (168) a (170). As eqs.(168) são também n equações diferenciais, mas cada uma delas,a equação j-ésima, por exemplo, tem a forma

Mj1u1 +Mj2u2 + . . .+Mjnun +Kj1u1 +Kj2u2 + . . .+Kjnun = 0. (172)

Envolve portanto todas as n funções incógnitas uk(t), sendo um sistema do tipo acoplado. Já em(170) a j-ésima equação envolve apenas uma única função incógnita, ηj(t). Desta forma, cada umadas equações (170) pode ser resolvida separadamente. Quando todos os ηj(t) forem determinados,usamos a transformação (169) e temos a solução u(t). Note que os coeÞcientes ηj(t) são apenascoeÞcientes, funções temporárias no processo do cálculo, não tendo signiÞcado físico deÞnido,como tem u(t).Passemos aos detalhes. Cada equação desacoplada em (170) tem solução já determinada, uma

vez que corresponde ao problema de vibrações livres não-amortecidas de um grau de liberdade. Asolução pode ser tomada de (29) para amortecimento nulo, ξ = 0:

ηj(t) = Aj cos(ωjt− φj), (173)

onde Aj e φj são constantes a serem determinadas de forma que a função satisfaça as condiçõesiniciais do sistema. Considere o sistema mostrado na Figura 19. Ali se mostra uma das muitasmaneiras de modelar matematicamente o comportamento dinâmico de um veículo. No caso, podem-se obter informações úteis por uma simulação simpliÞcada, onde o veículo é modelado como umasimples viga apoiada sobre duas molas que representam toda a ßexibilidade dos pneus e do sistemade suspensão. Aplicado um carregamento impulsivo sobre o veículo, cada parte dele se põe a mover.No esquema da Figura 19b indicamos uma modelagem de elementos Þnitos de cinco elementos deviga para modelar o comportamento do modelo inicial. Note que temos dois modelos: o veículo físicofoi primeiro modelado pela viga e molas. O segundo modelo, o de elementos Þnitos visa aproximar ocomportamento dinâmico do primeiro modelo, que por sua vez pretende-se que dê informações sobreo comportamento do veículo em si.


1 2 3 4 5 6

Ka

u (t)1

bK

Figura 19: (a) Esquema de um veículo, e (b) uma possível representação esquemática.

Buscamos então as seis funções, u1(t) a u6(t) em cada um dos seis nós da malha da Figura 19b,por exemplo. Para a resolução de qualquer problema diferencial que envolve o tempo, é necessárioque se tenha um problema de valor inicial, dito bem posto. Isto signiÞca que devemos tera(s) equações diferenciais que regem o problemas. Mas elas não são suÞcientes. É necessário queconheçamos as condições de contorno, que em geral dizem como o sistema está vinculado. Alémdisso, é necessário que se conheça as condições iniciais do sistema, isto é, devemos conhecer osdeslocamentos e velocidades de cada ponto do sistema no instante inicial, t = 0. No exemplo daFigura 19b isto signiÞca conhecer os valores de deslocamento em cada um dos nós no instante t = 0.Frequentemente a origem da medida de deslocamentos é escolhida como a conÞguração do sistemano instante inicial, ou em sua conÞguração de equilíbrio. Nestes casos as condições iniciais seriam

u1(0) = u2(0) = . . . = un(0) = 0. (174)

Observemos também que no instante inicial cada ponto do sistema pode estar se movendo comuma certa velocidade. Estes valores de deslocamentos e velocidades iniciais, quando não nulos, fazemàs vezes de excitação no sistema. Observe o caminhão do exemplo. Considere que ele está na Figura19a, em sua posição de equilíbrio estático, imovel. Se nada for feito com ele, ele continuará paradocomo é previsto pela segunda lei de Newton. Agora se aplica uma força vertical sobre ele, mas semque se faça medição ou se conheça esta força. Ampliamos o valor da força até que o deslocamentomedido em cada um dos seis nós atinja certos valores, por exemplo, u01, u02, u0n. Neste momentoremovemos a carga e começamos a contar o tempo. Podemos então considerar que o sistema possuiuma distribuição de deslocamentos iniciais ui(0) = uoi, ou, em forma vetorial,

u(0) = uo. (175)

É de se esperar que o sistema continue a mover-se, devido a interações entre as forças de inércia eas elásticas. Isto é o que explica por que a equação do movimento (168) pode ter solução não nula senela não aparece carregamento. A excitação aparece em termos de deslocamentos iniciais. O mesmopode ser dito sobre as velocidade iniciais. Podemos ter a cada nó valores conhecidos de velocidadeinicial:

úu(0) = vo, (176)

isto é, úu1(0) = υ01, úu2(0) = υ02, . . . , úun(0) = υ0n.


Note que, com uso da transformação (169), transformamos o conjunto de equações diferenciais demovimento (168), escrito em termos de deslocamentos, nas equações diferenciais (170), escritos emtermos das funções nj(t). Como desejamos resolver primeiro (170), devemos transformar tambémas condições iniciais (175) e (176) em termos de nj(t). Aplicamos então a transformação (169) àscondições iniciais:

u(0) = Φ η(0) ⇒ u0 = Φ η0,

úu(0) = Φ úη(0) ⇒ v0 = Φ úη0.(177)

Pré-multiplicamos as equações do lado direito pela matriz ΦTM obtendo

ΦTMuo = ΦTMΦηo,

ΦTMvo = ΦTMΦ úηo.

Devido à ortonormalidade dos autovetores em Φ temos

ηo = ΦTMuo,

úηo = ΦTMvo,

(178)

o que nos dá os valores iniciais das funções η(t), a serem usados na solução do problema (170).Aplicamos (172) no instante inicial:

ηj(0) = Aj cosφj,

úηj(0) = Ajωj senφj.(179)

Observe que cada ηj(0) e úηj(0) é já conhecido de (178). Então (179) forma um sistema de duasequações para cada j,e duas incógnitas, Aj e φj. Se dividimos a segunda pela primeira equação temos

tanφj =úηj(0)

ωjη(0). (180)

Obtendo φj, de (179)1 obtemos Aj como:

Aj =ηj(0)

cosφj. (181)

Estas constantes são calculadas para cada uma das equações desacopladas, de forma que todas asfunções ηj(t) são conhecidas de (173). Conhecidas estas funções a solução dos deslocamentos nodaisdo sistema não amortecido sob excitação inicial vem de (169), isto é,

u(t) =nPj=1

φjAj cos¡ωjt− φj

¢. (182)

Esta solução pode também ser calculada numa forma alternativa. Usamos a solução (37) doproblema de vibração livre unidimensional. Fazendo o amortecimento nulo naquelas equações, ζ = 0,obtemos as simpliÞcações ωd = ωn e φ = 0, o que resulta na solução do problema desacoplado como:

ηj(t) = ηj(0) cosωjt+úηj(0)

ωjsenωjt, para j = 1, . . . , n, (183)


onde ηj(0) e úηj(0) são obtidos resolvendo (178). Então os deslocamentos físicos vem de (169) como:

u(t) =nPj=1

φj∙ηj(0) cosωjt+

úηj(0)

ωjsenωjt

¸. (184)

8.3.1 Exemplo 7 Resposta para Deslocamento Inicial pelo MEF

Considere uma barra como a do Exemplo 2, Figura 14 na página 28, de comprimento L = 1, 0m, áreaA = 1, 0 cm2, densidade ρ = 7.800 kg/m3 e módulo E = 2 ·105MPa. Inicialmente se aplica uma forçaaxial F = 105N como na Figura 20. Esta força é aplicada de forma quasi-estática, isto é, lentamente,sem gerar acelerações apreciaveis, até que se atinja o deslocamento máximo na barra. Neste instantea força é removida subiamente. O sistema então começa vibrar axialmente. Determine esta respostadinâmica do sistema. Discretize a barra em três elementos Þnitos.

1 2 3 4

L , E , A , ρF4

Figura 20: Barra discretizada por três elementos, sob deslocamento inicial provocado pela remoçãoda carga F .

Solução:Após o instante em que a carga é liberada temos um problema de vibração livre com condição

de deslocamento inicial prescrito não nulo. Primeiramente então devemos calcular os deslocamentosiniciais nodais uo.No Exemplo 2 na página 28, tinhamos já obtido as matrizes K eM para uma modelagem de três

elementos, e do Exemplo 5 temos as freqüências e modos naturais de vibração na eq.(157). Essesvalores são os seguintes:

K =ρAL

12

⎡⎣ 2 1 01 4 10 1 2

⎤⎦ ; M =2EA

L

⎡⎣ 1 −1 0−1 2 −10 −1 1

⎤⎦ ;Φ =

r18

ρAL

⎡⎣ 0, 1705 1/√6 0, 2710

0, 2953 0 −0, 46950, 3410 −1/√6 0, 5422

⎤⎦ , (185)

e as freqüências são:

[ω1; ω2; ω3] = [8045, 0; 26311, 5; 47733, 4] s−1. (186)

O problema de determinação de uo é um problema estático, deÞnido por Kuo = F, isto é,


ρAL

12

⎡⎣ 2 1 01 4 10 1 2

⎤⎦⎡⎣ uo1uo2uo3

⎤⎦ =⎡⎣ 0

0105

⎤⎦N.A solução é: uo =

£1; 2; 3

¤m/600 =

£1, 667; 3, 333; 5, 0

¤m. As condições de velocidade

inicial são vo = 0. Passamos a seguir à determinação do sistema desacoplado de equações domovimento, eq.(170). Como a matriz modal já está normalizada pela massa, as eqs.(185) resultamem ⎡⎣ η2(t)η3(t)

η4(t)

⎤⎦+⎡⎣ 6, 47 · 107 6, 92 · 108

2, 28 · 109

⎤⎦⎡⎣ η2(t)η3(t)η4(t)

⎤⎦ =⎡⎣ 000

⎤⎦ . (187)

A determinação das condições iniciais ηo e úηo é feita por (178), o que resulta em:

ηo = ΦTMuo =

⎡⎣ 25, 31−2, 8321, 146

⎤⎦ 10−4, úηo = ΦTMvo = 0. (188)

De (180) os ângulos de fase em cada modo são: φ1; φ2; φ3 = 0; 0; 0, e de (181) as correspondentesamplitudes são:

Aj =ηj(0)

cosφj −−−−→ A = 25, 31;−2, 832; 1, 14610−4.

A resposta do sistema é dada por (182) que Þca na forma:

u =3Xj=1

φjAj cosωjt,

= φ1A1 cosω1t+ φ2A2 cosω2t+ φ

3A3 cosω3t.

isto é,

⎡⎣ u2(t)u3(t)u4(t)

⎤⎦ =

⎡⎣ 2, 0733, 5904, 146

⎤⎦ 10−3 cos(8045, 0t) +⎡⎣ −101

⎤⎦ 5, 554 · 10−4 cos(26331, 5t)+

⎡⎣ 1, 492−2, 5842, 984

⎤⎦ 10−4 cos(47733, 4t). (189)

Observe que no instante inicial, t = 0, esta solução dá:

u(0) =

⎡⎣ 1, 6673, 3335, 000

⎤⎦ 10−3 = uo, (190)

como era de se esperar.


0.0000 0.0020 0.0040 0.0060 0.0080

-0.0040

-0.0020

0.0000

0.0020

0.0040

t

u (t)4

solução com 2 modos

solução com 3 modos

Figura 21: Curva do movimento vibratório da extremidade livre da barra ao longo do tempo, comresposta usando os três modos na análise modal e apenas os dois primeiros modos.

No Exemplo 6, Figura 18, vimos que a aproximação do terceiro modo com uma malha de trêselementos de barra é bastante pobre. Sua inclusão na análise modal é então desaconselhavel, e pode-se obter melhores resultados usando apenas os dois primeiros modos em (189), como visto na Figura21. A contribuição dos modos mais altos, neste caso, consiste apenas em gerar irregularidades nascurvas. Observe na eq.(189) que a amplitude da contribuição do terceiro modo no deslocamento donó 4, u4(t), é pequena quando comparada às demais, e sua freqüência é alta. Sua eleiminação entãonão afeta sensivelmente as amplitudes de u4(t), mas apenas a suavidade da curva pela eliminação daparcela de alta freqüência.

8.3.2

Exemplo 8 Solução Analítica para Barra sob Deslocamento InicialUse a solução analítica obtida no Exemplo 6 para os modos naturais de uma barra engastada

numa das extremidades, e obtenha a solução analítica da resposta dinâmica da barra devido aodeslocamento inicial prescrito, dado no Exemplo 7, eq.(190). (Observe que este problema usa umasérie de resultados da teoria de equações diferenciais parciais. Se o estudante não se sentir confortavelcom o assunto pode simplesmente pular para a solução, eq.(197).)

Solução:O deslocamento u(x, t) da barra é a solução do problema (159) visto no Exemplo 6. Alí usamos

o método de separação de variáveis, isto é, supusemos que a solução possa ser colocada na formau(x, t) = X(x)T (t), o que resultou em duas equações diferenciais ordinárias, eqs.(161), uma em xoutra no tempo. A primeira equação foi já resolvida no Exemplo 6, onde obtivemos as freqüênciasnatuais ωj e os modos Xj(x) mostrados nas eqs.(164) e (165).Temos agora que resolver o problema no tempo, a segunda das eqs.(161). Para isto temos que

primeiro identiÞcar as condições iniciais na barra. Se aplicarmos uma força F na extremidade,


temos um problema estático, onde o deslocamento u(x) em cada ponto é u(x) = Fx/EA. A veloci-dade inicial é nula. Então o problema no tempo é deÞnido por:

d2T (t)

dt2= ω2T (t) = 0, para t > 0. (191)

A solução deste problema é conhecida e tem a forma

T (t) = C senωt+D cosωt. (192)

Da solução do problema em x temos que inÞnitas freqüências naturais ω satisfazem o problema, comovisto em (164). Então inÞnitas soluções existem para (192), e cada uma delas tem a forma:

Tn = Cn senωnt+Dn cosωnt. (193)

A solução do problema original da barra tem então a forma u(x, t) =P∞

n=1Xn(x)Tn(t), isto é,

u(x, t) =∞Xn=1

sen(2n− 1)πx

2L

ÃCn sen

(2n− 1)π2L

sE

ρt+Dn cos

(2n− 1)π2L

sE

ρt

!. (194)

As constantes Cn e Dn devem ser determinadas de forma a fazer cm que u(x, t) satisfaça ascondiçòes iniciais, que são: ⎧⎪⎪⎨⎪⎪⎩

u(x, 0) =Fx

EA,

du(x, t)

dt

¯t=0

= 0.

(195)

A única forma de (194) satisfazer a segunda condição é que todos os Cns sejam nulos. A primeiracondição resulta no seguinte:

u(x, 0) =∞Xn=1

Dn sen(2n− 1)πx

2L=Fx

EA. (196)

A forma de determinar os Dns consiste em multiplicar esta equação por sen (2n − 1)πx/2L eintegrar no intervalo 0 ≤ x ≤ L. Observe que:Z L

x=0

sen(2m− 1)πx

2Lsen

(2n− 1)πx2L

dx =

½0 se m 6= n,L(−1)m2π

se m = n.

Então para um dado m, = 3, por exemplo, o somatório em (196) reduz-se a um único termonão-nulo, o termo n = m, no caso 3. Então,Z L

x=0

Dn sen2(2n− 1)πx

2Ldx =

Z L

x=0

Fx

EAsen2

(2n− 1)πx2L

dx,

o que resulta em

Dn =8FL(−1)n−1π2EA(2n− 1)2 .

A solução (194) para o movimento da barra devido ao deslocamento inicial então Þca:


u(x, t) =8FL

π2EA

∞Xn=1

(−1)n−1(2n− 1)2 sen

(2n− 1)πx2L

cos(2n− 1)π2L

rE

ρt. (197)

8.4 Análise Modal Geral

Consideramos aqui a situação mais geral que (127), o sistema de equações diferenciais de movimentode um sistema com amortecimento viscoso, carregado, com n graus de liberade e devidas condiçõesiniciais dadas por ¯

¯ Mu(t) +C úu(t) +Ku(t) = F(t),u(0) = uo,

úu(0) = vo,

(198)

onde C é a matriz de amortecimento do sistema, de ordem n× n, simétrica. Na seção a seguirveremos métodos para a determinação da matriz de amortecimento de um sistema. A princípiopodemos tentar aplicar o mesmo processo de análise modal usado no caso de vibrações livres não-amortecidas. Primeiro deÞnimos e resolvemos o problema de autovalor£

K− ω2jM¤φj = 0. (199)

Conhecidos os n autovetores formamos amatriz modal Φ. Usamos essa matriz para transformaros deslocamentos em novas funções η(t), como na eq.(169). Substituímos essa em (198) e pré-multiplicamos o resultado por ΦT . Usamos as relações de ortogonalidade (148) e (151), a eq.(198)Þca na forma:

η(t) +ΦT C Φ úη(t) +Λ2 η(t) = ΦT F(t). (200)

Lembremos que o objetivo da transformação da equação do movimento com o uso da matriz modalé o de obter um conjunto de equações desacoplado que possa ser resolvido uma a uma. Em (200) asmatrizes coeÞcientes de η e η são diagnoais, porém se C for uma matriz qualquer, o resultado deΦT C Φ será uma matriz n × n não-diagonal, o que não facilita em nada a resolução do sistema.Uma forma de contornar este problema consiste em usar uma matriz de amortecimento C construídade forma especial de tal maneira a se saber previamente que ΦT C Φ será diagonal, ou pelo menosquase diagonal (isto é, com uma banda bastante estreita, como as matrizes tridiagonais por exemplo).Consideramos nesse ponto que sabemos como determinar C tal que ela possa ser diagonalizada. Deforma geral, deÞne-se a matriz transformada c e o vetor força transformado f(t) como:

c = ΦT C Φ e f(t) = ΦTF(t). (201)

O sistema de equações de movimento transformado (200) então Þca

η(t) + c úη(t) +Λ2η(t) = f(t). (202)

Se tivermos c diagonal, este é um sistema desacoplado, do tipo do sistema (170), com n equaçõesna forma


ηj(t) + cj úηj(t) +ω2jηj(t) = fj(t), j = 1, 2, ..., n, (203)

onde cj = cjj. Novamente, cada uma destas n equações diferenciais é idêntica à equação do problemade um grau de liberdade. Comparando com (47) reescrevemos cada cj em termos da freqüêncianatural ωj e da taxa de amortecimentos ζj, isto é,

cj = 2ζjωj, (204)

tal que a matriz c tem a forma

c =

⎡⎢⎢⎢⎣2ζ1ω1

2ζ2ω2. . .

2ζnωn

⎤⎥⎥⎥⎦ .A solução do problema desacoplado é composta pela sobreposição de duas parcelas. Uma parte

consiste na solução do problema de vibrações livres sob carregamento inicial, mais a solução doproblema de vibração forçada sob condições iniciais nulas. Em outras palavras a solução do problema(182) é equivalente à adição das soluções dos seguintes problemas lineares:

⎧⎨⎩ Mu(t) +C úu(t) +Ku(t) = 0,u(0) = uo,úu(0) = vo,

e

⎧⎨⎩ Mu(t) +C úu(t) +Ku(t) = F(t),u(0) = 0,úu(0) = 0,

Para um grau de liberdade as soluções destes problemas aparecem na equação (77). Essaequação dá então a solução de cada uma das equações de movimento desacopladas em (203). Asolução é

ηj(t) =1

ωd j

Z t

o

fj (τ) e−ζjωj(t−τ) senωdj(t− τ) dτ

+ηj(0)q1− ζ2j

cos¡ωdjt− φj

¢+úηj(0)

ωd jsenωdjt, j = 1, . . . , n,

(205)

onde ηj(0) e úηj(t) são calculados por (178) usando os valores conhecidos de uo e vo , f(t) =

ΦTF(t), ωdj = ωj

q1− ζ2j de (30), e de (33) se tem tanφj = ζ/

p1− ζ2. ωj e φj são os auto-

pares do problema de autovalor associado (199). Obtidos os valores de η(t) num dado instante, asolução do problema de vibrações forçadas amortecidas é dada por (169):

u(t) = Φη(t). (206)

8.4.1 Exemplo 9 Solução pelo MEF de barra sob Carga Variável no Tempo

Considere a barra do Exemplo 7 sob condições iniciais nulas, isto é, uo = vo = 0, submetida a umaforça dinâmica aplicada na extremidade de 105sen 4.000t, em Newtons. Determine a resposta dosistema.


Solução:Devemos formar o sistema desacoplado de equações do movimento, eqs.(203). Com o auxilio de

(185) no Exemplo 7, basta calcularmos as forças generalizadas usando de (201):

f(t) = ΦTF(t) =

r18

ρAL

⎡⎣ 0, 170518 1/√6 0, 271037

0, 295345 0 −0, 4695360, 341035 −1/√6 0, 542173

⎤⎦T ⎡⎣ 00105

⎤⎦ senΩt, (207)

onde Ω = 4.000 s−1 é a freqüência de excitação. O sistema desacoplado é obtido usando as freqüênciasnaturais já obtidas no Exemplo 3, página 37:

⎡⎣ η2(t)η3(t)η4(t)

⎤⎦+⎡⎣ 6, 47 · 107 6, 92 · 108

2, 28 · 109

⎤⎦⎡⎣ η2(t)η3(t)η4(t)

⎤⎦ =⎡⎣ 163, 811−196, 093260, 464

⎤⎦ senΩt. (208)

As condições iniciais dos deslocamentos generalizados vem de (178), que resultam em: ηo = úηo =0. A solução então para cada equação vem da integral de Duhamel em (205). Por exemplo, paraη2(t), com amortecimento nulo, ξj = 0, temos:

η2(t) =1

ω2

Z t

o

f2 (τ) senω2(t− τ) dτ,

=1, 64 · 10526.306

Z t

o

senΩt sen26.306(t− τ) dτ.

Para os demais modos o processo é o mesmo. Podemos deÞnir o vetor carregamento temporal comof(t) = R senΩt. A solução analítica para um ηj qualquer é:

ηj(t) =Rj

ωj(Ω2 − ω2j)[Ω senωjt− ωj senΩt]. (209)

Observe que apenas em casos bastantes simples a integral de Duhamel poderá ser feita analitica-mente como feito aqui. Em geral ela será estimada por integração numérica.A solução do sistema a cada instante em termos de deslocamentos nodais físicos vem de (206),

u(t) = Φη(t): ⎡⎣ u2(t)u3(t)u4(t)

⎤⎦ = Φ⎡⎣ 33, 63 sen 8.045t− 16, 72 sen 8.045t−2, 90 sen 26.311t+ 0, 441 sen 26.311t1, 151 sen 47.733t− 0, 0965 sen 47.733t

⎤⎦ 10−4.A Figura 22 ilustra o movimento descrito pelo nó 4 ao longo do tempo.

8.5 Resumo do Método de Análise Modal

De um ponto de vista teórico o processo de cálculo descrito acima, o chamado método de análisemodal, é capaz de fornecer a resposta dinâmica de um sistema. Consideremos porém a situação maiscomum nos dias de hoje. Busca-se a modelagem de sistemas complexos de forma tão detalhada, que aúnica ferramenta possível para a modelagem do sistema é o método de elementos Þnitos, tanto pelas


0.000 0.002 0.004 0.006 0.008

t

-20.0

-10.0

0.0

10.0

20.0

Figura 22: Resposta no tempo do movimento do nó 4 na extremidade da barra sob carga harmônica.(Deslocamento dividido por

p18/ρAL.)

características do método quanto, principalmente, pelo nível de detalhamento que ele permite. Ocorreque, via de regra, chega-se a matrizes estruturais (K,M e C) de ordens elevadas, em modelos quepodem variar desde menos de mil até algumas centenas de milhares de graus de liberdade. Quandose observa a transformação (169) nota-se que todo o processo de cálculo descrito acima, é baseadona prévia solução de um problema de auto-valor associado, exigindo a determinação de todos osn auto-pares do problema. Na prática esta determinação completa é basicamente inviável. Numaseção a seguir examinaremos alguns métodos para a solução de problemas de autovalor de grandeporte provenientes de modelagens de elementos Þnitos. Tornar-se-á evidente o alto custo e tempocomputacional envolvidos, mesmo nos melhores métodos.Por outro lado consideremos que os auto-pares em Λ2 e Φ foram organizados em ordem crescente,

isto é, ω21 < ω22 < . . . < ω

2n e os correspondentemente em Φ. O processo de transformação (169) das

equações de movimento pode ser visto como uma expansão das funções força e solução em termosda base de vetores deÞnida por Φ no espaço vetorial n-dimensional deÞnido por K eM, da seguinteforma:

u(t) = Φη(t) = η1(t)φ1 + η2(t)φ

2 + . . .+ ηn(t)φn,

F(t) = Φg(t) = g1(t)φ1 + g2(t)φ

2 + . . .+ gn(t)φn.

Cada termo no somatório representa a contribuição da freqüência correspondente no valor de u(t)e F(t). Para a maioria dos carregamentos as contribuições das várias freqüências geralmente sãomaiores para os baixos modos e tendem a decrescer para altas freqüências. Isto signiÞca que aimportância dos termos η1(t)φ

1e g1(t)φ1 na composição do deslocamento e da força aplicada são

muito mais impostantes que as contribuições do último modo, ηn(t)φn e gn(t)φ

n.Ao mesmo tempo, a modelagem de um sistema complexo por qualquer método, como o de elemen-

tos Þnitos, sempre produzirá auto-pares com precisão decrescente para os modos mais altos. Frquente-mente, mesmo que se pague o preço de uma determinação completa dos n autopares, os resultadosobtidos pelos últimos 2/3 dos modos são geralmente classiÞcados apenas como lixo numérico.Esses dois fatos permitem então que a análise modal possa ser feita usando não todos os n modos,

mas apenas os m primeiros. Para a determinação de apenas os m primeiros autopares, m¿ n, doismétodos são mais utilizados, com suas muitas variações disponíveis: o método da interação sub-espacial e o de Lanczos, que serão descritos em seções posteriores.As etapas de cálculo do método geral de análise modal são descritas abaixo.


ETAPA 1 Determinação das matrizes da equação de movimento, eqs. (198) e vetor força F(t),

Mu(t) +C úu(t) +Ku(t) = F(t),

e identiÞcação das condições iniciais e condições de contorno do problema. Para componentes com-plexos e estruturas podem-se usar elementos Þnitos:

ETAPA 2 Formar o autoproblema (130):£K− ω2jM

¤φj = 0.

Usar um dos métodos disponíveis, da interação sub-espacial ou Lanczos, por exemplo, e deter-minar os m primeiros autopares. Colocá-los em ordem crescente. Em caso de autovalores repetidos,usar o método de Gram-Schmidt para ortogonalizar os autovalores. Normaliza-los pela massa, istoé, gerar as matrizes reduzidas

Λ = diag pω21 ω22 . . . ω2my, m×m,Φ =

£φ1φ2 . . . φm

¤, n×m,

tal que

ΦTMΦ = Im×m,

ΦTKΦ = Λm×m.

ETAPA 3 Obter o problema desacoplado (203) agora em sua forma reduzida, com m equaçõese m incógnitas ηj(t):

ηj(t) + 2ζjωj úηj(t) + ω2jηj(t) = fj(t), j = 1, . . . ,m,

onde, de (201),

fj(t) = φjTF(t),

2ζjωj = φjTCφj.

ETAPA 4 Transformar as condições iniciais usando (178):

ηo = ΦTMuo,

úηo = ΦTMvo.

ETAPA 5 Para cada modo j, onde j = 1, 2, . . . ,m, calcular:

ωdj = ωj

q1− ζ2j de (30),

tanφj =ζj√1−ζ2j

de (33).

ETAPA 6 Em cada instante t calcular, de (205):


ηj(t) =1

ωdj

Z t

o

fj (τ) e−ζjωj(t−τ) senωdj(t− τ) dτ

+ηj(0)q1− ζ2j

cos¡ωdjt− φj

¢+ηj(0)

ωdjsenωdjt, j = 1, . . . ,m.

Em geral, em caso de carregamento genérico, a integral de Duhamel que aparece aqui deve serintegrada numericamente.

ETAPA 7 Em cada instante calcular os n deslocamentos físicos nodais de (169):

u(t) = Φη(t).

9 Determinação do Amortecimento

9.1 Um Grau de Liberdade

Consideremos a equação de movimento para um sistema de um grau de liberdade na forma (15) comsua solução (29). O termo exponencial é responsável pela redução na amplitude das oscilações aolongo do tempo. Quanto maior o quociente de amortecimento ζ, mais rápida a atenuação. Pode-sebuscar uma relação entre o quociente das amplitudes em distintos picos com o valor de ζ.Tomemos a eq.(29). O p-ésimo pico ocorre no instante tp. Os picos podem ser identiÞcados pelos

instantes em que cos θp atinge valor +1, isto é, θp = 2πp. Então,

ωdtp − φ = 2πp −→ tp =2πp+ φ

ωd, (210)

como esquematizado na Figura 3, página 8. (Note que a posição destes picos de fato é alterada pelotermo exponencial em (29).) Tomemos agora a amplitude xp no pico p e a amplitude xp+q, q picos afrente. A relação entre eles é obtida usando (29) e simpliÞcando:

xpxp+q

=e−ζωntp

e−ζωntp+q=

e− ζωn

ωd(2πp+φ)

e− ζωn

ωd(2(p+q)π+φ)

= e2qπζωnωd .

DeÞne-se decaimento logaritmico δ como:

δ = lnxpxp+q

=2qπζωnωd

=2qπζp1− ζ2

. (211)

Então δ é a porcentagemde decaimento nas amplitudes após q ciclos. Para ζ = 8% por exemplo, odecaimento é de δ = 50% em apenas 1 ciclo. Para um decaimento de 50% pode-se obter uma relaçãoentre o amortecimento ζ e o número de ciclos q necessários

q =

p1− ζ24πζ

. (212)

Esta relação é plotada na Þgura 23.O signiÞcado do termo quociente de amortecimento para ζ pode ser entendido da seguinte

forma. Para baixos valores de amortecimento, o último termo de (211) pode ser aproximado por


cnxpxp+m

≈ 2qπζ −→ xpxp+m

≈ e2qπζ .

0.00 0.05 0.10 0.15 0.20ζ

0.00

2.00

4.00

6.00

q

Figura 23: Número de ciclos necessários para um decaimento δ de 50% para dado valor de amortec-imento ζ.

O exponencial pode ser expandido em série de Taylor, tomando apenas os dois primeiros termos,os lineares. Então, entre dois ciclos subsequentes,

xqxq+1

≈ 1 + 2πζ −→ ζ ≈ xq − xq+12πxq+1

. (213)

Então ζ tem aproximadamente o signiÞcado da variação da amplitude sofrida entre dois ciclos sub-sequentes, dividido pela amplitude Þnal.

A equação (211) é a base de um dos métodos experimentais mais simples e mais usados na deter-minação do quociente de amortecimento ζ. É o chamado método do decaimento em vibraçõeslivres, e consiste no seguinte. Coloca-se a estrutura sob vibrações livres usando um procedimentoadequado e mede-se a amplitude do movimento em um intervalo de m ciclos. De (211) calcula-seentão o amortecimento:

ζ =δq

δ2 + (2πm)2, δ = ln

xpxp+m

. (214)

Uma série de outros métodos experimentais são disponíveis, como o da ampliÞcação resso-nante, o método da largura de banda, da perda de energia por ciclo e do amortecimentohisteretico, e podem ser vistos em textos padrão de dinâmica.

9.2 n-Graus de Liberdade Elementos Finitos

A matriz de amortecimento C em geral não pode ser modelada por elementos Þnitos com a mesmafacilidade com que se obtém as matrizes massa e rigidez. Formalmente podemos seguir o mesmo


procedimento usado na determinação dos demais termos da equação do movimento. Consideremosque o amortecimento é em última instância um fenômeno relacionado às intereções entre os compo-nentes da estrutura molecular do material. Desta forma é um parâmetro distribuído por unidade demassa ou volume, como a densidade. Numa região de material de volume inÞnitesimal dV , a forçade amortecimento é dada por

dfa = c (x, y, z) úu (xyz) dV. (215)

Este elemento de força pode ser incorporado diretamente à expressão do princípio do trabalho virtual,junto às forças de corpo e de inércia, levando às equações discretizadas de movimento do MEF. Porexemplo para o elemento de barra da seção 7, a equação do P.T.V. (94) se tornaria

AE

Z L

o

∂u

∂x(x, t)

du(x)

dx

Z L

o

(b (x, t)− ρu− c (x) úu) u(x) dx−Af(t)u(L) = 0 ∀u ∈ V ar. (216)

O tratamento de todos os termos é o mesmo visto naquela seção restando o termo do amorteci-mento. Expandindo úu(x, t) e u(x) por (96) e (97) em termos de valores nodais e funções de interpo-lação, o termo de amortecimento Þca:

δIa = A

Z L

o

c (x) úu (x, t) u(x) dx =Xe

A

Ze

c (x) úu (x, t) u(x) dx.

Em cada elemento e esta parcela Þca

δIea = Ace

Ze

úu (x, t) u(x) dx = Ace

Ze

Ã2Xi=1

úui(t)ϕei (x)

!Ã2Xi=1

uiϕei (x)

!dx,

onde o amortecimento c(x) foi tomado como constante em cada elemento e. Isto resulta na matrizde amortecimento do elemento e deÞnida por:

Ceij = Ace

Ze

ϕeiϕej dx. (217)

Essa expressão então se junta às demais associadas K eM na eq.(114).

9.3 Métodos Experimentais

A prática de modelagens numéricas, entretanto, raramente utiliza o procedimento acima. A determi-nação do amortecimento em cada região da estrutura, para cada grau de liberdade não é algo simplesde ser feito como o é a determinação de densidade e propriedades elásticas do material. Além disto,de fato, em estruturas metálicas, o amortecimento de material é em geral desprezivel, ea fonte importante de dissipação de energia se dá nas interfaçes das junções entre diferentes partes.Desta forma, em lugar do procedimento visto na seção acima, o quociente de amortecimentoζj para alguns modos j são escolhidos ou determinados experimentalmente. Numa etapa deprojeto, onde evidentemente não há ainda uma estrutura construída, os amortecimentos podem serescolhidos. Para a maioria das estruturas metálicas ou de concreto ζ Þca na faixa 0,01 a 0,20.Resultados experimentais de modelos semelhantes já construidos podem ser usados como referência.A determinação experimental para uma estrutura real que apresenta inÞnitos modos de vibração,

evidentemente não é simples, mesmo que se busque apenas os amortecimentos para os poucos primiros


modos. Os métodos experimentais para um grau de liberdade comentados no início da seção podemser usados para a estrutura tridimensional para a determinação de um valor representativo de ζ,como ilustrado no exemplo seguinte.

9.3.1 Exemplo 10 Determinação Experimental de Amortecimento

Considere um edifício como na Figura 24, idealizado como um pêndulo amortecido composto por umamassa m, uma mola e um amortecedor, ambos de ßexão. Um teste de vibrações livres é realizado,onde o topo do edifício é deslocado lateralmente, por um macaco hidráulico, por exemplo, que ésubitamente removido. O deslocamento inicial foi medido como sendo u = 5, 0mm para uma forçaaplicada de F = 100 kN. Após a liberação a máxima amplitude no ciclo seguinte foi de u2 = 4, 0mm eo período foi T = 1, 5 s. Obtenha uma estimativa para o quociente de amortecimento ζ da estrutura.

u

F

m

c k

u

(a) (b) (c)

Figura 24: (a) Representação física de um edifício, (b) modelo experimental de ßexão de um grau deliberdade, (c) força e deslocamento inicial aplicados.

Solução:Trabalhemos na identiÞcação dos parâmetros do modelo na Figura 24b. Como o modelo é de um

grau de liberdade, a constante da mola ßexural se relaciona à força e deßexão por

k =F

u=100.000 N5 · 10−3 m = 20 · 106N

m.

Note que aqui não é de interesse a massa real do ediÞcio, mas um valor diferente, associado à ßexãoda estrutura modelada por um grau de liberdade. Este valor pode ser estimado a partir dos valoresobtidos experimentalmete, o período e a rigidez:

T =2π

ω= 2π

rm

k−→ m = k

µT

2π

¶2= 20 · 106

µ1.5

2π

¶2= 1, 14 · 106 kg,

e a freqüência de vibração é

ω =2π

T=2π

1, 5= 4, 19 rad/s.

O decaimento logaritmico vem de (211):


δ = lnu

u2= ln

5, 0

4, 0= 0, 223,

e o quociente de amortecimento vem de (214)

ζ =δq

δ2 + (2π)2=

0, 223q0, 2232 + (2π)2

= 0, 0355.

Quando se usa o método de análise modal, a determinação prévia da matriz C de amortecimentonão é essencial, uma vez que os valores de ζj estimados para cada modo são aplicados diretamentenas equações de movimento desacopladas (203). Por outro lado, quando se deve usar um métodocomo o de integração direta por exemplo, necessário num problema não linear, devemos ter algumaforma de estimar a matriz C a partir dos quocientes de amortecimento usados ζj. A determinaçãode C é o assunto seguinte, onde apresentaremos dois métodos clássicos.

9.4 Método Analítico 1 Rayleigh

Para entender o método de Rayleigh para a determinação de C, suponha que se tenha conhecidos osvalores de ζj para alguns dos poucos primeiros p modos do sistema, digamos 1, 2 ou 3. A questãoé: como obter a matriz C correspondente, isto é, aquela que, diagonalizada por (201) resulte noamortecimento previsto ζ1, ζ2, . . . , ζp nestes p modos e nos demais modos apresente amortecimentoscompatíveis (embora não iguais) a estes?Rayleigh apresentou uma solução a este problema. Compare a aproximação de elementos Þnitos

em (217) para a matriz de amortecimento Ce com a deÞnição da matriz massaMe em (104). Nota-seclaramente que ambas as matrizes são proporcionais, na forma Ce = αMe. De forma geral isto não éexato, devido à complexidade dos diversos processos de amortecimento. Numa primeira aproximaçãoentretanto, pode-se considerar que C teria também uma contribuição da rigidez K, isto é

C = aoM+ a1K, (218)

que é a expansão de Rayleigh. Estamos buscando uma estimativa para C, de preferência umaque seja diagonalizavel. De fato esta expansão é diagonalizavel. Para veriÞca-lo basta levar estaexpansão a (201), e usar as relações de ortonormalidade da matriz modal Φ, o que resulta

c = ΦTCΦ = aoI+ a1Λ2, (219)

isto é, c é diagonal. Como existem duas constantes a determinar, ao e a1, é possível impor quocientespara dois modos, ζ1 e ζ2.

Existe uma forma mais geral que (218), diagonalizável, que permite a incorporação de tantosvalores de ζ quantos se queira ou possua. Para p valores disponíveis de ζ deve-se determinar as pconstantes ab em:

C =MP

b ab [M−1K]b =

PbCb, para b = . . .−2,−1, 0, 1, 2,| z

p valores

. . . (220)

Esta forma contém (218) como caso particular onde b = 0 e 1. Os valores de b podem assumirquaisquer valores entre −∞ e +∞, mas os melhores resultados são obtidos tomando os valores emtorno de 0.


Para cada modo j o amortecimento generalizado é dado por (201) e (204),

cj = φjTCφj = 2ζωj, para j = 1, 2, ...,m, (221)

onde m é o quantidade de modos φj usados na análise. Mas como C vem de (220) em termos dosab, a parcela b de cj é

cjb = abφjTM

£M−1K

¤bφj, para j = 1, 2, ...,m. (222)

Observe que, se não se conseguir uma maneira de evitar a inversão de M a formulação acima éinútil, devido ao esforço computacional inadmissível envolvido na inversão completa de uma matrizde alta ordem. Passemos pois a um artifício para contornar o problema. Considera-se o problemado autovalor KΦ =MΦΛ2. Se o pré-multiplicamos porM−1 temos£

M−1K¤1Φ = ΦΛ2. (223)

Pré-multiplicamos agora por [M−1K] e usando (223) obtemos£M−1K

¤2Φ = ΦΛ4,

e pré-multiplicando novamente por [M−1K] e usando (223) obtemos£M−1K

¤3Φ = ΦΛ6.

Observando o padrão destas três equações temos por indução que

[M−1K]bΦ = ΦΛ2b . (224)

Agora pré-multiplicamos por ΦTM e obtemos

ΦTM£M−1K

¤bΦ = Λ2b. (225)

Mas o lado esquerdo é justamente aquele da expansão de C em (220). Observe que agora não maisé necessário inverterM, basta usar as freqüências naturais em Λ2. Tomamos apenas o modo j comoem (222):

cjb = abφjTM [M−1K]bΦj,

cjb = abω2bj , para j = 1, 2, ...,m.

Então, de (220) e (221),

cj =Xb

abω2bj = 2ζjωj, para j = 1, 2, ...,m.

Então,

ζj =1

2ωj

Xb

abω2bj . (226)

Esta expressão deÞne p equações algébricas em termos dos p valores incógnitos ab, conhecidos pvalores de amortecimento ζj. Consideremos alguns casos particulares.


Caso 1 - Suponha que tenhamos conhecido apenas um valor de amortecimento, o do modo 1, ζ1.Neste caso buscaremos o valor de ao. De (226) temos uma única equação:

ζ1 =1

2ω1ao −→ ao = 2ζ1ω1.

Levando b = 0 em (220) temos a expansão de C :

C = aoM. (227)

Caso 2 - Se tivermos valores de ζ para os dois primeiros modos ζ1 e ζ2. Buscaremos ao e a1 de(226) escrevendo duas equações:

1

2

∙ 1ω1

ω11ω2

ω2

¸ ∙aoa1

¸=

∙ζ1ζ2

¸−→

∙aoa1

¸=2

4

⎡⎣ ω2 −ω1− 1

ω2

1

ω1

⎤⎦∙ ζ1ζ2

¸, (228)

onde 4 é o determinante, 4 = ω2/ω1 − ω1/ω2. A solução deste problema dará aoe a1, e a expansãode C em (220) Þca

C = aoM+ a1K. (229)

Caso 3 - Caso tenhamos ζ1, ζ2 e ζ3 os valores de a−1, aoe a1 vem de (226) como a solução de

1

2

⎡⎢⎢⎢⎢⎢⎣1

ω21

1

ω1ω1

1

ω22

1

ω2ω2

1

ω23

1

ω3ω3

⎤⎥⎥⎥⎥⎥⎦⎧⎨⎩ a−1aoa1

⎫⎬⎭ =

⎧⎨⎩ ζ1ζ2ζ3

⎫⎬⎭ . (230)

A expansão de C Þca

C = aoM+ a1K+ a2M£M−1K

¤2. (231)

Devemos buscar mais uma identidade vetorial antes de prosseguir. Tomamos a relação de ortonor-malidade ΦTMΦ = I. Então, ΦTMΦ = Φ−1Φ. Então temos a inversa da matriz modal

Φ−1 = ΦTM. (232)

Pós-multiplicando esta relação em (224) temos£M−1K

¤b= ΦΛ2bΦTM.

Podemos agora retornar a (231) que Þca

C = aoM+ a1K+ a2MΦΛ4ΦTM. (233)

Este procedimento pode ser usado mesmo que a matriz modal não tenha sido completamentedeterminada, isto é, se tiver sido determinado apenas os m primeiros autopares, de forma que em vezde Λ2 e Φ temos em (233) Λ

2e Φ, matrizes m×m e n×m. Porém a eq. (232) perde o signiÞcado.


VeriÞca-se que no caso 1, onde C = aoM, o amortecimento será maior nos primeiros modos eserá mínimo nos mais altos. Caso, C = a1K ocorre o inverso, o amortecimento será maior nos modosmais altos.Observe que em (220), a matriz C será não-simétrica se houver valores de b maiores que 1, uma

vez que M−1K é uma matriz não-simétrica, mesmo que M e K o sejam. Isto explica porque emgeral são usados apenas dois termos, como em (228).

9.5 Método Analítico 2

Um segundo método, diferente daquele baseado no método de Rayleigh é disponível. Consideremosas eqs. (201) e (204)

c = ΦTCΦ = 2

⎡⎢⎢⎢⎣ζ1ω1

ζ2ω2. . .

ζnωn

⎤⎥⎥⎥⎦ . (234)

Se os ζ 0s forem disponíveis, C pode ser obtido invertendo a matriz modal:

C = Φ−TcΦ−1. (235)

A inversa de Φ é tomada de (232), o que dá:

C =MΦcΦTM. (236)

Note que esta matriz contém amortecimento apenas nos modos especiÞcados. Se apenas ζ1, ζ2 eζ3 forem especiÞcados em (234), C não apresentará nenhum amortecimento nos demaismodos, diferentemente do método anterior, eq. (220).

9.5.1 Exemplo 11 Determinação Experimental da Matriz de Amortecimento

Considere a barra do Exemplo 7, página 47, modelada por três elementos, com as matrizes de rigidezK, massaM e modal Φ dadas em (185). a) Determine a matriz de amortecimento C usando ométodo de Rayleigh, usando (218) e as constantes ao e a1 para ζ1 = ζ2 = 0, 010 = ζ; b) Determine amatriz de amortecimento diagonalizada (219).

Solução:As constantes ao e a1 vem de (228). As duas primeiras freqüências naturais foram obtidas no

Exemplo 2: ω1 = 8045 s−1 e ω2 = 2634, 5 s−1. Então,

1

2

∙ 1ω1

ω11ω2

ω2

¸ ∙aoa1

¸=

∙ζ1ζ2

¸−→

∙aoa1

¸=

∙123, 22

5, 8214 · 10−3¸.

A matriz de amortecimento vem de (218), usando K eM de (185), na página 47. Então, C =aoM+ a1K:

C =

⎡⎣ 91, 215 −29, 589 0−29, 589 98, 216 −29, 589

0 −29, 589 45, 608

⎤⎦ . (237)


A matriz diagonalizada é dada por (219), c = ΦTCΦ, isto é, c = aoI+ a1Λ2:

c =

⎡⎣ 160, 9 526, 221449, 6

⎤⎦ . (238)

Entretanto, de (204), cada termo de amortecimento cj se relaciona com o correspondente fator ξje freqüência ωj na forma cj = 2ζjωj. Pode-se comparar que a matriz de amortecimento C gerada em(237) corresponde a ξ1 = ξ2 = 0, 010, conforme tinha sido inicialmente imposto, e também obtemosξ3 = 0, 0152.

9.5.2 Exemplo 12 Vibração Amortecida de Barra sob Deslocamento Inicial

Refaça o Exemplo 7, de vibrações livres de uma barra sob deslocamento inicial prescrito usando ométodo de elementos Þnitos, incluindo agora os efeitos de amortecimento. Faça com que ζ1 = ζ2 =0, 010 e use a matriz C obtida no Exemplo 11.

Solução:Por comodidade repetimos aqui os dados já obtidos para este problema, que tem sido desenvolvido

em diversos exemplos desde o Exemplo 2 (ver eqs.(185), (186) e (238) nas páginas (47) e (64)):

K =ρAL

12

⎡⎣ 2 1 01 4 10 1 2

⎤⎦ ; M =2EA

L

⎡⎣ 1 −1 0−1 2 −10 −1 1

⎤⎦ ;Φ =

r18

ρAL

⎡⎣ 0, 170518 1/√6 0, 271037

0, 295345 0 −0, 4695360, 341035 −1/√6 0, 542173

⎤⎦ ; (239)

Λ2 =

⎡⎣ 8.045, 02 26.311, 52

47.733, 42

⎤⎦ s−2, c =

⎡⎣ 160, 9 526, 221449, 6

⎤⎦ .onde os valores de E, A, ρ e L são dados no Exemplo 2. Os valores dos coeÞcientes de amortecimento,associados a c foram obtidos no Exemplo 11 como: ζ1 = ζ2 = 0, 010 e ζ3 = 0, 0152. Os deslocamentose velocidades generalizadas foram obtidas no Exemplo 7 como ηo = [25, 31;−2, 812; 1, 146]T e úηo = 0.Naquele problema os ângulos de fase φj eram todos nulos uma vez que o amortecimento era nulo.Agora, usando (30) e (30),

ωdj =q1− ζ2jωj −−→ . ωd = [8.045; 26.310; 47.727]

T ,

tanφj =ζjq1− ζ2j

−−→ . φd = [0, 010; 0, 010; 0, 0152]T ,

(240)

Observe que os ângulos de fase são de ordem de meio grau, uma vez que o amortecimento é baixo.Isto faz também com que as freqüências amortecidas ωdj pouco se diferenciem de ωj. A solução doproblema generalizado é dado apenas pelo terceiro termo em (205), que Þca:


ηj(t) = Aj cos¡ωdjt− φj

¢, com Aj =

ηj(0)q1− ζ2j

j = 1, . . . , n.

As amplitudes são:

A = [25, 31; −2, 832; 1, 146]T 10−4. (241)

Uma vez que os deslocamentos nodais físicos são dados por u(t) = Φη(t), a solução obtida porelementos Þnitos é a seguinte:

104

⎡⎣ u2(t)u3(t)u4(t)

⎤⎦ =

⎡⎣ 20, 7335, 9041, 46

⎤⎦ e−80,45t cos(8.045t− 0, 010) +⎡⎣ −101

⎤⎦ 5, 554e−263t cos(26.310t− 0, 01)+

⎡⎣ 1, 492−2, 5842, 984

⎤⎦ e−477t cos(47.727t− 0, 0152) (242)

A Figura 25 mostra o histórico do deslocamento do nó na extremidade da barra, a aproximaçãoda função u4(t) obtida por elementos Þnitos a partir da eq.(242). É visível a atenuação da amplitudeao longo dos ciclos.

0.002 0.004 0.006 0.008 0.01 0.012t

-0.004

-0.002

0.002

0.004

Figura 25: Vibração livre amortecida da extremidade da barra do Exemplo 12: aproximação porelementos Þnitos do deslocamento no nó 4 ao longo do tempo, u4 × t.

9.5.3 Exemplo 13 Vibração Forçada Amortecida pelo MEF

Considere a barra do Exemplo 12, agora sob condições iniciais de deslocamentos e velocidade nulas,uo= vo= 0, submetida a uma força aplicada na extremidade variando no tempo conforme:105 senΩt,com Ω = 4.000 s−1. Determine a aproximação de elementos Þnitos para a resposta do sistema.(Observe que este é o mesmo Exemplo 9, porém incluindo agora o amortecimento). Os dados dosistema estão já sumarizados nos Exemplos 9 e 12. A modelagem de elementos Þnitos é feita comtrês elementos de barra como na Figura 14 na página 14.

Solução:


Os dados do sistema estão sumarizados nas eqs.(239), enquanto que as forças generalizadas emcada nó j, fj(t) do sistema desacoplado foram calculados na eq.(207) e valem:

f(t) = ΦTF(t) =

r18

ρAL

⎡⎣ 0, 170518 1/√6 0, 271037

0, 295345 0 −0, 4695360, 341035 −1/√6 0, 542173

⎤⎦T ⎡⎣ 00105

⎤⎦ senΩt, (243)

=

⎡⎣ 163.811−196.093260.469

⎤⎦ senΩt.Temos assim todos os termos da equação completa do movimento desacoplada, eq.(203). A

solução é dada por (205), onde apenas o termo da integral de Duhamel é não nulo, uma vez que ascondições iniciais são nulas. Cada um dos termos tem a seguinte forma:

ηj(t) =1

ωd j

Z t

o

ajsenΩτ · e−ζjωj(t−τ) · senωdj(t− τ) dτ, (244)

onde aj é a amplitude da força generalizada em cada nó dada em (243). A solução para o históricodos deslocamentos físicos nodais é dada por ui(t) =

Pnj=1Φijηj(t). O deslocamento no nó 4, por

exemplo, na extremidade da barra, é dado por:

u4(t) = −7, 527 · 10−5co + 7, 274 · 10−5e−26,62tc1 + 1, 7697 · 10−6e−263tc2 (245)

+7, 6920 · 10−7e−726tc3 + 6, 375 · 10−3so − 2, 737 · 10−3e−80,5ts1−8, 6424 · 10−5e−263ts2 − 2, 5116 · 10−5e−726ts3

co = cosΩt, c1 = cos 8.044t, c2 = cos 26.309, 7t, c3 = cos 47.727, 5t, so = senΩt, so = senΩt,c1 = sen 8.044t, c2 = sen 26.309, 7t, c3 = sen 47.727, 5t.A Figura 26 mostra a plotagem de u4(t) conforme (245), no intervalo de tempo de 0 a 0,06 s.

Pode-se visualizar um leve decaimento da amplitude devida ao amortecimento. Da equação nota-seque, no limite para t→∞, todos os termos que contem expenencial tendem a zero, deixando apenasos termos em co e so, de forma que u4(t) tende a

u4(t)|∞ = −7, 527 · 10−5 cosΩt+ 6, 375 · 10−3 senΩt,isto é, um movimento com amplitude constante, harmônico, com a mesma freqüência do carrega-mento, Ω = 4.000 s−1.

10 Lista de ExercíciosSeções 1 a 4.

1.1 Deduza a eq.(4) a partir de (5) e (6).

1.2 Deduza as eqs.(8) e (9). (Dica: para obter (9), use (6)1 com sen2 φ+ cos2 φ = 1.)

1.3 Deduza a eq.(30).


0.01 0.02 0.03 0.04 0.05 0.06t

-0.0075

-0.005

-0.0025

0.0025

0.005

0.0075

Figura 26: Solução de elementos Þnitos para o movimento do nó na extremidade da barra no Exemplo13, sob vibração forçada amortecida.

1.4 Deduza as eqs.(38)1 e (38)2 usando (24) e (24), respectivamente, além de (28).

1.5 Deduza a solução geral de um sistema não-amortecido para um carregamento F (t) = kA cos(ωnt+ψ), sob condições iniciais nulas (xo = vo = 0). (a) usando (54); (b) usando (51). (Solução:x(t) = −A cosψ cosωnt+

βA(1−β2)sen ψ sen ωnt+

A(1−β2) cos(ωt+ ψ).)

1.6 Mostre a relação entre as partes reais de (51) e (54).

Seção 6

1.7 Pesquise em livros de equações diferenciais, a demonstração da propriedade (74). (Dica: façasubstituição de variáveis, por exemplo, ρ = t− τ).)

1.8 Considere um sistema de um grau de liberdade, não amortecido, sujeito a um carregamento dotipo:

F (t) = Fo senωt para t > 0,

F (t) = 0 para t < 0.

Calcule a resposta do sistema usando integral de Duhamel. x(0) = úx(0) = 0. Compare com asolução (53).

1.9 Considere um sistema de um grau de liberdade, amortecido, sujeito a um carregamento do tipo:

F (t) = u (t− a) =(0 se t < a,

f se t > a.

(u (t− a) é a chamada função degrau unitário, ou função Heaviside.) Calcule a respostado sistema para a = 0 e x(0) = úx(0) = 0 (Resposta: x(t) = 1

k

h1− e−ζωnt

³cosωdt+

ζωnωdsenωdt

´iF (t)

1.10 Use a solução do problema 1.9 e obtenha a solução para o carregamento da Figura 27.(Resposta:x(t) = t

k[1− cosωn(t+ T )]u(t+ T )− [1− cosωn(t− T )], onde u(α) é a função degrau do

Exercício 1.9.)


F(t)

f

0 T - T t

Figura 27: Carregamento temporal do Exercício 1.10.

1.11 Obtenha a eq.(80).

1.12 Integre as deÞnições das matrizes massa, rigidez e força do elemento de barra e obtenha (106).

1.13 Integre as equações (112) usando as funções de interpolação lineares (95), e determine as ex-pressões das energias no elemento de barra, eqs.(113).

1.14 Use a integral de convolução para obter a solução transiente de um sistema de um grau deliberdade amortecido sob carga f(t) = At, com condições iniciais uo = vo = 0.

1.15 Idem para uma força na forma da Figura 28.

F(t)

0

F

T t

Figura 28: Carregamento para o Exercício 15.

1.16 Considere o caso não amortecido onde o carregamento é impulsivo na forma da Figura 29, dadopor F (t) = kA senωt para t < t1 e F (t) = 0 para t > t1.

a) Determine a resposta do sistema para t < t1.

b) Use os valores de u(t1) e úu(t1) como valores iniciais no intervalo t > t1 e obtenha a respostado sistema em t > t1.

c) Diferencie a resposta a) e determine o instante e o valor da deßexão máxima devida aoimpulso aplicado.(Resposta: a) ua(t) = A[senωt − β senωnt]/(1 − β2); b) u(t) = ua(t1)

ωnsenωn(t − t1) +

ua(t1) cosωn(t− t1); c) ωt = 2π/(1 + β).)


F(t)

t 1

F(t) = KA sen ωt

Figura 29: Carregamento do Exercício 1.16.

1.17 Considere o caso não amortecido onde o carregamento impulsivo é na forma indicada na Figura30. Resolva os ítens a) e b) do problema anterior.

(Resposta: a) ua(t) = A(1−cosωnt); b) ub(t) = úua(t1)ωn

(1− senωn(t− t1))+ua(t1) cosωn(t− t1).)

F(t)

t

kA

t1

Figura 30: Carregamento do Exercício 17.

Seção 8

1.18 Observe as matrizes da barra do Exemplo 2 para dois e três elementos Þnitos. Determine opadrão ou regra de deÞnição dos termos das matrizes e mostre como seriam para um númeroarbitrário de elementos.

1.19 Use os resultados do Exercício anterior e determine as freqüências naturais para a barra embalanço no Exemplo 3, para 4, 5,..., 10 elementos. Complete a Tabela 1 e plote os erros relativospara cada modo versus o número de nós do modelo.

(Solução analítica: ωn = (2n− 1)πpE/ρ/2L para o modo n.)

1.20 Considere o sistema da Figura 31.

a) Mostre que as matrizes do sistema são:

K = k

∙2 −1

−1 2

¸; M = m

∙1 00 2

¸.


2 m = m

1 m = 2m

2

(t) 2 u (t) u 1

k = k k k = 1

k = k 3

Figura 31: Sistema do Exercício 1.20.

b) Calcule as freqüências naturais do sistema;

c) Calcule os modos de vibração.

(Resposta: w1,2 = [3(1∓1/√3)/2]1/2

pk/m; φ1 = [1, 0; 1, 366025]T , φ2 = [1, 0; −0, 366025]T .)

1.21 Calcule os autopares deÞnidos pelas seguintes matrizes:

K =

∙3 −1

−1 2

¸; M = m

∙1 00 0, 5

¸.

(Resposta: w21 = 2; w22 = 5; φ

1 = [1, 0; 1, 0]Tp2/3; φ2 = [1, 0; −2, 0]Tp1/3.)

1.22 Considere o sistema do problema anterior com as seguintes condições iniciais: uo = [1/3; 1]T ,vo = [0; 0]

t. Obtenha a resposta do sistema sob vibração livre não amortecida.

1.23 Considere um sistema não-amortecido com as matrizes

K =

⎡⎣ 1 −1 0−1 3 −20 −2 5

⎤⎦ ; M = m

⎡⎣ 1 0 00 1, 5 00 0 2

⎤⎦ .em vibração livre sob as seguintes condições iniciais: uo = [0, 5; 0, 4; 0, 3]T , vo = [0; 9; 0]T .

a) Determine os autopares do problema;

b) Determine as equações desacopladas de movimento e sua solução.(Resposta: Λ2 = Diag[14, 52; 31, 12; 46, 12]; φ1 = [1, 0; 0, 644; 0, 3]T/

√1, 181; φ2 =

[1, 0; −0, 601; −0, 676]T/√2, 455; φ3 = [1, 0; −2, 57; 2, 47]T/√23, 1; η1 = 0, 332 senω1t+0, 592 cosω1t; η2 = −0, 106 senω2t −0, 108 cosω2t; η3 = −0, 033 senω3t +0, 019 cosω3t.)

1.24 Considere o sistema do Exercício 21 com condições iniciais nulas, mas com um carregamentodado por F(t) = [0; 5]T . (Resposta: η1 = 2, 5

p2/3(1 − cosω1t); η2 = −2/√3(1 − cosω2t);

u(t) = Φη.)

1.25 Considere a seguinte equação do movimento:

m

⎡⎣ 1 0 00 1 00 0 1

⎤⎦⎡⎣ u1u2u3

⎤⎦+ k⎡⎣ 2 −1 0−1 2 −10 −1 1

⎤⎦⎡⎣ úu1úu2úu3

⎤⎦ =⎡⎣ 000

⎤⎦ .


Determine oa autopares do autoproblema associado. (Resposta: Λ2 = Diag[0,1982; 1,552,3,252]k/m; φ1 = [0, 328; 0, 591; 0, 737]T/

√m; φ2 = [0, 737; 0, 328; −0, 591]T/√m; φ3 =

[0, 591; −0, 737; 0, 328]T/√m.)1.26 Plote e compare a solução analítica da barra em vibração livre sob deslocamento inicial pre-

scrito, eq.(197) obtida no Exemplo 8, com a solução de elementos Þnitos, eq.(189), obtida noExemplo 7 com três elementos.

Seção 9

1.27 Obtenha as constantes ao e a1 para o amortecimento de Rayleigh. As freqüências naturais são[14, 5; 31, 1; 46, 1]s−1 (do Exercício 23). Imponha quocientes de amortecimento ao primeiro eterceiro modos de ζ1 = ζ3 = 0, 05. Determine o valor previsto para ζ2. (Resposta: ao = 1, 10,a1 = 0, 00165, ζ2 = 0, 0433.)

1.28 Resolva o problema 23 com amortecimentos ζ1 = ζ3 = 0, 05, submetido a um carregamentoF(t) = [0; 0;A senωt]T . Use os resultados do Exercício anterior.

1.29 Resolva o problema 25 com amortecimento ζ1 = ζ2 = 0, 05, sob vibração livre. Considere queos deslocamentos iniciais são aqueles devidos a uma carga estática F(t) = [0; 0;Ak]T . Comparecom a solução não amortecida.

1.30 Considere a treliça mostrada na Figura 32. A estrutura foi modelada numericamente e foideterminado que seu segundo modo de vibrações tem a forma esboçada na Þgura (b). Aplicou-se um carregamento estático nos pontos a e b da estrutura como ilustrado, de forma a simular osegundo modo de vibrações. A carga foi subitamente eliminada e a estrutura pos-se a mover emvibrações livres. Os deslocamentos foram monitorados e as amplitudes a cada ciclo registradas.VeriÞcou-se que, após dois ciclos a amplitude passou de 1,9 para 1,0. Determine uma estimativapara o quociente de amortecimento do segundo modo, ζ2. (Resposta: ζ2 = 0, 051.)

a

bFb

Fa

(t)bx

Figura 32: Treliça do Exercício 30.

REFERÊNCIAS 72

Referências[1] BATHE, K.-J. Finite element procedures in engineering analysis, Prentice-Hall, Inc., N.J., USA,

1982.

[2] BILLINGTON, E.W. The physics of deformation and ßow, McGraw-hill Publ. Co., 1981.

[3] BRUSH, D.O., ALMROTH, B.O. Buckling of bars, plates and shells, McGraw-Hill.

[4] CLOUGH, R.W., PENZIEN, J. Dynamics of structures, McGraw-Hill Book Co., Auckland,1982.

[5] COOK, R.D., MALKUS, D.S., PLESHA, M.E. Concepts and applications of Þnite elementanalysis, John Wiley & Sons, N.Y., 1988.

[6] CULLUM, J.K., WILLOUGHBY, R.A. Lanczos algorithms for large symmetric eigenvalue com-putations, Vol I Theory, Vol II Programs, Birkhauser, Boston, 1985.

[7] DAVIS, P.J. RABINOWITZ, P. Methods of numerical integration, 2nd. ed., Academic Press,Inc., Orlando, USA, 1984.

[8] DYM, C.L., SHAMES, I.H. Solid mechanics - a variational approach, McGraw-Hill Co., 1973.

[9] GOLUB, G.H., van LOAN, C.F.matrix computations, 2nd. ed., Johns Hopkins University Press,USA, 1989.

[10] HACKBUSH, W. Iterative solution of large sparse systems of equations, Springer-Verlag, N.Y.,1994.

[11] HUGHES, T.H. The Þnite element method - linear static and dynamic Þnite element analysis,Prentice-Hall, Inc., N.J., 1987.

[12] JOHNSON, C. Numerical solution of partial differential equations by the Þnite element method,Cambridge University Press, 1992.

[13] KREIDER, D., KULLER, R.G., OSTBERG, D.R., PERKINS, F.W. An introdution to lienaranalysis, Addison-Wesley publ. Inc., 1966.

[14] KREYSZIG, E. Advanced engineering mathematics, 1st. ed., John Wiley & Sons, Inc., N.Y.,1967.

[15] LANGHAAR, H.L. Energy methods in applied mechanics, John Wiley and Sons, Inc., NY, 1962.

[16] MALVERN, E.L. Introduction to the mechanics of a continuous medium, prentice-hall, Inc.,N.J., 1969.

[17] MARGUERRE, K., WOERNLE, H.-T. Elastic plates, Blaisdell Publ. Co., Mass., 1969.

[18] MEIROVITCH, L. Elements of vibration analysis, McGraw-Hill Kogakusha, Ltd., 1975.

[19] ORTEGA J.M Matrix theory - a second course, Plenum Press, N.Y., 1987.

[20] ORTEGA, J.M., RHEINBOLDT, W.C. Iterative solution of nonlinear equations in severalvariables, Academic Press, N.Y., 1970.

REFERÊNCIAS 73

[21] PRESS, W.H., FLANNERY, B.P. , TEVKOLSKY, S.A. Numerical recipes - the art of scientiÞccomputing, Cambridge University Press, Cambridge, 1986.

[22] SZABÓ, B., BABUSKA, I. Finite element analysis, John Wiley & Sons, Inc., N.Y., 1991.

[23] TIMOSHENKO, S. P., KRIEGER, W.S. Theory of Plates and Shells, McGraw-Hill Book Com-pany, New York, 1959.

[24] TIMOSHENKO, S.P., GERE, J.M. Theory of elastic stability, 2nd Ed., McGraw-Hill Interna-tional Book Co. 1985.

[25] TIMOSHENKO, S.P., GOODIER, J. N. Theory of Elasticity, McGraw-Hill, 3rd. Ed., 1970.

[26] TSE, F.S., MORSE, I.E., HINKLE, R.T. Mechanical vibrations - theory and applications, 2nd.ed., Allyn and Bacon, Inc., Boston, USA, 1978.

[27] VOYIADJIS, G.Z., KARAMANLICHS, D. Advances in the theory of plates and shells, Elsevier,Amsterdam, 1990.

[28] YANG, T.Y. Finite element structural analysis, prentice Hall, Inc., N.J., USA, 1986.

[29] ZAHAVI, E. The Þnite element method in machine design, Prentice-hall, N.J., 1992.

[30] ZIENKIEWICZ, O.C., TAYLOR, R.L. The Þnite element method - Basic formulation and linearproblems, Vol. 1, 4th.ed., McGraw-Hill Int. Ed., 1989.

[31] ZIENKIEWICZ, O.C., TAYLOR, R.L. The Þnite element method - Solid and ßuid mechanics,dynamics and non-linearity, Vol. 2, 4th.ed., McGraw-Hill Int. Ed., 1991.

REFERÊNCIAS 74

10.0.4 Exemplo 14 Problema não-linear - Pêndulo duplo

O problema do pêndulo duplo (Figura 33) é um problema clássico que serve tanto para ilustrar umsistema não-linear quanto para a aplicação do método de Lagrange.Nota-se que, para grandes ângulos de deßexão, θ1 e θ2, o problema é não-linear, e a

determinação das equações do movimento é feita de forma mais simples pelas equações de Lagrange(que exige apenas a determinação das velocidades das massas no cálculo da energia cinética) queusando o método newtoniano, que necessita das forças e acelerações.

Figura 33: Pêndulo duplo.

Solução:Escolhem-se θ1 e θ2 como as coordenadas generalizadas do problema (posteriormente identiÞ-

caremos melhor o que sejam as coordenadas generalizadas). Então, as velocidades ao quadrado dasmassas são:

⎧⎪⎪⎪⎨⎪⎪⎪⎩v21 = L

21úθ2

1,

v22 =

⎡⎣⎛⎝L1 úθ1 sen θ1| z ab

+ L2 úθ2 sen θ2| z be

⎞⎠2

+

⎛⎝L1 úθ1 cos θ1| z de

+ L2 úθ2 cos θ2| z cd

⎞⎠2⎤⎦ . (246)

A energia cinética é simplesmente

T =1

2m1v

21 +

1

2m2v

22, (247)

e a energia potencial é calculada em relação à posição de equilíbrio θ1 = θ2 = 0:


U = m1gL1(1− cos θ1) +m2g [L1(1− cos θ1) + L2(1− cos θ2)] . (248)

Como temos apenas dois graus de liberdade, θ1 e θ2, as equações de Lagrange tornam-se:⎧⎪⎪⎨⎪⎪⎩d

dt

µ∂T

∂ úθ1

¶− ∂T

∂θ1+∂U

∂θ1= 0,

d

dt

µ∂T

∂ úθ2

¶− ∂T

∂θ2+∂U

∂θ2= 0.

(249)

O lado direito é nulo, uma vez que não temos forças aplicadas às massas além das forças mássicasassociadas à energia potencial U . Substituindo (246) e (248), e após algumas simpliÞcações, temosque as equações do movimento são:

⎧⎨⎩ L1

hL1(m1 +m2)θ1 + g(m1 +m2) sen θ1 + L2m2θ2 cos(θ1 − θ2) + L2m2

úθ2

2 sen (θ1 − θ2)i= 0,

L2m2

hL2θ2 + L1θ1 cos(θ1 − θ2) + g sen θ2 − L1 úθ21 sen (θ1 − θ2)

i= 0.

(250)Esse é um sistema não-linear de equações diferenciais, de segunda órdem, homogêneo. Caso

as oscilações sejam pequenas, (θ . 3, por exemplo), o sistema pode ser linearizado. Isso é feitoconsiderando que sen θ ≈ θ, cosθ ≈ 1, e que todo termo envolvendo produtos de incógnitas édesprezado em comparação aos demais. Issoresulta em⎧⎨⎩ L1

hL1(m1 +m2)θ1 + L2m2θ2 + g(m1 +m2) θ1

i= 0,

L2m2

hL1θ1 + L2θ2 + g θ2

i= 0.

Em forma matricial, a equação Þca

∙L21(m1 +m2) L1L2m2

L1L2m2 L22m2

¸½θ1θ2

¾+

∙L1g(m1 +m2) 0

0 L2gm2

¸½θ1θ2

¾=

½00

¾,

i.e.,M θ +K θ = 0.

11 Princípio de HamiltonEm diversos ramos da mecânica, como em análise dinâmica de rotores, é interessante utilizar umprimcípio o princípio de Hamilton, em lugar do princípio dos trabalhos virtuais, embora ambos osprincípios sejam relacionados e sejam, de fato, equivalentes em certa classe de problemas. Assim,primeiramente, serão apresentados os aspectos básicos do cálculo variacional em sua forma geral, eem seguida, o princípio de Hamilton. Será utilizado o problema de vibrações de vigas como pano defundo. Finalmente, o método de elementos Þnitos para vibrações de vigas será apresentado.

11.1 Cálculo variacional

Por conveniência, retornaremos à dedução que levou ao princípio da energia potencial mínima noCapítulo ??, mas aplicando ao problema de ßexão estática de viga. A energia potencial total deuma viga elástica em balanço (Figura 34), é


Figura 34:

V =1

2EI

Z L

0

(v)2 dx− Fv(L), (251)

onde desconsiderou-se a deformação cisalhante transversal. (A notação () signiÞca d2()/dx2). Oprincípio da energia potencial mínima aÞrma que, se a viga for submetida a qualquer deslocamentovirtual admissível3 suÞcientemente pequeno, o incremento ∆V na energia é positivo (∆V > 0). Issoé o mesmo que dizer que a conÞguração v(x) de equilíbrio é aquela que minimiza V , se v(x) for umafunção admissível.O deslocamento virtual é uma função, que pode ser posta na forma Eη(x), como na Figura 34,

onde E é um escalar unitário. Se η(x) é uma função admissível, então v(x) + Eη(x) também será. Aenergia potencial para o deslocamento v + Eη é

V +∆V =1

2EI

Z 1

0

(v+ Eη)2 dx− F (v(L) + Eη(L)). (252)

Expandindo o quadrado na integral e comparando com (251), tem-se

V +∆V =1

2EI

Z 1

0

(v)2 dx− Fv(L)| z V

+ EEI

Z 1

0

vηdx− EFη(L)| z δV

+1

2E2EI

Z 1

0

(η)2 dx| z δ2V

. (253)

Tem-se então que a variação total de energia, ∆V , é

∆V = δV +1

2δ2V. (254)

3Nesse problema, um deslocamento admissível é aquele descrito por uma função v(x) que satisfaz duas condições:(a) é uma função contínua em x ∈ (0;L); (b) satisfaz às condições de contorno essenciais (geométricas), no casov(0) = v0(0) = 0.


δV e δ2V são as denominadas primeira e segunda variação de V , e, por deÞnição, são propor-cionais a E e E2, respectivamente.Nota-se que, para qualquer F, L e η(x), caso δV 6= 0, tem-se que |δV | > 1

2

¯δ2V

¯, bastando que

E seja suÞcientemente pequeno. Então, o termo δV deÞne o sinal da variação total ∆V . Uma vezque E pode ser positivo ou negativo, mas ∆V deve ser sempre positivo, segue-se que δV = 0 é umacondição necessária para V ter um mínimo relativo. Como δ2V é sempre positivo, δV = 0 étambém uma condição suÞciente para o mínimo de V .

Que o mínimo de V , nas condições do princípio, garantem a satisfação do equilíbrio, pode serfacilmente demonstrado nesse exemplo de viga. De fato, fazendo integração por partes em δV ,eq.(253), (duas vezes sucessivamente), tem-se

δV = EEI

Z 1

0

d4v

dx4η dx+ EEIvη0|L − EEI

d3v

dx3η

¯L

− EFη(L). (255)

Uma vez que η(x) é admissível, η(0) = η0(0) = 0, essa equação se reduz a

δV = EEI

Z 1

0

d4v

dx4|z0

η dx− Eη(L)∙EId3v

dx3

¯L

+ F

¸| z

0

+ Eη0(L)EIv|L| z 0

. (256)

Para que δV = 0 para qualquer η(x), é necessário que¯¯¯d4v

dx4= 0, para ∀x ∈ (0;L),

EIv|L = 0,EId3v

dx3

¯L

+ F = 0.

(257)

Da teoria de vigas de Euler-Bernoulli, nota-se que (257)1 é a equação diferencial de equilíbrio(para carga distribuída nula), a segunda equação é a condição de momento ßetor nulo na extremidade

(M(L) = EIv(L) = 0) e a terceira condição é o esforço cortante, V (L) = EI d3vdx3

¯L= −F .

Observe que todo o procedimento demonstrado para o problema de viga pode ser repetido paraqualquer outro problema que possua um princípio de mínimo. Considere um problema mais geral,de determinar a função v(x) que minimiza a integral deÞnida

V =

Z x1

xo

F(x, v, v0, v, · · · , v(n)) dx, (258)

onde F é uma função dada, e v(x) é a função incógnita, que deve satisfazer certas condições decontinuidade e de contorno.Suponha que v(x) seja um mínimo relativo de V , Vmin. Consideremos uma outra função admis-

sível, v = v(x) + Eη(x). Então, V = Vmin + ∆V > Vmin, isto é, ∆V > 0, para E suÞcientementepequeno. Então, para F = F(x, v, v0, v), por exemplo, (estaremos aqui limitanto F a apenas duasderivadas, para simpliÞcar a exposição), temos que

∆V = V − Vmin=

Z x1

xo

F(x, v + Eη, v0 + Eη0, v+ Eη) dx−Z x1

xo

F(x, v, v0, v) dx. (259)


Para F suÞcientemente diferenciavel no intervalo, é possível expandi-lo em torno do ponto(x, v, v0, v) usando série de Taylor:

F(x, v + Eη, v0 + Eη0, v + Eη) = F(x, v, v0, v) + ∂F∂vEη +

∂F∂v0Eη

0 +∂F∂vjEη

+1

2

∙∂2F∂v2

(Eη)2 +∂2F∂v02 (Eη

0)2 +∂2F∂vj2

¡Eη¢2

2∂2F∂v∂v0E

2ηη0 + 2∂2F∂v∂vj

E2ηη + 2∂2F∂v0∂vjE

2η0η¸+O

¡E3¢.(260)

Logo, (259) Þca

∆V = E

Z x1

xo

µ∂F∂vη +

∂F∂v0η0+

∂F∂vjη¶dx| z

δV

+1

2!E2Z x1

xo

µ∂F∂vη +

∂F∂v0η0+

∂F∂vjη¶2dx| z

δ2V

+O¡E3¢. (261)

Os termos δV e δ2V são chamados primeira e segunda variação de V . Da mesma forma que nocaso da viga, para que ∆V > 0 para E suÞcientemente pequeno, é necessário que δV = 0.

Da mesma forma que no caso da viga, é possível mostrar qual equação diferencial (naquele casoa equação de equilíbrio) a minimização de V representa. Para isso, basta realizar integrações porpartes em δV . O objetivo é obter o integrando sem dirivadas em η. Então, no segundo e terceirotermos de δV , fazemos uma e duas integraçòes por partes, respectivamente. Por exemplo, para osegundo termo,

R x1xo

∂F∂v0η0 dx, deÞnem-se as variáveis auxiliares r = ∂F/∂v0 e ds = η0 dx, tal que a

integral se tgorna ∂F∂v0η

¯x1xo− R x1

xoddx

¡∂F∂v0¢η dx. Usando o mesmo procedimento para o terceiro temo

de δV , tem-se

δV = E

Z x1

xo

∙∂F∂v

− d

dx

µ∂F∂v0¶+d2

d2x

µ∂F∂v

¶¸| z

0

η dx+ E

∙η

µ∂F∂v0 −

d

dx

µ∂F∂v

¶¶+ η0

µ∂F∂v

¶¸x1xo| z

0

Para que δV seja nulo para todo η, é necessário que o termo no contorno seja nulo, o que obrigaque v(x) deve satisfazer a certas condições no contorno, as chamadas condições de Newmann, ounaturais. Em seguida, se o termo no contorno é nulo, a integral também deve sê-lo, para qualquerη admissível. Isso só será possível se o colchete for nulo em cada ponto do domínio. Assim, tem-setrês condições locais satisfeitas pela solução δV = 0:¯

¯¯∂F∂v

− d

dx

µ∂F∂v0¶+d2

d2x

µ∂F∂v

¶= 0, para ∀x ∈ (0, L),

∂F∂v

− d

dx

µ∂F∂v

¶= 0 ou η = 0 em x = 0 e x = L,

∂F∂v

= 0 ou η0 = 0 em x = 0 e x = L,

(262)

As duas últimas são as chamadas condições de contorno naturais do problema (na viga são os esforçosnas extremidades), e a primeira equação é conhecida como equação de Euler (no caso da viga, eraa equação diferencial de equilíbrio).


11.1.1 Operador variação

Em aplicações práticas de deduções envolvendo cálculo de variações, freqüentemente a função incre-mental Eη(x) é denotada por δv(x), isto é,

δv(x) = Eη(x). (263)

De fato, observa-se que δ pode ser visto como um operador, de forma similar ao operador d/dx.Quando de faz uma derivação, usa-se o operador d/dx de forma maquinal, em vez de usar a deÞniçãodf/dx = lim∆x−→0∆f/∆x. Analogamente, deseja-se determinar um procedimento automático quepermita trabalhar a variação através do operador δv(x), em vez de usar a forma Eη(x). Primeiramene,observa-se que o valor de uma função v(x) pode ser incrementado de duas formas: (a) por umdiferencial dv devido a um incremento dx (pontos AB na Figura); (b) por uma variação δv(x),dado pela adição de uma nova função: v(x) −→ v(x) + δv(x), (pontos AC na Figura). No caso devariação, a coordenada x não muda, apenas o valor da função é alterado. Na derivação, o valor dafunção muda devido à mudança de x.Como, por deÞnição, δv(x) = Eη(x), segue-se que

δv0(x) = Eη0(x),δv(x) = Eη(x), etc.

Como d(δv)/dx = Eη0, segue-se que

δdv

dx=d(δv)

dx, δ

d2v

dx2=d2(δv)

dx2, etc. (264)

Isso signiÞca que os operadores δ e d/dx são intercambiaveis.O uso do operador variação permite trabalhar com a variação de uma função de forma simbólica.

Por exemplo, a função F(x, v, v0, v) em (258) tem a primeira variação dada no primeiro integrandode (261)

δF = ∂F∂v(Eη)|zδv

+∂F∂v0(Eη0)|z

δv0

+∂F∂vj(Eη)| z δv

. (265)

O processo automatizado de realizar a operação, sem o uso da adição explícita do incrementoEη consiste em usar a mesma regra da derivada total de funções de várias variáveis: dado f =f(x, y, z), por exemplo, a derivada total é

df =∂f

∂xdx+

∂f

∂ydy +

∂f

∂zdz.

Nota-se a similaridade com a notação na equação (265), que se torna:

δF = ∂F∂vδv +

∂F∂v0δv0+

∂F∂v

δv. (266)

A k-ésima variação é deÞnida por δkF = δ ¡δk−1F¢.


11.2 Energia cinética de um sistema e coordenadas generalizadas

Considere um sistema de coordenadas retagular (X,Y, Z), Þxado a um sistema de referência newto-niano, e considere um sistema de coordenadas generalizadas (x1, x2,· · · , xn) deÞnido para um dadosistema mecânico. O problema pode ser melhor visualizado tomando o sistema discreto do pênduloduplo da Figura 33, onde θ1 e θ2 são as cordenadas generalizadas. Para uma dada partícula p, as coor-denadas retangulares são (Xp, Yp, Zp), que podem ser determinadas pelas coordenadas generalizadaspor

Xp = Xp(x1, x2, · · · , xn),Yp = Yp(x1, x2, · · · , xn),Zp = Zp(x1, x2, · · · , xn). (267)

A energia cinética do sistema é a soma das energias de cada partícula, com as velocidades dadasem componentes do sistema retangular de coordenadas:

T =1

2

NpXp=1

mp

³úX2p + úY 2p + úZ2p

´, (268)

onde mp é a massa da partícula p, e Np é o número de partículas no sistema. Diferenciando (267) notempo, obtém-se

úXp =nXq=1

∂Xp∂xq

úxq, úYp =nXq=1

∂Yp∂xq

úxq, úZp =nXq=1

∂Zp∂xq

úxq, (269)

que, substituídas em (268) resulta na energia cinética do sistema na forma:

T = 12

nXi=1

nXj=1

Mij úxi úxi. (270)

Os coeÞcientes Mij são dependentes das coordenadas generalizadas xq, e formam uma matrizsimétrica.Essa equação signiÞca que a energia cinética de um sistema com um número Þnito de graus de

liberdade pode ser obtida diretamente das velocidades generalizadas, desde que essas determinemperfeitamente as coordenadas retangulares de cada pertícula; relações do tipo (267) devem ex-istir para cada partícula, embora não precisemos ter sua representação explícita pararealizar o cálculo da energia cinética. Por exemplo, no problema do pêndulo duplo, Figura 33não precisamos determinar T em termos de ( úXp, úYp) de cada partícula, e também não precisamosconhecer as relações Xp = Xp(θ1,θ2) e Yp = Yp(θ1,θ2), p = 1, 2. Pode-se determinar diretamente Tem termos das coordenadas generalizadas, T = T (θ1,θ2), como foi feito em (247).

11.3 Princípio de Hamilton

O prncípio de Hamilton é uma generalização do princípio dos trabalhos virtuais para sis-temas discretos, pela introdução das forças de inércia através do princípio de DAlembert.O princípio de Hamilton permite a dedução da segunda lei de Newton e, para sistemas concerva-tivos, gera um princípio de mínimo. Adicionalmente, as equações de movimento de Lagrange são asequações de Euler associadas ao princípio de mínimo.


Não apresentaremos a dedução do princípio de Hamilton (ver, por exemplo [15]), mas apenasmostraremos seu enunciado.Hamilton supôs que, num dado instante to, o sistema tem uma conÞguração representada pelas

coordenadas xo de todas as partículas, e que no instante t1 as coordenadas de cada ponto são x1.O princípio de Hamilton estabelece que o movimento real das partículas durante ointervalo to −→ t1 é uma função x(t). Esse percurso é tal queZ t1

to

(δT + δW ) dt = 0, (271)

para qualquer variação δx admissível.4 T é a energia cinética total e δW é o trabalho virtualdas forças não-inerciais. Seguindo discussões apropriadas, o princípio é também estendido a sistemasmecânicos contínuos.

A segunda lei de Newton para uma partícula pode ser deduzidada do princípio de Hamilton (defatro, a dedução pode ser expandidada para sistemas discretos multi-graus de liberdade e a sistemascontínuos). Se a partícula de massa m tem posição x = (x(t), y(t), z(t)), e está submetido a umaforça F(t), o trabalho virtual num dado instante é δW = F·δx, e a energia cinética é T = 1

2m úx · úx e

δT = m úx · δ úx. Logo, pelo princípio de Hamilton,Z t1

to

(δT + δW ) dt =

Z t1

to

(m úx · δ úx+ F·δx) dt = 0.

IntegrandoZ t1

to

δT dt por partes, tem-se queZ t1

to

m úx · δ úx dt = m úx · δx|t1to −Z t1

to

mx · δx dt. Como,pelo princípio de Hamilton, δx = 0 em to e t1, segue-se queZ t1

to

(F−mx) ·δx dt = 0.

Como o princípio diz que a integral deve ser nula para qualquer δx admissível, segue-se que, a cadainstante,

F−mx = 0,o que é a equação de movimento da partícula.

11.3.1 Equação de Lagrange para sistemas concervativos

Um sistema é dito concervativo se δW = −δV , onde W é o trabalho total das cargas externase V é a energia potencial do sistema. Nesse caso, o princípio de Hamilton gera um princípio demínimo: a partir de (271) pode ser deÞnido o funcional A tal que

δA = 0, onde A =

Z t1

to

L dt, e L = T − V . (272)

4A função trajetória é tal que x(to) = xo e x(t1) = x1, isto é, seus valores são Þxos. Uma variação δx é admissívelse for contínua no intervalo (to, t1), e δx(to) = δx(t1) = 0.


L é a função Lagrangeana, ou energia Lagrangeana, e A é denominada integral de ação. Aequação variacional δA = 0 gera o seguinte princípio de mínimo:

Entre todos os trajetos x(t) do sistema, entre as conÞguraçõe xo e x1,entre os instantes to e t1, aquele que realmente ocorre é o trajeto quegera um valor estacionário para A, isto é, faz com que δA = 0.

(273)

Prova-se que as equações de Euler associadas ao funcional A do princípio de Hamilton, paraum sistema discreto com n graus de liberdade, são precisamente as equações do movimento deLagrange em termos das coordenadas generalizadas:

d

dt

∂L

∂ úxp− ∂L

∂xp= 0, para cada g.l. p. (274)

A energia potencial V pode ser escrita como

V = U −We, (275)

onde U é a energia de deformação elástica e We é o trabalho das forças externas. Numsistema discreto,

We =

NpXq=1

Fqxq, (276)

Fq e xq são as forças e coordenadas generalizadas. Os xq são medidos em relação à posição em quetodos os F 0qs são nulos.Substituindo (275) e (276) em (274), tem-se que L = T − U +We e

d

dt

∂T

∂ úxp− ∂T

∂xp+∂U

∂xp= Fp. (277)

Como xp(t) = xp(0) + up(t), isto é, a posição atual da partícula é a posição original mais o desloca-mento sofrido até o instante t, a equação (277) pode ser posta como em (109):

d

dt

∂T

∂ úup− ∂T

∂up+∂U

∂up= Fp. (278)

11.4 Sistemas contínuos - vibrações livres em vigas

A extensão do princípio de Hamilton para sistemas contínuos será apenas apresentada através de umasituação exemplo, o caso de vibrações livres de viga de Euler-Bernoulli. Para isso, considere aviga sob ßexão o plano xy como na Figura 34, com o deslocamento numa seção de cordenada x dadopela função v(x, y). A energia de deformação elástica de ßexão e cisalhamento é

U =1

2EI

Z L

x=0

(v)2 dx+1

2GAκ

Z L

x=0

(θ − v0)2 dx, (279)

Essa expressão de energia contempla a parcela associada à deformação cisalhante transversal, atravésda teoria de Timoshenko. Essa teoria considera que o cisalhamento é uniforme ao longo da seção.Uma vez que, de fato, o cisalhamento tem uma variação similar à uma parábola, para quantiÞcar


corretamente a energia de deformação, é usado o fator de cisalhamento κ Freqüentemente é usado ovalor κ = 6/5 = 1, 2. Em (279), θ é o ângulo de rotação da seção, como na Figura 34.

Figura 35: Deslocamento transversa /l v(x) e rotação θ(x) da seção transversal de uma viga.

Pelas hipóteses cinemáticas da teoria de Timoshenko, o deslocamento axial de um ponto arbitráriop, de coordenadas (x, y), é dado por u(x, y, t) = −θ(x, t). A energia cinética de uma massa diferencialdm = ρdV (ρ e dV são a densidade e o elemento diferencial de volume) localizado em (x, y) édT = 1

2dm (u2 + v2). Logo, a energia cinética da viga é

T =1

2ρ

Z L

x=0

ZA

³úv2 + úθ

2y2´dA dx.

ComoRAy2 dA = I, o momento de inércia, temos que

T =1

2ρ

Z L

x=0

³A úv2 + I úθ

2´dx. (280)

O primeiro temo é a energia cinética translacional e o segundo a energia rotacional.A teoria de Euler-Bernoulli é obtida considerando desprezivel a deformação cisalhante transvesal,

de forma que θ(x) = v0. Isso signiÞca que as energias de deformação e cinética em (279) e (280)simpliÞcam-se para

U =1

2EI

Z L

x=0

(v)2 dx e T =1

2ρ

Z L

x=0

³A úv2 + I úθ

2´dx. (281)

Para vibrações livres o trabalo das forças externas We é nulo.A energia lagrangeana vem de (272)

L = T − U = 1

2

Z L

x=0

hρA úv2 + ρI úθ

2 −EI (v)2idx. (282)


Soluções aproximadas podem ser obtidas minimizando a integral de ação³A =

R t1toL dt

´. Para

isso, podem ser usados métodos como o método de Ritz, em que a solução v(x, t) é aproximadaem uma série de funções prescritas que satisfazem as condições de contorno essenciais. Entretanto,ilustramos a seguir a aproximação via método de elementos Þnitos.

11.5 Vibrações livres de vigas pelo MEF e princípio de Hamilton

Considera-se um elemento Þnito de viga reta, de dois nós. O deslocamento transversal é aproximadonuma seção arbitrária de coordenadas x por

v(x, y) =£ψ1(x), ψ2(x), ψ3(x), ψ4(x),

¤⎧⎪⎪⎨⎪⎪⎩v1(t)θ1(t)v1(t)θ2(t)

⎫⎪⎪⎬⎪⎪⎭ , isto é,v(x, y) = ψU(t). (283)

ψ(x) é o vetor de funções de interpolação e U(t) é o vetor de deslocamentos generalizados nodaisdo elemento, dependentes do tempo. v1 e v2, θ1 e θ2 são deslocamentos e rotações nos nós 1 e 2 doelemento. Para o modelo sem cisalhamento transversal, as funçõs de interpolação são:

ψ1(x) = 1− 3³xL

´2+ 2

³xL

´3,

ψ2(x) = x

∙1− 2x

L+³xL

´2¸,

ψ3(x) = 3³xL

´2− 2

³xL

´3,

ψ4(x) = L

∙³xL

´3−³xL

´2¸. (284)

Substituindo a expansão (283) (com o uso de (284)) na energia Lagrangeana (282), e realizandoa integração,obtém-se, após extenso trabalho,

Le =1

2úUTMe úU− 1

2UTKeU. (285)

Le é a energia lagrangeana do elemento, e U = U(t) é o vetor de deslocamentos nodais do elemento.Me e Ke são as matrizes massa e rigidez do elemento, deÞnidas como: Me =Me

T +MeR, e

MeT =

ρAL

420

⎡⎢⎢⎣156 22L 54 −13L

4L2 13L −3L2156 −22L

sim. 4L2

⎤⎥⎥⎦ ,

MeR =

ρAI

120

⎡⎢⎢⎣36 3L −36 3L

4L2 −3L −L236 −3L

sim. 4L2

⎤⎥⎥⎦ ,


Ke =EI

L3

⎡⎢⎢⎣12 6L −12 6L

4L2 −6L −2L212 −6L

sim. 4L2

⎤⎥⎥⎦ . (286)

MeT e M

eR são as matrizes de inércia translacional e rotacional, respectivamente, e K

e é a mesmamatriz de rigidez do elemento Þnito usado na análise estática.As equações do movimento de Lagrange para o elemento Þnito são obtidas substuindo a energia

Lagrangeana L de (285) em (278), o que resulta em

Me úU(t) +KeU(t) = 0. (287)

essas são as equações diferenciais de movimento de elementos Þnitos para um elemento arbitrário ede viga. O sistema para a estrutura é obtido realizando as operações de sobreposição das matrizesmassa e rigidez.

12 Lista de Exercícios1.31 Mostre, através do cálculo variacional, que o segmento de reta é o caminho mais curto entre dois

pontos no plano. (Dica: Considere a Figura 36, e considere um caminho arbitrário entre os doispontos, descrito pela curva y = y(x). Um comprimento diferencial de curva é ds =

pdx2 + dy2.

Use dy = ydx, tal que o comprimento total da curva é S =R x2x1

q1 + (y)2dx. Faça δS = 0,

integre δS por partes até obter δS na formaR x2x1[·] δy, faça [·] = 0, que é uma equação diferencial.

Resolva a equação e obtenha a solução: y = ax+ b.)

Figura 36:

1.32 Considere a máquina de Atwood, mostrada na Figura 37, sem atrito, com a polia de raio r emomento de inércia de massa I. Determine a posição x(t) usando as equações de Lagrange.(Dica: V = (m1 −m2) gx, T = 1

2(m1 +m2) úx

2 + 12I úθ

2, onde θ = x/r é o ângulo de rotação

da polia. Usando a equação de Lagrange, a equação de movimento é: a x + b = 0, onde


a = (m1+m2+ I/r2), e b = m1−m2. A solução é: x(t) = −b t2/a+ vot+ xo, onde vo e xo são

velocidade e posição inicial.)

Figura 37:

1.33 Deduza e explique o método de redução de Guyan: deduza a matriz de transformação;mostre a deÞnição do sistema reduzido. Qual a motivação para reduzir as matrizes?

1.34 Para o método de redução de Guyan, qual o processo de seleção automática das equaçõesretidas (que identiÞcam os graus de liberdade master?

1.35 Explique o processo de expanssão no método de redução de Guyan, e no de decomposiçãomodal?

1.36 Uma vez que se tenha um programa para solução de problemas dinâmicos lineares com car-regamentos harmônicos, sen ωt e cosωt, como ele poderia ser alterado para resolver problemasde carregamento periódicos, através de séries de Fourier?

Documents

Analise Dinamica Pelo Metodo Dos Elementos Finitos