DinaFEM

Prévia do material em texto

Análise Dinâmica pelo Método de Elementos
Finitos
Prof. Paulo de Tarso R. Mendonça, Ph.D.
Universidade Federal de Santa Catarina - UFSC
Grante - Depto. Engenharia Mecânica, CP 476,
Florianópolis, SC. CEP 88040-900,
mendonca@grante.ufsc.br
23 de Agosto de 2006
Resumo
O estudo que relaciona as forças que atuam sobre um corpo com o movimento, tanto do corpo
como 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 do
movimento 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ção
de equilíbrio. Este texto lida com a simulação numérica por elementos Þnitos de vibrações em
corpos sólidos. Nas primeiras seções desse capítulo e no próximo, apresentamos uma revisão
dos fundamentos de vibrações mas deve-se ter claro que este não é um substitutivo a um curso
formal no assunto. O objetivo destes tópicos aqui consiste apenas em homogeneizar o material
e a notação de forma a permitir o tratamento do problema pelo método dos elementos Þnitos.
Conteúdo
1 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 12
5.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 21
7.1 Princípio de D’Alembert . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
7.2 Princípio dos Trabalhos Virtuais em Barras . . . . . . . . . . . . . . . . . . . . . . . 22
1
CONTEÚDO 2
7.3 Matrizes de Elementos Finitos de Barras . . . . . . . . . . . . . . . . . . . . . . . . . 23
7.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 28
8.1 Vibrações Livres Não-Amortecidas . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30
8.2 Propriedades dos Autovetores e Autovalores . . . . . . . . . . . . . . . . . . . . . . . 34
8.2.1 Ortogonalidade . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34
8.2.2 Normalização e Ortonormalidade . . . . . . . . . . . . . . . . . . . . . . . . . 35
8.2.3 Exemplo 3 — Freqüências Naturais . . . . . . . . . . . . . . . . . . . . . . . . . 37
8.2.4 Exemplo 4 — Modos Naturais . . . . . . . . . . . . . . . . . . . . . . . . . . . 38
8.2.5 Exemplo 5 — Normalização de Autovetores . . . . . . . . . . . . . . . . . . . . 39
8.2.6 Exemplo 6 — Solução Analítica de Vibrações Livres em Barra . . . . . . . . . . 40
8.2.7 Autovetores Linearmente Independentes . . . . . . . . . . . . . . . . . . . . . 42
8.3 Análise Modal para Excitação Inicial - Sistema não-Amortecido . . . . . . . . . . . . 43
8.3.1 Exemplo 7 — Resposta para Deslocamento Inicial pelo MEF . . . . . . . . . . 47
8.3.2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49
8.4 Análise Modal Geral . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51
8.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 56
9.1 Um Grau de Liberdade . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56
9.2 n-Graus de Liberdade — Elementos Finitos . . . . . . . . . . . . . . . . . . . . . . . . 57
9.3 Métodos Experimentais . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58
9.3.1 Exemplo 10 — Determinação Experimental de Amortecimento . . . . . . . . . 59
9.4 Método Analítico 1 — Rayleigh . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60
9.5 Método Analítico 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63
9.5.1 Exemplo 11 — Determinação Experimental da Matriz de Amortecimento . . . . 63
9.5.2 Exemplo 12 — Vibração Amortecida de Barra sob Deslocamento Inicial . . . . 64
9.5.3 Exemplo 13 — Vibração Forçada Amortecida pelo MEF . . . . . . . . . . . . . 65
10 Lista de Exercícios 66
10.0.4 Exemplo 14 — Problema não-linear - Pêndulo duplo . . . . . . . . . . . . . . . 74
11 Princípio de Hamilton 75
11.1 Cálculo variacional . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75
11.1.1 Operador variação . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79
11.2 Energia cinética de um sistema e coordenadas generalizadas . . . . . . . . . . . . . . 80
11.3 Princípio de Hamilton . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80
11.3.1 Equação de Lagrange para sistemas concervativos . . . . . . . . . . . . . . . . 81
11.4 Sistemas contínuos - vibrações livres em vigas . . . . . . . . . . . . . . . . . . . . . . 82
11.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 Liberdade
Consideremos 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 de
corpos 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ão
do comportamento deste sistema é a peça básica sobre a qual são construídas as formulações dos
sistemas 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 é considerada
sem massa e o amortecedor linear da constante c é considerado sem massa ou rigidez. Estes são
evidentemente idealizações. A rigidez da mola signiÞca que Fk = k δk onde δk é o deslocamento entre
as extremidades da mola provocado pela força Fk. O amortecedor é tal que, Fd = c úδd, onde úδd é a
velocidade 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, obtemos
a 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ção
corresponde à situação onde a mola está descarregada. Substituindo as expressões para as forças
obtemos a equação do movimento na forma:
mx¨(t) + c úx(t) + k x(t) = F (t). (2)
Paulo de Tarso R Mendonça - Grante - EMC - UFSC, 2004 4
Essa é uma equação diferencial linear, ordinária de coeÞcientes constantes m, c e k, que deÞnem
as 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 de
forç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 movimento
durante 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, a
solução deste problemafornece subsídeos para a solução de problemas excitados, o que constitui a
outra 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, ω
2
n =
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 do
deslocamento 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 =
A1
cosφ
=
A2
senφ
, tanφ =
A2
A1
(6)
As novas constantes A e φ tem signiÞcado físico mais evidentes que A1 e A2, são a amplitude
e â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 com
a 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
ωn
senωnt. (8)
A amplitude e o ângulo de fase são
Paulo de Tarso R Mendonça - Grante - EMC - UFSC, 2004 5
A =
s
x2o +
µ
υo
ωn
¶2
, tanφ =
υo
xoωn
. (9)
Além das formas (4) e (5) existe ainda uma terceira forma para a solução do problema (3), dada
por
x(t) = A est, (10)
onde A e s são constantes a serem determinadas. Derivando (10) é possível ver que ela realmente
satisfaz a equação diferencial (3). Fazendo a substituição obtemos
s2A est + ω2nA e
st = 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) = A1e
iωn +A2e
−iωn. (14)
3 Vibrações Livres de Sistemas Amortecidos
Quando o amortecimento do sistema não é nulo, a equação de movimento (2) é reescrita para uma
forma análoga a (3):
x¨(t) + 2ζωn úx(t) + ω
2
n x(t) = 0 (15)
onde
ζ =
c
2mωn
. (16)
é o chamado quociente de amortecimento viscoso. ζ tem signiÞcado físico deÞnido e será visto
na 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+ ω
2
n = 0. (17)
Dois valores de s satisfazem esta equação:
s1
s2
¾
=
µ
−ζ ±
q
ζ2 − 1
¶
ωn. (18)
Paulo de Tarso R Mendonça - Grante - EMC - UFSC, 2004 6
Cada raiz produz uma solução em (10). Da teoria de equações diferenciais lineares, temos que a
solução do problema é uma combinação linear das soluções:
x(t) =
h
A1e
√
ζ2−1ωnt +A2e−
√
ζ2−1ωnt
i
e−ζωn t. (19)
Os coeÞcientes A1 e A2 podem ser complexos. Claramente as raízes s1 e s2 podem ser reais ou
complexas 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 = ωn
p
1− ζ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 parte
imaginária de x(t). De (23) nota-se que amas são deÞnidas pela mesma base de funções, sen ωdt e
cosωdt. (Que essas funções são solução pode ser provado simplesmente substituindo-as na equação
diferencial.) 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,com
comprimento A e ângulo φ em relação ao vetor eiωdt. Assim, em vez da forma (20), a solução pode
ser 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ária
tem-se:
Paulo de Tarso R Mendonça - Grante - EMC - UFSC, 2004 7
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φ =
c2
c1
A2 = c21 + c
2
2.
(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Þcado
das constantes, mas elas permanecem sempre duas, ou A1 e A2, ou A e φ, ou c1 e c2. Cada uma destas
diferentes formas para a solução é mais adequada a diferentes tipos de interpretações e manipulações
algébricas, algumas das quais faremos uso ao longo do texto. Esta última forma é obtida usando as
deÞ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/ωd
e a solução Þca
x(t) =
vo
ωd
e−ζωnt senωdt, ωd = ωn
p
1− ζ2. (30)
Essa é a resposta do sistema à velocidade inicial vo, denominada solução transiente do
problema. O termo “transiente” se refere ao fato de que ela consiste em uma função periódica, o
seno, que por si tem amplitude constante igual a 1 para todo t > 0. Porém o amortecimento no termo
exponencial faz com que o fator multiplicando seno decresça ao longo do tempo. Assim, as oscilações
vão diminuindo de amplitude, como ilustrado na Figura 3. (Ali mostramos o caso de deslocamento
inicial, 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 inicial
x(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)
Paulo de Tarso R Mendonça - Grante - EMC - UFSC, 2004 8
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ção
com uso da segunda, o que resulta em
tanφ =
ζωn
ωd
e A =
xo
cosφ
. (32)
Tomando ωd = ωn
p
1− ζ2 de (30), tem-se
tanφ =
ζp
1−ζ2
. (33)
De trigonometria temos que, da relação acima, cosφ =
p
1− ζ2. (considere um triângulo retângulo
com catetos ζ e
p
1− ζ2). Então, da segunda das eqs. (32),
A =
xop
1− ζ2
(34)
Então a solução do problema de vibração livre amortecida com x(0) = xo e úx(0) = 0 é obtida
levando (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) + ω
2
n x(t) = 0,
x(0) = uo,
úx(0) = vo.
(36)
Paulo de Tarso R Mendonça - Grante - EMC - UFSC, 2004 9
Uma vez que o problema é linear podemos simplesmente sobrepor a solução (30) obtida para
x(0) = 0 e úx(0) = vo, com a solução (35) obtida para x(0) = uo e úx(0) = 0. Então a solução
completa é
x(t) =
∙
uo√
1−ζ2
cos (ωdt− φ) + voωd senω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
ωd
senω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ζωn
xoωd
.
4 Carregamento Harmônico
Consideremos agora a solução particular do problema (2). O caso mais simples de carregamento é o
chamado “carregamento harmônico”, isto é, um que varia segundo um seno ou cosseno ao longo do
tempo:
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) + ω
2
n x(t) = ω
2
nA 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) obtemos
a 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 fazer
o 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.
Paulo de Tarso R Mendonça - Grante - EMC - UFSC, 2004 10
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 carregamento
considerado é 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+ ψ), ou
na parte imaginária, Im(f(t)) = sen (ω t+ ψ). Da solução complexa tomaremos também a parte
real 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 de
cosseno. Na presente dedução, pretende-se aumentar a generalidade, permitindo que o carregamento
tenha uma fase, o que permite que o seno (ou cosseno) se inicie num valor diferente de sua posição
natural. 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) + ω
2
n x(t) = ω
2
nA e
i(ωt+ψ). (47)
Paulo de Tarso R Mendonça - Grante - EMC - UFSC, 2004 11
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 forma
X = 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ω) = ω
2
nA sen (ψ + φ),
que podem ser resolvidas para a amplitude X¯ e para a fase φ:
tan(ψ + φ) =
2ζω/ωn
1− (ω/ωn)2
e 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 a
soluçã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 conjugado
complexo 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)
Paulo de Tarso R Mendonça - Grante - EMC - UFSC, 2004 12
onde β = ω/ωn. O termo real, o primeiro dentro do colchete, corresponde à solução transiente
amortecida 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 resulta
na solução:
x(t) =
A
(1− β2) [cosωt− cosωnt]. (56)
5 Resposta a Carregamentos não-Periódicos
Estamos interessados em obter a resposta do sistema a um carregamento qualquer como o ilustrado
na Figura 5. Vários métodos existem para estimar a solução deste problema, e nos concentraremos
aqui 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 ilustrada
na 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)
Paulo de Tarso R Mendonça - Grante - EMC - UFSC, 2004 13
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 pseudo
função (ou função generalizada), denominadafunção Delta de Dirac δ (t− a) , como
δ (t− a) = lim
7→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 (decorrente
de (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 integral
igual ao valor desta função no ponto a. Esta característica pode ser entendida com a ajuda da Figura
7. 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)lim
7→0
ga(t) dt = lim
1
E
Z a+7
a
G(t) dt = lim
1
E
G(a)E = G(a),
que resulta (62). Lembremos que as operações acima são apenas formais, e que resta provar algumas
delas, como a passagem do limite para fora da integral. Aquela relação constitui a principal utilidade
da função Delta de Dirac.
Consideremos agora um outro conceito, o impulso, deÞnido como
Impulso −→ I = R∞
t=0
F (t) dt. (63)
Paulo de Tarso R Mendonça - Grante - EMC - UFSC, 2004 14
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 ga
ilustrada na Figura 6 quando E → 0. Apesar da duração da força tender a zero, desejamos que seu
impulso, 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:
lim
7→0
R 7
0
ˆFδ(t) dt = ˆF,
lim
7→0+
R 7
0
mx¨ dt = lim
7→0+
(m úx)|70 = lim
7→0+
m ( úx (E)− úx (0)) = m úx(0+),
lim
7→0+
R 7
0
c úx dt = lim
7→0+
(cx)|70 = lim
7→0+
c (x (E)− x (0)) = 0,
lim
7→0+
R 7
0
kx 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 apenas
x (0) = 0, mas também x (0+) = 0.
O resultado das integrações em (66) é
úx
¡
0+
¢
=
ˆF
m
. (67)
Paulo de Tarso R Mendonça - Grante - EMC - UFSC, 2004 15
Isto mostra que a aplicação de uma força impulsiva ˆF produz uma mudança instantânea de
velocidade. Isto é interessante porque estamos interessados em obter a resposta do sistema a ˆF . Mas
devido a (67) podemos usar a solução que já dispomos para a resposta do sistema a uma velocidade
inicial. Usamos então a eq. (37) com vo = ˆF/m e uo = 0:
x (t) =
ˆF
mωd
e−ζωnt sinωdt, ωd = ωn
¡
1− ζ2¢1/2 , t > 0. (68)
Fazendo ˆF = 1 temos a chamada resposta impulsiva
h(t) =
1
mωd
e−ζω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 instante
inicial τ = 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 instante
t = 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 por
uma seqüência de retângulos como aquele hachurado, cada um iniciando num instante τ , de duração
4τ , 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 o
instante do impulso e o instante t. Isto porque, obviamente, o corpo não sofre nenhum efeito do
impulso antes dele ter sido aplicado, isto é em t < τ . Mais ainda, a resposta do sistema no instante
t depende apenas do lapso de tempo desde a aplicação do impulso, isto é, da extensão de tempo
decorrido t − τ . Então, se o impulso foi aplicado no instante τ , a solução em t pode ser obtida
simplesmente usando a solução (70) substituindo t por (t− τ), isto é
4x(t) = F (τ) 4τ h (t− τ) = F (τ)
mωd
e−ζωn(t−τ) senωd (t− τ)4τ . (71)
Mas observe que no instante t o valor do deslocamento não é apenas devido a este intervalo de
carregamento aplicado entre τ e τ + 4τ . Existem também parcelas provenientes dos impulsos de
duração 4τ aplicados desde o instante 0 até t que compõem a curva F (t). Então a resposta total em
t é
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)
Paulo de Tarso R Mendonça - Grante - EMC - UFSC, 2004 16
Esta integral aparece em diversas áreas das ciências físicas e é logicamente objeto de estudo
matemático em busca de suas propriedades. É a chamada integral de convolução. Uma das
propriedades 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 é chamada
solução particular ou solução de regime permanente. Lembremos que esta é apenas parte da
soluçã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
−ζωntp
1− ζ2
cos (ωdt− φ) + vo
ωd
e−ζωnt senωdt+
1
mωd
Z t
0
F (τ)e−ζωn(t−τ) senωd(t− τ) dτ,
ωd = ωn
p
1− ζ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 um
carregamento 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τ.
Paulo de Tarso R Mendonça - Grante - EMC - UFSC, 2004 17
0 2 4 6
t
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 movimento
da massa não é oscilante em torno do ponto de equilíbrio x = 0, mas sofre um movimento oscilante
onde 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.0
t
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), usando
os 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. Observe
a deÞnição em (77). Consideremos o caso em que o carregamento seja dado como uma combinação
de dois outros carregamentos:
Paulo de Tarso R Mendonça - Grante - EMC - UFSC, 2004 18
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 para
t < t1 e F1(t) = senωt e F2(t) = senω(t − t1) para t ≥ t1. Podemos deÞnir uma função F (t) como
F (t) = F1(t)−F2(t), como ilustrado na Figura 10c. Da integral de convolução temos a resposta para
o 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 de
uma função arbitrária, em geral ela é realizada numericamente. Entretanto, algumas propriedades
da integral de convolução permitem alargar um pouco o leque de situações passíveis de ter solução
analítica. Por exemplo, suponha que se tenha conseguido obter a solução para um carregamento
F1(t) como mostrado na Figura 10(a) e (b). Se transladarmos F1(t) em t1 e deÞnirmos assim a
função F2(t), a solução x2(t) é a solução x1(t) transladada, isto é, x2(t) = x1(t − t1) para t > t1 e
x2(t) = 0 para t < t1. Assim, a solução associada a F = aF1 + b F2 (Figura 10c)é a combinação das
duas soluções:
x(t) = a x1(t) + b x2(t). (83)
Paulo de Tarso R Mendonça - Grante - EMC - UFSC, 2004 19
Esta possibilidade de combinação é resultante da linearidade da equação do movimento usada. Estas
operações também necessitam que as condições iniciais uo e vo sejam combinadas da mesma forma
atravé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 como
descrito nas seções acima, composto um bloco rígido de massa m, ligado a uma base por uma mola
e um amortecedor. Estamos interessados principalmente em determinar o comportamento dinâmico
dos 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ções
acima é 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 onde
o 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 componentes
e 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, a
maneira hoje padrão de se tratar destes problemas consiste em abrir mão do desejo de obter uma
soluçã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 uma
coleção de massas discretas unidas por molas e amortecidedores entre si. A forma de realizar este
processo de discretização não é óbvio, e existem diversos métodos, dentre os quais o próprio método
que 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 n
massas discretas. Cada uma dessas massas pode ser considerada um corpo rígido, de forma que
podemos aplicar a ela a segunda lei de Newton. A Figura 11c representa um diagrama de corpo livre
de uma massa genéricami. Sobre ela atuam uma força externa Fi(t) e as forças internas provenientes
dos 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 constante
Ci e Ci+1 dos amortecedores. Pela segunda lei de Newton, a resultante de todas estas forças deve ser
igual à força de inércia mix¨i. 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) = miu¨i. (84)
Podemos rearranjar os termos colocando a parte conhecida, a força externa Fi(t), do lado direito:
miu¨i − 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:
Paulo de Tarso R Mendonça - Grante - EMC - UFSC, 2004 20
(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)1 u (t)i+1u (t)i nu (t)
i+1mmi-1
k 1 ik
i+1c
nk
1c ic nc
k i+1
mi
ik 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 corpo
livre da massa mi indicando força externa aplicada, força de inércia, forças elásticas de mola e forças
de amortecimento.
⎡⎢⎢⎢⎢⎢⎣
m1
m2
m3
. . .
mn
⎤⎥⎥⎥⎥⎥⎦
⎡⎢⎢⎢⎢⎢⎢⎢⎣
u¨1
u¨2
u¨3
u¨4
...
u¨n
⎤⎥⎥⎥⎥⎥⎥⎥⎦
+
⎡⎢⎢⎢⎢⎢⎣
(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
. . .
⎤⎥⎥⎥⎥⎥⎥⎥⎦
⎡⎢⎢⎢⎢⎢⎢⎢⎣
u1
u2
u3
u4
...
un
⎤⎥⎥⎥⎥⎥⎥⎥⎦
=
⎡⎢⎢⎢⎢⎢⎢⎢⎣
F1(t)
F2(t)
F3(t)
F4(t)
...
Fn(t)
⎤⎥⎥⎥⎥⎥⎥⎥⎦
(86)
Paulo de Tarso R Mendonça - Grante - EMC - UFSC, 2004 21
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 as
chamadas matriz massa ou de inérciaM, matriz de amortecimento C e a matriz de rigidez K, todas
simétricas. Embora neste exemplo a M seja diagonal, de forma geral isto não é assim.
7 Elementos Finitos em Dinâmica
Nas próximas seções trataremos de métodos para determinar ou estimar a solução do problema
algébrico (87). Antes disso daremos uma amostra do processo geral de como aquelas matrizes são
determinadas pelo método de elementos Þnitos na discretização de um corpo sólido, isto é, um sistema
contí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 trabalho
virtual (PTV), onde no caso dinâmico fazemos uma alteração em seu enunciado pelo uso do chamado
princípio de D’Alembert, descrito a seguir. Um outro procedimento a ser apresentado, além do
PTV, é a obtenção das equações de matriciais de movimento pelo uso das equações de movimento
de Lagrange. Estas equações, porém, são apenas uma forma derivada do mesmo PTV aplicado à
dinâmica.
7.1 Princípio de D’Alembert
Julgando-se apenas peloseu enunciado, este princípio é de uma simplicidade enorme. Sua utilidade
entretanto é também enorme na engenharia. Considere a equação do movimento de uma partícula
de massa m, dada pela segunda lei de Newton:
nX
i=1
Fi +mb = ma, (88)
isto é, a resultante de todas as n forças externas aplicadas Fi, incluindo a força do corpomb, deve ser
igual a força da inércia, dada pela massa vezes a aceleração a desenvolvida pela massa m; b é uma
força de corpo por unidade de massa. Quando as forças são tais que a aceleração é nula, as forças
estão em equilíbrio e esta equação é chamada equação de equilíbrio. Obviamente, o tratamento
de problemas de equilíbrio é mais simples que o tratamento de problemas dinâmicos. D’Alembert,
de certa forma, realizou uma operação bastante “simples”. Ele transferiu a força de inércia do lado
direito de (88) e passou-a para o lado esquerdo obtendo
nX
i=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 o
que existe desenvolvido para os problemas estáticos pode ser adaptado para os problemas dinâmicos.
O princípio de D’Alembert então aÞrma que as forças de inércia podem ser incorporadas às forças de
corpo e o problema pode ser tratado como uma equação estática, embora, agora, as forças de corpo
associadas à inécia sejam desconhecidas, por incorporarem as acelerações.
Paulo de Tarso R Mendonça - Grante - EMC - UFSC, 2004 22
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 tipos
de componentes para o comportamento estático. Com o uso do princípio de D’Alembert as mesmas
expressõ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
dˆu
dx
dx−A
Z L
o
bˆu 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 extremindade
x = L. Considere um elemento diferencial de volume numa barra, como ilustrado na Figura 12. A
massa 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 um
elemento diferencial numa barra sob solicitação dinâmica.
Observe na Figura 13 o comportamento dinâmico de uma barra sob carga axial. O elemento
diferencial inicialmente encontra-se a uma distância x da origem. Num outro instante t a posição
daquela porção de material encontra-se a uma posição X da origem. A posição atual é função da
posiçã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 deslocamento
sofrido. Como a posição inicial não se altera, diferenciando (91) temos a velocidade e a aceleração:
Paulo de Tarso R Mendonça - Grante - EMC - UFSC, 2004 23
⎧⎪⎪⎨⎪⎪⎩
ú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 diferencial
de 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 corpo
assiciada à inércia. Assim, a força de inércia pode ser incluída no P.T.V. da equação (90) substituindo
a força de corpo estática b(x) por (b(x, t)−binércia), isto é, por (b(x, t)− ρu¨(x, t)), levando à expressão
do PTV dinâmico
AE
Z L
o
∂u(x, t)
∂x
dˆu(x)
dx
dx−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 extremidade
e forças de corpo b distribuídas ao longo de sua extensão, como ilustrado na Figura 12, mas agora as
forç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 carregamentos
tem valores deÞnidos b(x, t¯) e f(t¯). Tem-se então o PTV estático neste instante, onde se deve buscar
a solução também estática, u(x, t¯). O tratamento por elementos Þnitos consiste ainda em discretizar o
corpo 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 − x
Le
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 deslocamento
axial no elemento e é então interpolado por
u(x, t) =
2X
i=1
ui(t)ϕ
e
i (x). (96)
Observe que os deslocamentos nodais em (96), os u0is, agora variam com o tempo. De fato, essa
expanssã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) =
2X
i=1
ˆuiϕ
e
i (x). (97)
Entretanto, os valores nodais ˆui da função peso não variam no tempo. De (96), temos claramente a
aproximação da aceleração no elemento:
Paulo de Tarso R Mendonça - Grante - EMC - UFSC, 2004 24
u¨(x, t) =
2X
i=1
u¨i(t)ϕ
e
i (x). (98)
Levamos essas expressões ao P.T.V. da eq. (94):
AE
Z
e
Ã
2X
i
ui(t)
dϕei
dx
(x)
!
| {z }
u(x,t)
Ã
2X
i
ˆui
dϕei
dx
(x)
!
| {z }
ˆu(x)
dx−
A
Z
e
Ã
b(x, t)− ρ
2X
i
u¨i(t)ϕ
e
i (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 dois
somatórios na primeira integral e transfere-se o termo independente de x para fora da integral (os
coeÞcientes ui(t) e ˆui). Primeiramente obtém-se a forma expandida dos somatórios
AE
R
e
µ
u1(t)
dϕe1
dx
+ u2(x)
dϕe2
dx
¶
| {z }
u(x,t)
µ
ˆu1
dϕe1
dx
+ u2
dϕe2
dx
¶
| {z }
ˆu(x)
dx−A R
e
b(x, t) (ˆu1ϕ
e
1 + ˆu2ϕ
e
2) dx
−ρA R
e
(u¨1(t)ϕ
e
1 + u¨2(t)ϕ
e
2) (ˆu1ϕ
e
1 + ˆu2ϕ
e
2) dx−Af(t)ˆu2 = 0 ∀ˆu1, ˆu2.
(100)
Como esta igualdade deve ser satisfeita para quaisquer valores de ˆu1 e ˆu2, podemos fazer ˆu1 = 1
e ˆu2 = 0 e obter uma equação algébrica. Em seguida pode-se fazer ˆu1 = 0 e ˆu2 = 1, obtendo uma
segunda equação:
=⇒
∙
AE
Z
e
dϕe1
dx
dϕe1
dx
dx
¸
| {z }
Ke11
u1(t) +
∙
AE
Z
e
dϕe1
dx
dϕe2
dx
dx
¸
| {z }
Ke12
u2(t)+
∙
ρA
Z
e
ϕe1ϕ
e
1 dx
¸
| {z }
Me11
u¨1(t) +
∙
ρA
Z
e
ϕe1ϕ
e
2 dx
¸
| {z }
M212
u¨2(t) = A
R
e
b(x, t)ϕe1 dx,
=⇒
∙
AE
Z
e
dϕe2
dx
dϕe1
dx
dx
¸
| {z }
Ke21
u1(t) +
∙
AE
Z
e
dϕe2
dx
dϕe2
dx
dx
¸
| {z }
Ke22
u2(t)+
∙
ρA
Z
e
ϕe2ϕ
e
1 dx
¸
| {z }
Me21
u¨1(t) +
∙
ρA
Z
e
ϕe2ϕ
e
1 dx
¸
| {z }
Me22
u¨2(t) = A
R
e
b(x, t)ϕe2 dx+Af(t).
(101)
Pode-se reconhecer as primeiras duas integrais em cada equação como os termos da matriz de
rigidez 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 ou
matriz 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
Paulo de Tarso R Mendonça - Grante - EMC - UFSC, 2004 25
(
Ke11u1(t) +K
e
12u2(t) + M
e
11u¨1(t) +M
e
12u¨2(t) = Fe
1 (t),
Ke21u1(t) +K
e
22u2(t) + M
e
21u¨1(t) +M
e
22u¨2(t) = F
e
2 (t),
(102)
ou, em notação matricial,
Meu¨e(t) +Keue(t) = Fe(t). (103)
Comparando com o caso estático vemos que agora existe uma força de inércia, Meu¨e(t), e que
agora o carregamento pode variar no tempo. Mas mesmo que o carregamento seja estático, a presença
do 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
Z
e
dϕei
dx
dϕej
dx
dx,
Meij = ρA
Z
e
ϕeiϕ
e
j dx,
F ei (t) = A
Z
e
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 acima
podem 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ção
analítica da matriz de inércia é a seguinte¯¯¯¯
¯¯¯¯
¯¯¯¯
¯¯¯
Me11 =
ρA
L2e
Z xe2
xe1
(Le − x)2 dx = ρALe
3
,
Me21 = M
e
12 =
ρA
L2e
Z xe2
xe1
(Le − x)x dx = ρALe
6
,
Me22 =
ρA
L2e
Z xe2
xe1
x2 dx =
ρALe
3
.
(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 =
ρALe
6
∙
2 1
1 2
¸
; Fe = Af
∙
0
1
¸
. (106)
Como no caso estático, essas são apenas as matrizes de um elemento, e devem ser rotacionadas e
sobrepostas nas matrizes globais para gerar as equações discretas de movimento que representam
o 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 simplesmente
resolvida por inversão de matriz. É um sistema de equações diferenciais e deve ser integrado para
dar a resposta do sistema, após a aplicação das condições de contorno e condições iniciais do sistema.
Paulo de Tarso R Mendonça - Grante - EMC - UFSC, 2004 26
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 das
equações de Lagrange. Considere que podemos expressar a energia de deformação U e a energia
ciné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 que
tentando 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
¶2
dx. (110)
A energia cinética pode ser obtida da seguinte forma. Lembre que a energia cinética de uma
massa pontual m é, por deÞnição, Ec = mv2/2. Agora considere o elemento diferencial de barra de
comprimento dx das Figuras 12 e 13. Este elemento tem massa diferencial dm = ρAdx e velocidade
axial úu(x, t). Então sua energia cinética é úu(x, t)2dm/2, isto é, ρA úu(x, t)2dx/2. A energia cinética
da 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 em
elementos, o que signiÞca simplesmente particionar o intervalo de integração nas deÞnições de W
e T em uma soma de integrais realizadas em cada elemento. Em cada elemento interpolamos o
deslocamento e velocidade usando (96). Então as energias em cada elemento se tornam:
W e =
AE
2
Z
e
µ
u1(t)
dϕe1
dx
+ u2(t)
dϕe2
dx
¶2
dx,
T e =
ρA
2
Z
e
( úu1(t)ϕ
e
1 + úu2(t)ϕ
e
2)
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
12
Le
∙
úu1(t)
úu2(t)
¸T ∙
2 1
1 2
¸ ∙
úu1(t)
úu2(t)
¸
.
(113)
1O expoente T indica transposto de vetor ou matriz.
Paulo de Tarso R Mendonça - Grante - EMC - UFSC, 2004 27
Se compararmos estas expressões a (106) vemos que as matrizes acima são proporcionais a rigidez
e massa do elemento respectivamente. DeÞnindo o vetor de deslocamentos nodais do elemento
como ue(t) = [u1(t);u2(t)]
T , a equação acima pode ser posta na forma½
W e = 1
2
ueT (t)Ke ue(t),
T e = 1
2
úueT (t)Me úue(t).
(114)
Tendo então as expressões discretizadas para as energias interna e cinemática, podemos passar
ao uso das equações de Lagrange. Para um dado elemento a equação de Lagrange em (109) reduz-se
a duas equações, para i = 1 e i = 2, correspondentes aos dois graus de liberdade do elemento. As
equaçõ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
[Ke11u
2
1 +K
e
12u1u2 +K
e
21u1u2 +K
e
22u
2
2,
T e = 1
2
[Me11 úu
2
1 +M
e
12 úu1 úu2 +M
e
21 úu1 úu2 +M
e
22 úu2 úu2.
Fazendo as derivações indicadas em (115) e recolocando o resultado em forma matricial obtemos∙
Me11 M
e
12
Me21 M
e
22
¸ ∙
u¨1(t)
u¨2(t)
¸
+
∙
Ke11 K
e
12
Ke21 K
e
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ípio
de D’Alembert 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 1
0 1 2
⎤⎦⎡⎣ u¨1(t)u¨2(t)
u¨3(t)
⎤⎦+ 2EA
L
⎡⎣ 1 −1 0−1 2 −1
0 −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 0
1 4 1 0
0 1 4 1
0 0 1 2
⎤⎥⎥⎦
⎡⎢⎢⎣
u¨1(t)
u¨2(t)
u¨3(t)
u¨4(t)
⎤⎥⎥⎦+ 3EAL
⎡⎢⎢⎣
1 −1 0 0
−1 2 −1 0
0 −1 2 −1
0 0 −1 1
⎤⎥⎥⎦
⎡⎢⎢⎣
u1(t)
u2(t)
u3(t)
u4(t)
⎤⎥⎥⎦ =
⎡⎢⎢⎣
F1(t)
F2(t)
F3(t)
F4(t)
⎤⎥⎥⎦ .
Paulo de Tarso R Mendonça - Grante - EMC - UFSC, 2004 28
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 Modal
Considere o caso em que normalmente se considera como estático, onde o carregamento não varia
com o tempo. Devemos lembrar que certamente houve um período inicial em que a carga teve que
ser 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 serem
desprezadas e aanálise pode ser feita como estática, sem o primeiro termo de (107), onde força e
deslocamento 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 de
um 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ão
oscilaria várias vezes sobre a suspensão. Devido ao amortecimento, essas oscilações gradualmente se
reduziriam 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 problema
dinâmico. Por exemplo, considere um carregamento cíclico com baixa freqüência. Novamente, se a
freqüência de carregamento é baixa, as acelerações, que também são cíclicas, serão baixas. Isto pode
ser visto de (48). Então as acelerações do sistema podem ser desprezadas em (107), resultando num
sistema 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 a
freqüência de excitação for menor que aproximadamente um terço da menor freqüência natural do
Paulo de Tarso R Mendonça - Grante - EMC - UFSC, 2004 29
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 s
ou to > 10Tmax, onde Tmax = 2π/ω1 e ω1 é a menor freqüência natural); (b) Aplicação com resposta
dinâmica, usuamente para tempos de carregamento da ordem de to < 0, 01 s.
sistema, isto é,
ω . ω1
3
, (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ção
de onda ocorrem em situações de impacto ou de explosões ou de acústica entre outros, onde tanto o
carregamento quanto a resposta são de alta freqüência e o período de duração da análise é em geral
curto, da ordem de um período da onda que cruza a estrutura. Por outro lado, quando a freqüência
de carregamento não é alta, no sentido de que é da mesma ordem, ou apenas algumas vezes maior
que 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 e
modos 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 abaixo
da primeira freqüência natural, ou pelo menos interessa evitar que a freqüência de excitação Þque
pró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 modos
naturais, comentados no parágrafo anterior, enquanto o de integração direta faz uma discretização
de diferenciais Þnitos no tempo na equação diferencial do movimento, (107), e faz uma integração
numérica. O método de análise modal é um método baseado fundamentalmente na linearidade do
sistema. Por outro lado, quando o sistema físico é modelado matematicamente, por elementos Þnitos
por exemplo, levando em conta efeitos nãolineares, como plasticidade emmetais, grandes deformações
Paulo de Tarso R Mendonça - Grante - EMC - UFSC, 2004 30
como em processos de conformação, ou grandes deslocamentos, o processo adequado a ser usado é o
de integração direta no tempo das equações de movimento, embora existam formas de circunscrever
as 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 processo
de 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 foram
já discretizados. Este tipo de equação é bastante conhecida e estudada em matemática, uma vez
que toda uma série de fenômenos físicos é modelada por sistemas de equações diferenciais ordinárias
desse tipo. Uma classe de solução tem a seguinte forma
u(t) = φf(t). (119)
Paulo de Tarso R Mendonça - Grante - EMC - UFSC, 2004 31
Consideremos o signiÞcado desta função. A Figura 17 ilustra o caso de uma viga modelada por
três elementos Þnitos e quatro nós. A eq.(119) representa os deslocamentos de cada nó ao longo do
tempo: ¯¯¯¯
¯¯¯¯ 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, por
exemplo, isto signiÞca, por esta hipótese que os todos os deslocamentos nodais no instante, t = t2
são f2/f1 vezes maior que no instante t = t1. Estes valores nodais são sempre proporcionais entre
si 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 no
tempo.
Observe que (119) não é a solução de (118), mas apenas sua forma geral. Antes desconhecíamos
os valores nodais da função do tempo u(t). Agora temos u(t) expresso em termos de um perÞl
de deslocamentos nodais φ, independente do tempo, e de um fator comum, f(t), ambos também
desconhecidos. A diferenca é que antes tínhamos n funções do tempo a determinar, agora as n
incógnitas são constantes, as componentes de φ, e apenas uma função incógnita dependente do
tempo, 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:
Paulo de Tarso R Mendonça - Grante - EMC - UFSC, 2004 32
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 se
obtém uma nova equação, em termos apenasde f(t). Se f¨(t)/f(t) = −λ, de (122) temos os dois
problemas:
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ção
caracterí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 é nulo
e o sistema é dito conservativo, isto é, não possui dissipação de energia, amortecimento. A única
forma deste sistema se mover é simplesmente continuar ummovimento iniciado anteriormente através
de um impulso rápido aplicado no instante inicial. O movimento deve ser tal que a quantidade de
energia total do sistema deve permanecer constante. Isto não permite que a amplitude do movimento
cresç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) = A1e
iω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)
Paulo de Tarso R Mendonça - Grante - EMC - UFSC, 2004 33
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 do
problema (123)2. Uma vez determinados φ e ω de (123), e usando f(t) de (128), a solução do
problema 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 com
freqüê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ção
matricial é um sistema algébrico de n equações e n incógnitas, os φi, i = 1, n. De álgebra linear
sabe-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) para
determiná-las, bastando que procuremos os valores de ω para os quais o determinante de A seja
nulo. 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 chamado
determinante característico. Se tivermos um sistema pequeno, 2 × 2 por exemplo, as raízes do
polinômio podem ser obtidas analiticamente. Para sistemas da ordem de milhares, como é comum em
elementos Þnitos, existem alguns métodos numéricos e dezenas de variações, que estimam algumas
ou 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 um
sistema de n equações reais possui um polinômio característico de grau n e possui n raízes. O
problema (130) é denominado problema de autovalor ou autoproblema, enquanto as raízes do
polinômio característico, os ω0is, são chamados autovalores do problema. Calculamos então os n
autovalores λj = ω2j do problema. A cada autovalor substituído em (130) poderemos resolver e obter
um distinto vetor solução φj, isto é£
K− ω2jM
¤
φj = 0 j = 1, 2, · · · , n. (134)
Paulo de Tarso R Mendonça - Grante - EMC - UFSC, 2004 34
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 par
de 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 modos
naturais 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 segunda
por exemplo, e o sistema Þca
φsTKφr = ω2rφ
sTMφr e φsTKTφr = ω2sφ
sTMTφr. (138)
Como K e M são matrizes simétricas, os termos se tornam idênticos entre as duas equações. Se
subtrairmos 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.
Paulo de Tarso R Mendonça - Grante - EMC - UFSC, 2004 35
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ção
de Gran-Schmidt.
Assim,deforma geral, considera-se sempre que se tem o conjunto
dos 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 umamatriz
peso, no caso é usado a matriz massa, de forma que U · V ≡ UTMV, enquanto a forma mais
conhecida do produto interno é o chamado produto euclidiano, dado por U ·V ≡ UTV. Observe
que se dois distintos vetores modais φr e φs são M-ortogonais, isto é, satisfazem (140), eles são
també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émchamado produto interno entre dois vetores, temos
entã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ção
perde 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 positiva
deÞnida. Uma matriz A é dita positiva deÞnida se°°°° VTAV ≥ 0, para qualquer V, eVTAV = 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). Suponha
que 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 cada
autovetor fazemos sua normalização, isto é, calculamos sua norma por (142) e fazemos
φr =
1
kφrkφ
r. (144)
Paulo de Tarso R Mendonça - Grante - EMC - UFSC, 2004 36
Visivelmente, agora φr tem norma unitária, isto é, φrTMφr = 1. Como os diversos autovetores
sã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 se
r = 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 conjunto
ortonormal de vetores, isto é, cada um é normalizado para norma unitária e também é ortogonal
a todos os demais. A relação (145) está colocada na chamada forma indicial, com os índices r e s
podendo assumir valores entre 1 e n. Essa relação pode também ser colocada numa forma matricial
completa. 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. Em
vez de representar um autopar de solução a cada vez, como em (136)1, aplicamos todos os autopares
simultaneamente. 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 φ
2
1 · · · φn1
φ12 φ
2
2
...
...
. . .
...
φ1n φ
2
n φ
n
n
⎤⎥⎥⎥⎦
=
⎡⎢⎢⎢⎣
M11 M12
M21 M22
. . .
Mnn
⎤⎥⎥⎥⎦
⎡⎢⎢⎢⎣
φ11 φ
2
1 · · · φn1
φ12 φ
2
2
...
...
. . .
...
φ1n φ
2
n φ
n
n
⎤⎥⎥⎥⎦
⎡⎢⎢⎢⎣
ω21
ω22
. . .
ω2n
⎤⎥⎥⎥⎦ .
Paulo de Tarso R Mendonça - Grante - EMC - UFSC, 2004 37
Note que (149) representam as n equações associadas a cada autovetor. Então tem-se de fato
n× n equações algébricas. Podemos em seguida pré-multiplicar (149) por ΦT obtendo:
ΦTKΦ = Φ
t
MΦΛ2.
Mas com a ortonormalidade dos autovetores em relação à matriz massa, eq. (148), o lado direito
simpliÞ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 de
rigidez é o quadrado do correspondente autovalor. A seção ?? descreve também outras propriedades
das 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. Use
E = 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 0
0 −1 2 −1
0 0 −1 1
⎤⎥⎥⎦− αj
⎡⎢⎢⎣
2 1 0 0
1 4 1 0
0 1 4 1
0 0 1 2
⎤⎥⎥⎦
⎫⎪⎪⎬⎪⎪⎭
⎡⎢⎢⎣
φj1
φj2
φj3
φj4
⎤⎥⎥⎦ =
⎡⎢⎢⎣
0
0
0
0
⎤⎥⎥⎦ . (152)
onde αj = ω2j
ρAL
18
.
Deve-se primeiramente aplicar as condições de contorno para vincular a barra. Uma vez que ela
está 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Þca
eliminar a primeira coluna de cada matriz junto com o termo φj1. Em seguida eliminamos a primeira
linha, Þcando então com matrizes 3× 3. Quanto às constantes multiplicativas, dividimos a equação
por 3EA/L e deÞnimos
αj = ω
2
j
ρAL/18
3EA/L
= ω2j
ρL2
54E
(153)
O polinômio característico então Þca
p(αj) = det
⎧⎨⎩
⎡⎣ 2 −1 0−1 2 −1
0 −1 1
⎤⎦− αj
⎡⎣ 4 1 01 4 1
0 1 2
⎤⎦⎫⎬⎭ = 0, (154)
que pode ser simpliÞcado para:
Paulo de Tarso R Mendonça - Grante - EMC - UFSC, 2004 38
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
1
4L
p
E/ρ = 1.265, 9 1.395,8 (10,3 %) 1.298,3 (2,56 %) 1.280,4 (1,10 %)
ω2
3
4L
p
E/ρ = 3.797, 8 4.536,5 (19,5 %) 4.187,6 (10,3 %)
ω3
5
4L
p
E/ρ = 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 nas
aproximaçõ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 modos
naturais 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 −1
0 −1 1
⎤⎦− αj
⎡⎣ 4 1 01 4 1
0 1 2
⎤⎦⎫⎬⎭
⎡⎣ φj2φj3
φj4
⎤⎦ =
⎡⎣ 00
0
⎤⎦ . (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 corresponde
a α1 = 0,04674m−2s−4. A eq.(156) para o modo j = 1 Þca⎡⎣ 1, 81302 −1, 04674 0−1, 04674 1, 81302 −1, 04674
0 −1, 04674 0, 90651
⎤⎦⎡⎣ φ12φ13
φ14
⎤⎦ =
⎡⎣ 00
0
⎤⎦ .
Triangularizando a matriz por fatorização gaussiana temos:
Paulo de Tarso R Mendonça - Grante - EMC - UFSC, 2004 39
⎡⎣ 1, 81302 −1, 04674 00 1, 20868 −1, 04674
0 0 0
⎤⎦⎡⎣ φ12φ13
φ14
⎤⎦ =
⎡⎣ 00
0
⎤⎦ .
Podemos fazer φ14 = 1, 0. Neste caso resolvemos φ
1
2 = 1, 5 e φ
1
3 = 1, 15471, isto é,
φ1 =
⎡⎣ 0, 50, 866
1, 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, 866
1, 0
⎤⎦ .
8.2.5 Exemplo 5 — Normalização de Autovetores
Considere o autoproblema do Exemplo 4. Normalize os autovetores, forme as matrizes massa e
rigidez.