Captulo 1

Introduo s Vibraes
Mecnicas com
MATLAB/Simulink: Simulao e
Animao
Diego A. Duarte

Neste captulo so apresentados conceitos bsicos para o estudo das vi-

braes mecnicas de partculas com o programa MATLAB/Simulink. Nas
primeiras sees introduzido o mtodo de diagrama de blocos, sua utilizao
como ferramenta para resolver problemas de engenharia e os procedimentos
necessrios para implementao destes diagramas no Simulink. Em seguida,
so apresentados alguns cdigos bsicos para animar o movimento de um
sistema massa-mola com vibrao forada. No nal do captulo, exemplos
mais complexos de vibrao mecnica so brevemente discutidos.

1.1 Anlise de sinais e sistemas com diagramas

de blocos

O diagrama de blocos uma estrutura lgica que representa uma relao

matemtica entre dados de entrada e sada de um determinado sistema fsico.
Quando um locutor fala em um microfone, por exemplo, a onda sonora emi-
tida pela sua voz transmitida para uma central de processamento que faz
a onda chegar at o ouvinte. Neste caso, a onda sonora emitida representa
o sinal de entrada e a onda captada o sinal de sada. Para que a voz tenha

1
2 CAPTULO 1. SIMULAO NUMRICA COM MATLAB

sonoridade, necessrio, na maioria das vezes, amplicar o sinal de entrada,

i.e., aumentar a amplitude da onda sonora. Esta etapa realizada por um
sistema que denido como um processo capaz de estabelecer uma relao
entre sinais [1]. Uma das maneiras mais bsicas de estabelecer uma relao
matemtica entre dois ou mais sinais por meio de equaes diferenciais.
Nas mais diversas reas da Cincia, vrios problemas podem ser mode-
lados e estudados com equaes diferenciais. Uma equao diferencial, con-
forme o prprio nome diz, uma equao formada por derivadas de alguma
funo. A equao apresentada abaixo um exemplo de equao diferencial.
Ela representa um circuito RLC em uma nica malha:

q + Rq + C 1 q = V (t)
L (1.1)

A constante L representa a indutncia do circuito, R a resistncia, C

a capacitncia, q a carga eltrica e V (t) a tenso eltrica de alguma fonte
externa. A soluo desta equao fornece o comportamento temporal da
carga eltrica armazenada no capacitor. Desta maneira, podemos compreen-
der como realizada a carga ou a descarga do capacitor e como a corrente
eltrica se comporta no circuito. O comportamento da tenso V (t) pode ser
na forma contnua, alternada, pulsada etc. Neste caso, V (t) representa o
sinal de entrada e q(t) o sinal de sada. O indutor, resistor e o capacitor
representam o sistema que vai processar o sinal V (t).
Uma equao diferencial pode ser resolvida por mtodos analticos ou
numricos. As solues analticas, conforme voc estudou em disciplinas da
clculo, envolvem diversos mtodos de resoluo. O problema que nem
sempre podemos resolver uma equao diferencial por mtodos analticos,
sendo necessria a utilizao de mtodos numricos. Isto acontece quando
trabalhamos com equaes diferenciais mais elaboradas e difceis de serem
resolvidas analiticamente. As equaes diferenciais que descrevem circuitos
eltricos mais complexos, por exemplo, no possuem solues analticas e
devem ser resolvidas numericamente.
Muitos problemas de fsica ou engenharia exigem um esforo aritmtico
considervel quando resolvido manualmente. Desta forma, aps compreender
os algoritmos de resoluo, natural utilizarmos computadores para acelerar
e renar a preciso destas solues. Para isso, precisamos conhecer lingua-
gens de programao.
Atualmente, existem diversas lnguas de programao e, basicamente,
no existe a melhor linguagem. Quem dene isso o programador, pois esta
concluso depende muito da sua empatia pela estrutura dos cdigos, layout
dos compiladores, objetivos do seu trabalho etc. Um dos critrios que pode
inuenciar na sua deciso a biblioteca de cdigos pr-denidos inseridos
no seu compilador, pois, dependendo do modelo que ser estudado, alguns
numericamente.

Muitos problemas de fsica ou engenharia, resolvidos numericamente, exigem

um esforo aritmtico considervel quando resolvido manualmente. Desta forma,
aps compreender os algoritmos de resoluo, natural utilizarmos computadores
para acelerar e refinar a preciso destas solues. Para isso, precisamos conhecer
linguagens de programao.
1.1.
existem diversas lnguas de programao e, basicamente, 3no
ANLISE DE SINAIS COM DIAGRAMA DE BLOCOS
Atualmente,
existe a melhor linguagem. Quem define isso o programador, pois esta concluso
depende muitoou
programas dalinguagens
sua empatia pela estrutura
possuem cdigosdos cdigos, layout
pr-denidos mais dos compiladores,
fceis de serem
objetivos do seu trabalho
implementados. Paraetc. Um dos
resolver umacritrios quediferencial
equao pode influenciar na sua deciso
pelo mtodo de Eu- a
biblioteca
ler, porde cdigos voc
exemplo, pr-definidos inseridos
pode escrever no seu compilador,
o mtodo pois, dependendo
numrico explcito ou apenasdo
modelo quea ser
digitar estudado,
equao algunse chamar
diferencial programas umaoufuno
linguagens possuem
que far cdigospor
o trabalho pr-
definidos mais fceis de serem implementados. Para resolver uma equao
voc. Na segunda opo, o programa aplicar a equao diferencial em um diferencial
pelocdigo
mtodo de Euler,e, por
pr-denido exemplo,apresentar
em seguida, voc pode diretamente
escrever o asmtodo numrico
solues. Se
voc j um expert no mtodo de Euler, no faz sentido escrever vrias linhas
explicitamente ou apenas digitar a equao diferencial e chamar uma funo que far
o trabalho
para implementao do mtodo. Basta chamar uma funo pr-denida. Po-em
por voc. Na segunda opo, o programa aplicar a equao diferencial
umrm,
cdigo pr-definido e, em seguida, apresentar diretamente as solues. Se voc j
se voc ainda aprendiz, fundamental digitar o cdigo explcito para
um expert no mtodo de Euler, no faz o mnimo sentido escrever linhas e linhas de
aprender programao e aguar o raciocnio lgico.
comandos. Basta chamar uma funo pr-definida. Porm, se voc ainda aprendiz,
Alm dos mtodos tradicionais de programao, existem ferramentas mais
fundamental digitar o cdigo explicitamente para aprender programao e aguar o
agradveis aos olhos do programador e que no necessitam, em princpio, de
raciocnio lgico.
um conhecimento profundo em linguagens de programao. Um destes mto-
dos a programao com diagrama de blocos. Esta tcnica consiste na repre-
Alm dos mtodos tradicionais de programao, existem mtodos visualmente
sentao grca de um determinado modelo fsico (e.g.: sistema massa-mola)
mais agradveis aos olhos do programador e que no necessitam, em princpio, de um
com diversos blocos em que cada um representa uma operao matemtica.
conhecimento profundo em linguagens de programao. Um destes mtodos a
Em um diagrama de blocos, a equao 1.1 pode ser representada pela es-
programao com diagrama de blocos. Esta tcnica consiste na representao grfica
trutura
de um da gura
determinado 1.1. Neste
modelo fsicocaso,
(e.g.: osistema
sistema (circuito RLC)
massa-mola) composto
com diversos porem
blocos
subsistemas que incluem todas as operaes matemticas necessrias
que cada um representa uma operao matemtica. Em um diagrama de blocos, a para
converter
equao (1) orepresentada
sinal de entrada V (t) noestrutura:
pela seguinte sinal de sada q(t).

ENTRADA V(t) SISTEMA SADA q(t)

FIGURA 1. REPRESENTAO ESQUEMTICA DO SINAL V(t) SENDO PROCESSADO POR UM CIRCUITO RLC.
Figura 1.1: Representao esquemtica de um diagrama de blocos.
em que o sistema (circuito RLC) composto por subsistemas que incluem todas as
operaes matemticas necessrias para converter o sinal de entrada V(t) no sinal de
sada q(t).
1.1.1 Representao de equaes diferenciais em dia-
grama de
1. Representao deequaes
blocosdiferenciais em diagrama de blocos

Para exemplicar a construo de um diagrama de blocos, vamos utilizar

Para exemplificar a construo de um diagrama de blocos, vamos utilizar uma
uma equao diferencial ordinria (EDO ou ODE do ingls ordinary die-
equao diferencial ordinria (EDO ou ODE do ingls ordinary differential equation)
rential equation ) com coecientes constantes, pois, estruturalmente, uma
com coeficientes constantes, pois, estruturalmente, uma das equaes diferenciais
das equaes diferenciais mais simples de estudar.
mais simples de estudar. Como mostra a equao (1), uma equao diferencial uma
Como mostra a equao 1.1, uma equao diferencial uma equao
equao composta por derivadas de alguma funo. Quando a derivada aplicada em
composta por derivadas de alguma funo. Quando a derivada aplicada em
relao uma nica varivel independente ( e.g.: tempo), esta equao recebe
o nome de equao diferencial ordinria (a palavra ordinria, na matemtica,
signica que h uma nica varivel independente na equao). Observe que
as suas derivadas so multiplicadas por constantes. Se estes valores perma-
necem xos, em toda faixa de tempo analisada, a equao chamada de
4 CAPTULO 1. SIMULAO NUMRICA COM MATLAB

equao diferencial ordinria com coecientes constantes. Se o sinal de en-

trada (tenso V (t)) zero, temos uma EDO homognea. Caso contrrio,
temos uma EDO no homognea.
A equao pode apresentar derivadas de vrias ordens. A derivada de
maior ordem dene a ordem da equao diferencial. A equao 1.1 uma
EDO de segunda ordem. Alm destas classicaes, uma equao diferencial
pode ser classicada como linear e no linear. Se a equao possui uma de-
rivada com potncia maior que 1 ( : q e.g.
2 ), chamada de no linear. Caso a
potncia da derivada seja 0 ou 1, a equao classicada como linear. Desta
forma, a equao 1.1 uma EDO de segunda ordem, linear, no homognea
e com coecientes constantes. Uma equao classicada desta forma pode
ser representada como:

n m
X di y X dk x
i = k (1.2)
i=0 dti k=0 dtk
em que i e k x(t) o sinal de entrada e y(t) o
so coecientes constantes,
sinal de sada. Os ndices n e m medem a ordem da EDO. Para compreender
como montar um diagrama de blocos com a equao 1.2, vamos particulariz-
la para uma EDO de primeira ordem e, em seguida, para um EDO de segunda
ordem. Assim:
1 y + 0 y = 1 x + 0 x (1.3)

Para resolver a equao 1.3 por mtodos numricos, podemos aplicar

qualquer uma das tcnicas que conhecemos. Porm, vamos transformar esta
equao em um diagrama de blocos e implement-lo no programa Simulink.
Este programa um dos recursos do MATLAB e nele basta chamar alguma
funo para resolver a equao diferencial por mtodos numricos. Entre as
tcnicas disponveis na biblioteca do MATLAB, h o mtodo de Euler (funo
ode1) e o mtodo de Runge-Kutta de quarta ordem (funo ode4). Mtodos
mais renados como Dormand-Prince (funo ode45) e Bogacki-Shampine
(funo ode23) tambm esto disponveis.
Cada termo da equao 1.3 representa uma parte do diagrama. Isolando
a derivada de maior grau do sinal de sada:

y = 1 (1 x + 0 x 0 y)
e integrando ambos os lados de 0 at t,
Z t Z t
0 )dt0 = y(t) y(0) = 1 1
y(t 0 ) + 0 x(t0 ) 0 y(t0 )]dt0
[1 x(t
0 0

obtemos a soluo:
Z t
1
y(t) = y(0) + 1 0 ) + 0 x(t0 ) 0 y(t0 )]dt0
[1 x(t
0
1.1. ANLISE DE SINAIS COM DIAGRAMA DE BLOCOS 5

que pode ser escrita como:

Z t Z t
1 0 0
y(t) = y(0) + 1 [1 x(t) + 0 x(t )dt 0 y(t0 )dt0 ] (1.4)
| {z } 0 0
| {z } | {z }
C1 C2 C3
que pode ser escrita como:
Para montar o diagrama, vamos utilizar a equao 1.4. Para isso, precisa-
remos de dois blocos
1 integradores,
quatro blocos amplicadores e dois blocos
) )
() = (0) + [

somadores. Acompanhe
1 1 () +
0 (
0 ( ]
pela0 gura 1.2. O 0bloco amplicador representar
(4)
1 2 as operaes 3 de soma e subtrao. Aps
as constantes e o bloco somador
o sinal x(t) entrar no sistema, haver uma bifurcao (representado por um
Para montar o diagrama, vamos utilizar a equao (4). Para isso, precisaremos
ponto
de doisescuro). Neste ponto,
blocos integradores, o blocos
quatro sinal passar por duas
amplificadores e doisvias.
blocosEm cada via
somadores.
oAcompanhe
sinal possui
pelaa figura
mesma intensidade
2. O do sinal
bloco amplificador original. as
representar Aps passarem
constantes pelos
e o bloco
somador as operaes
amplicadores 1 e 2,de ambos
soma e subtrao. Aps o sinal
sero somados x(t) entrar
no bloco no sistema,
somador haver
1. Observe
umano
que bifurcao (representado
bloco integrador por um
1 deve ponto
ser escuro). Neste
informada ponto, oinicial
a condio sinal passar porde
do sinal
duas vias. Em cada tvia
entrada no tempo
o sinal possui a mesma intensidade do sinal original. Aps
0 . O resultado do bloco somador 1 ser a soma das com-
passarem pelos amplificadores 1 e 2, ambos sero somados no bloco somador 1.
ponentes C1 e C2 da equao 1.4. Para somar a componente C3 adicionamos
Observe que no bloco integrador 1 deve ser informada a condio inicial do sinal de
mais umnobloco
entrada temposomador e umdo
t0. O resultado bloco
blocoamplicador
somador 1 seraps o bloco
a soma somador 1.
das componentes C1 O
1
bloco
e C2 amplicador
da equao (4).3 representar o termo 1 C3,da
Para somar a componente equao 1.4.
adicionamos maisEm
um seguida
bloco
somador e umoutra
adicionamos bloco bifurcao.
amplificador Uma
aps oviabloco somador
ser o sinal1.de
O bloco
sada amplificador 3
e a outra ser
representar o termo
o clculo da componente1/ 1 da equao (4). Em seguida adicionamos outra bifurcao.
C3 que possui um sinal negativo e que deve ser
Uma via ser o sinal de sada e a outra ser o clculo da componente C que possui
adicionado no diagrama aps o bloco amplicador 4. O termo3 y0 adicio-
um sinal negativo e que deve ser adicionado no diagrama aps o bloco amplificador 4.
nado no integrador 2. A estrutura representada pela gura 1.2 chamada de
O termo y0 adicionado no integrador 2. A estrutura representada pela figura 2
diagrama de blocos na forma direta 1 [1].
chamada de diagrama de blocos na forma direta 1 [1].

Bifurcao Bifurcao

Amplificador 1 Somador 1 Somador 2 Amplificador 3

Entrada Sada
x(t) 1 1-1 y(t)

Integrador 1 - Integrador 2

x(0) y(0)

0 0
Amplificador 2 Amplificador 4
FIGURA 2. DIAGRAMA DE BLOCOS NA FORMA DIRETA 1 PARA UMA EDO DE 1 ORDEM.

Figura 1.2: Diagrama de blocos na forma direta 1 para uma EDO de primeira
Similarmente, podemos representar um diagrama de blocos na forma direta 1
ordem.
para uma EDO de qualquer ordem. Para uma EDO de segunda ordem:

2 Similarmente,
+ 1 + 0 =podemos
2 + 1 representar
+ 0 um diagrama de blocos na forma di-
(5)
reta 1 para uma EDO de qualquer ordem. Para uma EDO de segunda ordem,

o diagrama representado na figura 3. Para obt-lo basta isolar a derivada de maior

grau do sinal de sada eseguir
2y + os1 ymesmos
+ 0 y procedimentos.
= 2 x + 1 x + 0 x (1.5)

A equao (5) a representao geral da equao (1). Assim, com as devidas

atribuies para as constantes do diagrama da figura 3, podemos escrever o diagrama
de blocos para o circuito RLC, conforme representado pela figura 4. Obviamente, os
amplificadores com ganho nulo podem ser retirados do diagrama, pois no haver
6 CAPTULO 1. SIMULAO NUMRICA COM MATLAB

o diagrama representado na gura 1.3. Para obt-lo basta isolar a derivada

de maior grau do sinal de sada e seguir os mesmos procedimentos.
A equao 1.5 a representao geral da equao 1.1. Assim, com as
devidas atribuies para as constantes do diagrama da gura 1.3, podemos
escrever o diagrama de blocos para o circuito RLC, conforme representado
pela gura 1.4. Obviamente, os amplicadores com ganho nulo podem ser
retirados do diagrama, pois no haver uxo de sinal por eles. Similarmente,
o amplicador com ganho unitrio tambm no precisa ser adicionado. O
fluxo de sinal por eles. Similarmente, o amplificador com ganho unitrio tambm no
diagrama simplicado apresentado na gura 1.5.
precisa ser adicionado. O diagrama simplificado apresentado na figura 5.
fluxo de sinalBifurcao
por eles. Similarmente, o amplificador com ganho unitrio Bifurcao tambm no
precisa ser adicionado. O diagrama simplificado
Amplificador 1 Somador 1 Somador 3 Amplificador 4 figura 5.
apresentado na
Entrada Sada
x(t)
Bifurcao 2 2-1 Bifurcao
y(t)

Entrada Amplificador
Integrador11 Somador 1
-
Somador 3 Amplificador
Integrador 3 4 Sada

x(t) 2 x(0)
Somador 2 y(0)-1
2 y(t)

Bifurcao
Integrador 11 - 1
Integrador 3
Bifurcao

Amplificador 2
x(0) Somador 4 Amplificador
y(0) 5
Somador 2
x(0) y(0)
Bifurcao 1
Integrador 2 1 4
Integrador Bifurcao

2
Amplificador
0
Somador 4 Amplificador
0
5

Amplificador 3
x(0) Amplificador
y(0) 6
Integrador
FIGURA23. DIAGRAMA DE BLOCOS NA FORMA DIRETA 1 PARA UMA EDO Integrador
DE 2 ORDEM. 4

Figura 1.3: Diagrama 0 de blocos na forma direta 1 para uma EDO de segunda
V(t) 0 L-10 q(t)
ordem. Amplificador 3 Amplificador 6 -
FIGURA 3. DIAGRAMA DE BLOCOS NA FORMA DIRETA 1 PARA UMA EDO DE 2 ORDEM.
V(0) q(0)

V(t) 0 0 RL-1 q(t)

-
V(0) q(0)
V(0) q(0)

0 1 C-1R
FIGURA 4. DIAGRAMA DE BLOCOS NA FORMA DIRETA 1 PARA UM CIRCUITO RLC.

V(0) q(0)
V(t) L-1 q(t)
-
1
V(0) V(0) q(0)
C-1
FIGURA 4. DIAGRAMA DE BLOCOS NA FORMA DIRETA 1 PARA UM CIRCUITO RLC.
R
Figura 1.4: Diagrama de blocos na forma direta 1 para um circuito RLC.
V(t) L-1 q(t)
q(0)
-
V(0) V(0) Cq(0)
-1

FIGURA 5. DIAGRAMA DE BLOCOS NA FORMA DIRETA 1 (VERSO SIMPLIFICADA) PARA UMA CIRCUITO RLC.
R

q(0)

-1
 V(0) q(0)

0 R

V(0) q(0)

1 C-1
1.2. RESOLUO DE DIAGRAMA DE BLOCOS COM SIMULINK 7
FIGURA 4. DIAGRAMA DE BLOCOS NA FORMA DIRETA 1 PARA UM CIRCUITO RLC.

V(t) L-1 q(t)

-
V(0) V(0) q(0)

q(0)

C-1
FIGURA 5. DIAGRAMA DE BLOCOS NA FORMA DIRETA 1 (VERSO SIMPLIFICADA) PARA UMA CIRCUITO RLC.

Figura 1.5: Diagrama de blocos na forma direta 1 para um circuito RLC

(verso simplicada).

A forma direta 1 um mtodo que utiliza 2N integradores (N a ordem

da EDO). H mtodos que reduzem a quantidade de integradores para a
quantidade mnima necessria para minimizar a utilizao de memria do
computador. Estas estruturas so chamadas de formas cannicas. Entre
os principais, esto a forma direta 2 e o modelo de espao de estados. No
entanto, os clculos que sero realizados no demandam um considervel
poder computacional. Assim, vamos utilizar apenas a forma direta 1. Se o
leitor tiver interesse por outros mtodos, veja a ref. [1].

1.2 Resoluo de diagrama de blocos com Si-

mulink

Na seo anterior vimos como representar uma equao diferencial em

diagramas de blocos. Nesta seo vamos resolver a equao diferencial e
encontrar o sinal da sada a partir de um determinado sinal de entrada. Para
isso, suponha a EDO de primeira ordem:

3y + 2y = 7x + 4x (1.6)

considerando x(0) = 0 e y(0) = 0. O diagrama na forma direta 1 apresen-

tado na gura 1.6.
Para resolver este diagrama, vamos utilizar o programa Simulink. Na
janela de comando do MATLAB, digite simulink para abrir a biblioteca de
blocos do programa. Em seguida, na biblioteca de blocos, clique em  File ,
 New  e  Model . Aparecer uma nova janela com fundo branco. neste
 2. Resoluo de diagrama de blocos com Simulink
Na seo anterior vimos como representar uma equao diferencial em
diagramas de blocos. Nesta seo vamos resolver a equao diferencial e encontrar o
sinal da sada a partir de um determinado sinal de entrada. Para isso, suponha a EDO
de primeira ordem:

3 + 2 = 7 + 4 (1)
8 CAPTULO 1. SIMULAO NUMRICA COM MATLAB

considerando x(0) = 0. O diagrama na forma direta 1 apresentado na figura 6.

x(t) -7 1/3 y(t)

-
x(0) y(0)

4 2
FIGURA 6. DIAGRAMA DE BLOCOS NA FORMA DIRETA 1 PARA A EQUAO (1).

Figura 1.6: Diagrama de blocos na forma direta 1 para a equao 1.6.

Para resolver este diagrama, vamos utilizar o programa Simulink. Na janela de
comando do MATLAB, digite simulink para abrir a biblioteca de blocos do
programa.
arquivo queEmmontaremos
seguida, na biblioteca de blocos,
o diagrama com os clique em File,
blocos. New e de
Na biblioteca Model.
blocos,
h uma lista de opes no lado esquerdo. Clique na pasta  Continuous . Em
Aparecer uma nova janela com fundo branco. neste arquivo que montaremos o
seguida, clique uma vez no bloco  Integrator  e arraste-o para a janela com
diagrama com os blocos. Na biblioteca de blocos, h uma lista de opes no lado
esquerdo. Clique na pasta Continuous. Em seguida, clique uma vez no bloco
fundo branco. A gura 1.6 mostra que precisamos de dois blocos integrado-
Integrator e arraste-o para a janela com fundo branco. A figura 6 mostra que
res, assim, repita
precisamos este procedimento.
de dois blocos integradores, assim, repita este procedimento.
Na lista de opes da biblioteca, clique na pasta  Math Operations . Pe-
gue doisNa lista de
blocos opes da ebiblioteca,
somadores clique na pasta Math
quatro amplicadores. EstesOperations. Pegue
blocos so chama-
doisde
dos Sum
blocos  e  Gain
somadores e quatro amplificadores. Estes
, respectivamente. blocos
Para so chamados
atribuir valores de Sum
numricos
e Gain, respectivamente. Para atribuir valores numricos para os
para os amplicadores, clique duas vezes sobre o bloco. O segundo somador
amplificadores,
clique duas vezes sobre o bloco. O segundo somador dever ter um sinal negativo.
dever ter um sinal negativo. Para isso, clique nele duas vezes. Uma janela
Para isso, clique nele duas vezes. Uma janela de propriedades abrir. Clique em List
de
of propriedades
signs e troque oabrir.
smbolo Clique |+.  List of signs  e troque o smbolo |++
|++ por em
por |+-.
ParaPara definir
denir o sinaldedeentrada,
o sinal entrada, clique
cliquenana
pasta Sources.
pasta  SourcesNela h diversas
. Nela, h di-
formas formas
versas para o sinal
parade oentrada.
sinal Vamos selecionar,
de entrada. por exemplo,
Vamos o bloco
selecionar, porSine Wave. o
exemplo,
bloco  Sine devemos
Neste caso, Wave . informar a amplitude, polarizao, frequncia angular, fase e o
Neste caso, devemos informar a amplitude, polarizao,
tempo de amostragem. Os valores padres so definidos como 1, 0, 1, 0 e 0,
frequncia angular, fase e o tempo de amostragem. Os valores padres so
respectivamente. Para mudar estes valores, basta clicar duas vezes sobre o bloco. A
denidos
polarizaocomo 1, 0, 1, 0 evertical
o deslocamento 0, respectivamente.
do sinal de entradaPara mudar
e o tempo estes valores,
de amostragem
basta clicar duas
a discretizao vezes
do sinal. sobre
Para o bloco.
visualizar de sada y(t), pegue
A polarizao
o sinal o bloco Scopeverti-
o deslocamento na
pasta
cal doSinks.
sinal de entrada e o tempo de amostragem a discretizao do sinal.
Para visualizar o sinal de sada y(t), pegue o bloco  Scope  na pasta  Sinks .
Com estes blocos podemos montar o diagrama da figura 6, conforme mostra a
Com estes blocos podemos montar o diagrama da gura 1.6, conforme
figura 7. Durante a montagem, ser necessrio rotacionar alguns blocos. Para isso,
mostra a gura 1.7. Durante a montagem, ser necessrio rotacionar alguns
clique uma vez com o boto direito, selecione Format, Rotate block e o sentido
blocos. Para isso, clique uma vez com o boto direito, selecione  Format ,
de rotao (horrio ou anti-horrio).
 Rotate block  e o sentido de rotao (horrio ou anti-horrio).
Para resolver a equao, devemos selecionar o mtodo numrico de reso-
luo. Para isso, clique em  Simulation , na janela de edio do diagrama
no Simulink, e, em seguida,  Conguration Parameters  ou  Model Congu-
ration Parameters . Abrir uma janela com as conguraes de simulao.
Nela, clique em  Solver  na lista do lado esquerdo. Em seguida, clique no
campo  Solver  e escolha o mtodo desejado. Como a equao diferencial
deste exemplo de primeira ordem, podemos escolher o mtodo de Euler.
1.2. RESOLUO DE DIAGRAMA DE BLOCOS COM SIMULINK 9

-7 1/3

Sine Wave Scope

Gain3 Gain

1 1
Integrator1 Integrator
s s
Gain1

4 2

Gain2

Figura 1.7: Diagrama de blocos na forma direta 1 para a equao 1.1 (verso
Simulink).

Para isso, clique no campo  Type  Fixed-step . Em seguida,

e selecione 
volte para o campo  Solver  e selecione  ode1 (Euler) . Logo abaixo, na
mesma tela, no campo  Fixed-step size (fundamental sample time) , fornea
um valor para o passo de integrao (e.g.: 0.1) e clique em  Apply . Para se-
lecionar o intervalo de integrao, clique no campo  Start time , para digitar
o tempo inicial, e  Stop time  para digitar o tempo nal. Os valores padres
so 0.0 e 10 s, respectivamente. Veja a gura 1.8.
Para inserir a condio inicial da equao 1.6, clique duas vezes sobre o
bloco integrador e digite zero no campo  Initial condition . Para resolver o
diagrama, clique em  Simulation  na janela de edio do diagrama no Simu-
Start  ou  Run . Para visualizar o sinal de sada, clique
link, e, em seguida, 
duas vezes no bloco  Scope . O resultado est representado na gura 1.9.
Alm do sinal de sada, podemos colocar mais um canal no osciloscpio
e visualizar o sinal de entrada. Para isso, clique duas vezes no bloco  Scope 
e, em seguida, no cone  Parameters . Nas opes deste cone, digite 2 no
campo  Number of axis  e clique em  Apply . Volte para o diagrama e faa
a conexo de um canal com o sinal de sada e o outro com o sinal de entrada,
conforme mostra a gura 1.10. Em seguida, faa uma nova simulao e veja
o resultado no osciloscpio. O resultado representado na gura 1.11.
Alm de obter os grcos diretamente pelo Simulink, podemos exportar
os dados para a tela de comando do MATLAB. Para isso, na biblioteca de
blocos, clique em  Sink , arraste o bloco  To workspace  e conecte-o no sinal
de sada. Clique duas vezes no bloco e coloque o nome da varivel e logo
abaixo, no campo  Save Format , selecione  Array . Para exportar o sinal
de entrada, repita o procedimento e, na tela de edio de blocos, clique em
 Start . O diagrama nesta forma est representado na gura 1.12.
Depois de exportados, os dados podem ser acessados na tela de comando
do MATLAB. Basta digitar o nome da varivel de interesse. Para represen-
10 CAPTULO 1. SIMULAO NUMRICA COM MATLAB

Figura 1.8: Conguraes do mtodo de resoluo do diagrama de blocos.

Figura 1.9: Sinal de sada gerado pelo diagrama da gura 1.7.

1.2. RESOLUO DE DIAGRAMA DE BLOCOS COM SIMULINK 11

-7 1/3
Scope
Sine Wave
Gain3 Gain

1 1
Integrator1 Integrator
s s
Gain1

4 2

Gain2

Figura 1.10: Diagrama de blocos na forma direta 1 com dois canais no osci-
loscpio (verso Simulink).

Figura 1.11: Sinais de entrada e sada para o diagrama da gura 1.10. O

grco superior representa o sinal de entrada e o grco inferior o sinal de
sada.

tar os dados na forma de grco, utilize os cdigos especcos do MATLAB

(comando plot). Para obter o grco da gura 1.9, por exemplo, digite
plot(tout,y). A varivel tout representa o intervalo de integrao. Para
exportar os dados para um arquivo de texto, digite, por exemplo, o cdigo
abaixo na janela de comando do MATLAB:
12 CAPTULO 1. SIMULAO NUMRICA COM MATLAB

To Workspace1

-7 1/3

Sine Wave Scope

Gain3 Gain

1 1 y
Integrator1 Integrator
s s
To Workspace
Gain1

4 2

Gain2

Figura 1.12: Diagrama de blocos na forma direta 1 com a exportao dos

dados para a janela de comando do MATLAB (verso Simulink).

arquivo = fopen('dados.txt','wt');
for i=1:length(y)
fprintf(arquivo,'%f %f %f\n',tout(i),x(i),y(i));
end
fclose(arquivo);

ou digite em um arquivo .m. Para isso, na tela de comando do MATLAB, cli-

que em  File ,  New  e  Blank M-File . Antes de digitar o cdigo, certique
que o arquivo est salvo.

1.3 Vibrao de uma partcula

Nesta seo aplicaremos os conceitos das sees anteriores para estudar

a vibrao livre e forada de uma partcula. Para isso, considere um corpo
de massa m sobre uma superfcie sem atrito e preso na extremidade de uma
mola com constante elstica k. A outra extremidade da mola est presa
em uma parede rgida, conforme mostra a gura 1.13. Pela segunda lei de
Newton, a equao do movimento ser:

d2 x
m = F (t) kx (1.7)
dt2
em que F (t) uma fora externa responsvel pela vibrao do sistema. O
termo kx representa a fora elstica. Ser considerado que o corpo est
inicialmente em repouso e na origem (posio de relaxamento da mola). A
equao 1.7 representa uma EDO de segunda ordem. Para obter a soluo,
sobre umaNesta seo aplicaremos
superfcie os conceitos
sem atrito e preso das sees
na extremidade anteriores
de uma para constante
mola com estudar a
vibraok. livre
elstica e forada
A outra de umadapartcula.
extremidade Para
mola est isso,emconsidere
presa um corpo
uma parede rgida,deconforme
massa m
sobre uma superfcie sem atrito e preso na extremidade de uma mola
mostra a figura 13. Pela segunda lei de Newton, a equao do movimento ser: com constante
elstica k. A outra extremidade da mola est presa em uma parede rgida, conforme
mostra
2 a figura 13. Pela segunda lei de Newton, a equao do movimento ser:
2 = () (1)
2
2 = () (1)
que F(t) uma fora externa responsvel pela vibrao do sistema. O termo kx
em
1.3. VIBRAO DE UMA PARTCULA 13
representa a fora elstica. Ser considerado que o corpo est inicialmente em repouso
eem que F(t)(posio
na origem uma fora externa responsvel
de relaxamento da mola). pela vibrao
A equao (1)do sistema. uma
representa EDOde
O termo kx
representa
vamos
segunda a fora elstica.
represent-la
ordem. Para obteremSer considerado
um diagrama
a soluo, que o corpo est
de blocos, conforme
vamos represent-la inicialmente em
mostra de
em um diagrama repouso
a gura
blocos,
e na origem (posio de relaxamento
1.14. mostra a figura 14.
conforme da mola). A equao (1) representa uma EDO de
segunda ordem. Para obter a soluo, vamos represent-la em um diagrama de blocos,
conforme mostra a figura 14.
F(t)
Parede k F(t)
fixa
Parede k m
fixa
m
FIGURA 13. SISTEMA MASSA-MOLA.

Figura 1.13: Sistema massa-mola.

FIGURA 13. SISTEMA MASSA-MOLA.

F(t) m-1 x(t)

F(t) - m-1 x(t)

F(0) F(0) - x(0)
F(0) F(0) x(0)

x(0)
x(0)
-1
k
k-1
FIGURA 14. DIAGRAMA DE BLOCOS NA FORMA DIRETA 1 PARA O SISTEMA MASSA-MOLA.

Para simular FIGURA

o 14. DIAGRAMA DE BLOCOS
movimento, NA FORMA
sero DIRETA 1 PARA Oos
considerados SISTEMA MASSA-MOLA.
valores arbitrrios: m = 1,0
kg eFigura
k = 25 1.14:
N.m. ODiagrama
sistema ser de analisado
blocos nade duas direta
forma formas:1 (i)paraem ovibrao
sistemalivre
massa-e (ii)
Para forada.
em vibrao
mola.
simular A o movimento,
vibrao livre sero considerados
representada pelaosaplicao
valores arbitrrios:
de uma fora m =que1,0
kg e ksobre
atuar = 25 N.m.
o corpo O sistema
apenas ser analisado
no incio de duas formas:
do movimento (i) emAvibrao
(impulso). vibraolivre e (ii)
forada
em vibrao forada.
Para simular
ser representada A
por ouma vibrao
movimento, livre
fora que sero representada
considerados
acompanhar pela
o corpo osemaplicao
valores de uma
movimento. m
todo arbitrrios: fora que
=
Como
aatuar
1,0 sobre
equao e k
kg (1) ocorpo
=uma 25 apenas
N.m. de
EDO nosegunda
O incio doordem,
sistema movimento
ser vamos(impulso).
analisado de duasoAformas:
resolver vibrao(i)com
diagrama forada
em o
ser representada
vibrao
mtodo por uma
livre e (ii)deem
de Runge-Kutta fora
a que
vibrao
4 ordem acompanhar
forada.
(ode4) com passo o corpo
A vibrao em todo movimento.
livre representada
fixo de 0,001. Como
a equao (1) uma EDO de segunda ordem, vamos
pela aplicao de uma fora que atuar sobre o corpo apenas no incio do
resolver o diagrama com o
a
mtodo de Runge-Kutta de 4 ordem (ode4) com passo
movimento (impulso). A vibrao forada ser representada por uma fora
fixo de 0,001.
que acompanhar o corpo em todo movimento. Como a equao 1.7 uma
EDO de segunda ordem, vamos resolver o diagrama com o mtodo de Runge-
Kutta de quarta ordem (ode4) com passo xo de 0,001.

1.3.1 Vibrao livre

A gura 1.15 mostra o diagrama de blocos no Simulink. Para representar

o impulso de uma fora, usamos o bloco  Step  que est na pasta  Sources 
14 CAPTULO 1. SIMULAO NUMRICA COM MATLAB

da biblioteca de blocos. Conforme mostra a gura 1.16(a), este bloco aplica

uma funo de Heaviside no sinal (funo degrau). Para acessar as congu-
raes da funo, basta clicar duas vezes no bloco.

Step F

To Workspace2

1 1
10 1
s s
Scope
Gain2 Integrator2 Integrator3 Gain1
1
Integrator x
Step1 s
To Workspace

1
Integrator1
s
Gain

25

Figura 1.15: Diagrama de blocos na forma direta 1 para o sistema massa-

mola com impulso unitrio (verso Simulink).

Para aplicar o impulso sobre o corpo, vamos utilizar dois blocos da fun-
o degrau e um bloco somador. Congure os dois primeiros blocos para
iniciarem o degrau em t e t + t, respectivamente. Em seguida, conecte-os
no bloco somador e congure para negativo a entrada do degrau que inicia
em t + t (gura 1.15). Desta forma, o bloco somador ir subtrair o pri-
meiro degrau (gura 1.16(a)) do segundo (gura 1.16(b)) e o sinal resultante
ser um pulso retangular com largura t (gura 1.16(c)). Para aumentar a
intensidade do sinal resultante, coloque um bloco amplicador aps o bloco
somador. No exemplo da gura 1.15, a amplicao 10 N e t = 0,1 s com
t= 1,0 s. Faa tambm a conexo de um bloco  To Workspace  no sinal de
entrada e outro no sinal de sada. Usaremos estes dados para programar a
animao do sistema massa-mola. A soluo apresentada na gura 1.17.
Observe que, aps a aplicao do impulso, o sistema oscila indenidamente
com uma amplitude de 0,2 m.
Durante a aplicao da fora, o sistema adquire energia cintica e energia
potencial elstica. No entanto, aps a fora ser retirada do sistema, a energia
mecnica conservada:

1 2 1
mvmax = kx2max
2 2
1.3. VIBRAO DE UMA PARTCULA 15

1 ,0
(a )
0 ,5 S in a l 1
In te n s id a d e d o s in a l

0 ,0
1 ,0

S in a l 2
(b )
0 ,5

0 ,0
1 ,0
(c )
0 ,5 S in a l 1 - S in a l 2

0 ,0
0 2 4 6 8 1 0
T e m p o (s )

Figura 1.16: Gerao de um pulso unitrio com duas funes de Heaviside.

Assim, a amplitude de oscilao ser:

1
xmax = vmax (1.8)
0
em que 0 a frequncia natural de oscilao:
s
k
0 = = 5 rad/s
m

A velocidade vmax pode ser calculada pela denio de impulso:

F t = mvmax

em que t = 0,1 s e F =
vmax = 1,0 m/s. Substituindo este re-
10 N. Logo,
sultado na equao 1.8, obtemos xmax = 0,2 m. Com a frequncia natural
de oscilao podemos calcular tambm a frequncia (f 0,8 Hz) e o perodo
(T 1,26 s) de oscilao.
 16 CAPTULO 1. SIMULAO NUMRICA COM MATLAB

0 ,2 (a )
Posio (m)

0 ,1
0 ,0
-0 ,1
-0 ,2
1 0 (b )
Fora (N)

8
6 I = F x t = r e a = 1 0 x 0 , 1 = 1 N s
4
2
0
0 ,0 0 ,5 1 ,0 1 ,5 2 ,0 2 ,5 3 ,0 3 ,5 4 ,0 4 ,5 5 ,0

T e m p o (s )

Figura 1.17: (a) Amplitude do sistema de massa-mola devido a aplicao de

um (b) impulso unitrio.

Podemos estudar tambm o comportamento da energia com os blocos

 Derivative  (disponvel na pasta  Continous ) e  To Workspace . Com es-
tes blocos posicionados no formato da gura 1.18, obtemos o comportamento
da velocidade do corpo. Para calcular as energias potencial elstica, cintica
e mecnica, basta digitar o comando abaixo na janela de comando do MA-
1
TLAB :

plot(tout,[0.5*25*x.*x,0.5*25*v.*v,0.5*25*x.*x+0.5*25*v.*v)

1 Quando
realizamos os clculos no Simulink, o programa armazena todas as variveis
em matrizes coluna. Cada elemento da coluna representa o valor daquela varivel em
um determinado instante de tempo. Ao calcularmos a energia potencial elstica da mola,
temos que obter o termo x2 . Na viso do MATLAB, esta matriz obtida pelo produto de
duas colunas com a multiplicao realizada termo a termo. Para informar ao MATLAB
este procedimento, devemos representar x2 como x.*x.
1.3. VIBRAO DE UMA PARTCULA 17

2
ou digitar o comando abaixo para gerar o arquivo de texto dados.txt :

dados=fopen('dados.txt','wt');
for i=1:length(tout)
K(i)=0.5*1.0*v(i)*v(i);
U(i)=0.5*25*x(i)*x(i);
fprintf(dados,'%f %f %f %f\n',tout(i),K(i),U(i),K(i)+U(i));
end
fclose(dados);

Step F

To Workspace2

1 1
10 1
s s
Scope
Gain2 Integrator2 Integrator3 Gain1
1
Integrator x
Step1 s
To Workspace

1
Integrator1 du/dt v
s

Gain Derivative To Workspace1

25

Figura 1.18: Diagrama de blocos na forma direta 1 para o sistema massa-

mola com a implementao do clculo da velocidade (verso Simulink).

O comportamento da energia cintica, potencial elstica e mecnica est

representado na gura 1.19. Devido ao impulso, o sistema adquire, inicial-
mente, energia cintica e energia potencial elstica. Aps a fora externa ser
retirada, a energia mecnica conservada (U + K = constante). Observe que
a fora e o intervalo de integrao podem ter os mais variados comportamen-
tos, desde que a rea abaixo da curva seja unitria. Neste caso, o impulso
poderia, inclusive, ser representado por uma funo delta de Dirac.

1.3.2 Vibrao forada

Para estudar a vibrao forada, aplicamos os mesmos procedimentos. A

nica diferena est no sinal de entrada. Para forar a vibrao, podemos
usar os blocos da pasta  Sources . Na seo anterior usamos o bloco  Sine
2 Neste caso o MATLAB interpreta o termo x2 como um produto simples, assim no
preciso declarar o produto como zemos no comando plot.
18 CAPTULO 1. SIMULAO NUMRICA COM MATLAB

Energia total

(a)
0 ,6
Energia (J)
0 ,4
Energia cintica
0 ,2 Energia potencial
elstica
0 ,0

1 0 (b)
8
Fora (N)

6 I = F x t = rea = 10 x 0,1 = 1 Ns
4
2
0
0 ,0 0 ,5 1 ,0 1 ,5 2 ,0 2 ,5 3 ,0

Tempo (s)

Figura 1.19: (a) Energia do sistema devido a aplicao de um (b) impulso

mecnico de 1 Ns.

Wave ; desta vez usaremos o bloco  Repeating Sequence  que dene o sinal
no formato de uma onda dente de serra. O diagrama para a vibrao forada
est na gura 1.20. Com criatividade, voc pode montar o seu prprio sinal
de entrada assim como zemos para montar um pulso com dois degraus. O
perodo de cada dente da serra foi denido como 2,0 s. O comportamento
da energia do sistema, amplitude de oscilao e a fora externa esto re-
presentados na gura 1.21. Como o sistema no conservativo, a energia
mecnica uma funo do tempo.

1.4 Animao

Nesta seo, estudaremos a animao de um sistema massa-mola com

os dados obtidos no Simulink. A animao ser criada com um loop de
vrios grcos. Cada gura possui apenas um ponto que representa a posio
instantnea da partcula. Assim, se a simulao realizada em um intervalo
1.4. ANIMAO 19

To Workspace2

1 1
10 1
s s
Scope
Repeating Gain2 Integrator2 Integrator3 Gain1
Sequence
1
Integrator x
s
To Workspace

1
Integrator1
s du/dt v
Gain
Derivative To Workspace1
25

Figura 1.20: Diagrama de blocos na forma direta 1 para a vibrao forada

de uma partcula (verso Simulink).

6 ,0
(a )
Energia (J)

E m e c
= K + U
4 ,0

K U
2 ,0

0 ,0 Posio (m)
(b ) 0 ,5

0 ,0

-0 ,5

(c )
Fora (N)

1 0 ,0

5 ,0

0 ,0

0 2 4 6 8 1 0
T e m p o (s )

Figura 1.21: (a) Energia e (b) posio de um sistema massa-mola excitado

por uma (c) fora externa perodica no formato de dente de serra.
20 CAPTULO 1. SIMULAO NUMRICA COM MATLAB

de tempo de 10 s com passo de 0,1 s, haver um loop com 100 grcos .

3

Aps realizar a simulao no Simulink e exportar os dados para a janela de

comando no MATLAB, abra um arquivo .m e salve-o. Nele, digite:

clc
figure
for i=1:length(tout)
plot(x(i),0,'o','MarkerSize',40,'MarkerFaceColor','r');
axis([min(x),max(x),-2,2]);
xlabel('Posio x (m)');
ylabel('Posio y (m)');
mov(i)=getframe(gcf);
end
movie2avi(mov,'myfirstmovie.avi','compression','None');

O comando length(tout) mede o comprimento da matriz tout. Assim,

se este vetor possui 100 linhas (100 instantes de tempo), o comando for rea-
lizar um loop com 100 ciclos e, portanto, 100 grcos. O comando plot est
representando a coordenada x(i) da partcula no eixo horizontal e a coorde-
nada y(i)=0 no eixo vertical. As demais funes deste comando representam
a forma geomtrica do ponto ('o') com dimetro de 40 unidades relativas
('MarkerSize',40) e face na cor vermelha ('MarkerFaceColor','r'). O
comando axis([x0,x,y0,y]) trava os eixos horizontal e vertical nos valo-
res declarados. Observe que no eixo horizontal, foram utilizados os comandos
min(x) e max(x). Estes comandos encontram o menor e o maior valor da
varivel x. Os comandos xlabel e ylabel permitem nomear os eixos do
grco. O cdigo responsvel por capturar as guras e transform-las em
uma animao o getframe. Dentro deste comando, a funo gcf signica
 get current gure  e ela a responsvel por gravar as informaes da gura
atual. Os detalhes de cada gura armazenada na varivel mov(i) que ser
utilizada para produzir a animao. Para gerar o vdeo, utilizamos o co-
mando movie2avi. Este comando converter o conjunto de grcos em um
vdeo com extenso avi. No primeiro campo deste comando, deve ser infor-
mado o vetor que armazena os grcos do loop (mov). Em seguida, deve ser
informado o nome do arquivo ('myfirstmovie.avi') e o mtodo de com-
presso do vdeo ('compression','None'). Se voc no possui uma boa
placa de vdeo e a animao possui muita informao, sugerido comprimir
os arquivos (visite a pgina ocial do MATLAB para maiores informaes).
No meu computador, o vdeo gerado possui 70 Mb (um notebook velhinho).
Em verses mais atuais do MATLAB, o comando movie2avi foi substitudo
por VideoWriter que permite a converso para, tambm, MP4. A gura

3 Para
criar a animao, no recomendado usar um passo de integrao pequeno,
exceto se a congurao do computador for adequada.
 1.4. ANIMAO 21

1.22 mostram dois quadros da animao. No foi necessrio utilizar com-

pressor, pois a animao possui apenas um ponto. Para deix-la mais real,
precisamos colocar a mola.

2 2

1.5 1.5

1 1

0.5 0.5
Posio y (m)

Posio y (m)
0 0

-0.5 -0.5

-1 -1

-1.5 -1.5

-2 -2
-0.2 -0.1 0 0.1 0.2 0.3 0.4 0.5 0.6 -0.2 -0.1 0 0.1 0.2 0.3 0.4 0.5 0.6
Posio x (m) Posio x (m)

Figura 1.22: Posies de uma partcula em um sistema massa-mola. As

guras mostram 2/100 quadros gerados na animao.

FIGURA 22. POSIES DE UMA PARTCULA EM UM SISTEMA MASSA-MOLA. AS FIGURAS MOSTRAM 2/100 QUADROS GERADOS NA ANIMAO.
Cada espira da mola ser representada por trs segmentos de reta (seg-
mentos s2 ,espira
Cada s3 e s4da
) conforme mostra a gura 1.23 com a mola na posio de
mola ser representada por trs segmentos de reta (segmentos
(xi = 0). Se a mola possui uma nica espira, sero utilizados
s2,relaxamento
s3 e s4) conforme mostra a figura 23 com a mola na posio de relaxamento (xi =
mais dois segmentos horizontais para xao (segmentos s1 e s5 ). As coor-
0). Se a mola possui uma nica espira, sero utilizados mais dois segmentos
denadas x
horizontais 5 e xfixao
para 4 so xas, pois a coordenada
(segmentos x5 est presax5em
s1 e s5). As coordenadas e xuma superfcie
4 so fixas, pois a
imvel. A partcula est na posio
coordenada x5 est presa em uma superfcie x s
i imvel. A partcula est na1 .posio
e conectada ao segmento As co-x e
i
ordenadas
conectada no dos
segmento s1. Ass1coordenadas
segmentos at s4 so dos
mveis, pois devero
segmentos s1 at s4acompanhar
so mveis, opois
movimento
devero da partcula
acompanhar durante
o movimento daapartcula
animao. Os comprimentos
durante a animao. Osde todos os
comprimentos
desegmentos so constantes.
todos os segmentos so constantes.

(-x3, y3)
Parede s4
fixa s5 s3 (-x1, 0) s1
(-x5, 0) (-x4, 0) s2 (xi, 0) x

(-x2, -y2)
FIGURA 23. ESPIRA DE UMA MOLA PROJETADA EM DUAS DIMENSES.
Figura 1.23: Espira de uma mola projetada em duas dimenses.

As coordenadas x e y dos pontos da figura 23 sero obtidas com o auxlio da

figuraAs Considerandoxuma
24.coordenadas e y dos xi gura
pontos da
deformao 0 na mola,
1.23 o
sero obtidas com total
seu comprimento o au-ser
D xlio
+ xi =da
2d gura
+ 4w. Assim, a cateto horizontal
1.24. Considerando umawdeformao xi 6=
dos segmentos s2, 0s3na
e s4mola,
dadoo por:
seu
comprimento total ser D + xi = 2d + 4w. Assim, a cateto horizontal w dos
1
= 4 ( 2 + ) (1)

e o cateto vertical ser:

= 2 2 (2)

em que H a comprimento do segmento. Para calcular h devemos definir um valor

mnimo para H que ser obtido quando a partcula atingir sua posio xi(mx). Quando
 22 CAPTULO 1. SIMULAO NUMRICA COM MATLAB

segmentos s2 , s3 e s4 dado por:

1
w = (D 2d + xi ) (1.9)
4
e o cateto vertical ser:

h= H 2 w2 (1.10)

em que H a comprimento do segmento. Para calcular h devemos denir

um valor mnimo para H que ser obtido quando a partcula atingir sua
posio xi(max) . Quando a mola estiver completamente tracionada (que ser
considerado, por ora, quando todos os segmentos estiverem na horizontal),
4
Hmin ser dado pela equao D + xi(max) = 2d + 4Hmin ou :

1
Hmin = (D 2d + xi(max) ) (1.11)
1 4
min = ( 2 + (mx) ) (3)
4

y
D xi

Parede
fixa
2h
x

h
d w 2w w d

FIGURA 24. DIMENSES DA ESPIRA DE UMA MOLA PROJETADA EM DUAS DIMENSES.

Figura 1.24: Dimenses da espira de uma mola projetada em duas dimenses.

A equao (3) no simula a condio de uma mola real, pois a mola, neste
caso, Ater a aparncia de uma deformao plstica quando todos os segmentos
equao 1.11 no simula a condio de uma mola real, pois a mola,
estiverem na horizontal. Para corrigir o problema, a equao (3) deve ser:
neste caso, ter a aparncia de uma deformao plstica quando todos os
segmentos estiverem na horizontal. Para corrigir o problema, a equao 1.11
1
deve
= 4ser:
( 2 + (mx) ) (4)
1
Hmin = C (D 2d + xi(max) ) (1.12)
em que C um fator de correo que 4trataremos mais adiante. Logo, com as equaes
em
(1), (2)que C eaum
e (4) fator deinicial
declarao correo que trataremos
das variveis D, d e omais
fator adiante. Logo,
de correo com
C, podemos
as equaes
calcular 1.9, 1.10 dos
as coordenadas e 1.12 e a declarao
segmentos s1 at sinicial das variveis
5 em funo D, xdi edao partcula
da posio fator
(veja a tabela C1).
de correo Com estes
, podemos dados as
calcular coordenadas
possvel inserir mola coms1uma
dosasegmentos s5 em na
at espira
funo da posio xi da partcula (veja a tabela 1.1). Com estes dados
animao:
possvel inserir a mola com uma espira na animao:
clc
figure4 A equao 1.11 foi obtida a partir da equao 1.9, pois quando todos os segmentos
for i=1:length(tout)
esto na horizontal,dew construo
% Parmetros = Hmin da mola
D=1.5; % Comprimento total (m)
d=0.5; % Comprimento dos conectores (m)
C=1.2; % Fator de correo
w=0.25*(D-2*d+x(i)); % Cateto horizontal de um segmento da espira
H=C*0.25*(D-2*d+max(x));% Comprimento de um segmento da espira
h=sqrt(H^2-w^2); % Cateto vertical de um cateto da espira
% Coordenadas dos segmentos da mola
x1=x(i)-d; y1=0;
x2=x1-w; y2=-h;
x3=x2-2*w; y3=h;
x4=-D+d; y4=0;
1.4. ANIMAO 23

clc
figure
for i=1:length(tout)
% Parmetros de construo da mola
D=1.5; % Comprimento total (m)
d=0.5; % Comprimento dos conectores (m)
C=1.2; % Fator de correo
w=0.25*(D-2*d+x(i)); % Cateto horizontal de um segmento da espira
H=C*0.25*(D-2*d+max(x));% Comprimento de um segmento da espira
h=sqrt(H^2-w^2); % Cateto vertical de um cateto da espira
% Coordenadas dos segmentos da mola
x1=x(i)-d; y1=0;
x2=x1-w; y2=-h;
x3=x2-2*w; y3=h;
x4=-D+d; y4=0;
x5=-D; y5=0;
% Representao grfica da mola
plot([x(i),x1],[0,0],'k','LineWidth',2); hold on; % s1
plot([x1,x2],[y1,y2],'k','LineWidth',2); % s2
plot([x2,x3],[y2,y3],'k','LineWidth',2); % s3
plot([x3,x4],[y3,y4],'k','LineWidth',2); % s4
plot([x4,x5],[y4,y5],'k','LineWidth',2); % s5
% Representao grfica da partcula
plot(x(i),0,'o','MarkerSize',40,'MarkerFaceColor','r'); hold off;
% Definio dos eixos horizontal e vertical
axis([-D,max(x),-2,2]);
xlabel('Posio x (m)');
ylabel('Posio y (m)');
% Captura da imagem
mov(i)=getframe(gcf);
end
% Converso das imagens no filme
movie2avi(mov,'mysecondmovie.avi','compression','None');

Tabela 1.1: Coordenadas dos pontos da gura 1.23.

Coordenada x Coordenada y
x1 = xi d y1 = 0
x 2 = x1 w y2 = h
x3 = x2 2w y3 = h
x4 = D + d y4 = 0
x5 = D y5 = 0

O resultado est representado na gura 1.25. Observe que temos seis

comandos plot no cdigo. Os primeiros cincos so responsveis por repre-
sentar a mola e o ltimo responsvel por representar a partcula. Em um
loop com apenas um comando plot, o MATLAB faz o grco solicitado pelo
comando atual e apaga a representao grca do comando anterior. Foi
assim que zemos a animao da gura 1.22. Porm, quando h mais de
24 CAPTULO 1. SIMULAO NUMRICA COM MATLAB

um comando plot durante um ciclo do loop, o MATLAB enxerga apenas

o ltimo comando plot e ignora os demais. Para evitar este problema, usa-
mos o comando hold on (ou hold all) logo aps o primeiro plot (veja o
cdigo). Com isto, todos os plot sero sobrepostos em um mesmo grco.
Aps o ltimo plot do ciclo obrigatrio inserir a funo hold off. Caso
contrrio, os grcos de todos os ciclos sero sobrepostos em uma nica ima-
gem (Se voc no tem uma boa placa de vdeo, no tente fazer isso).
Na organizao das linhas do cdigo, importante que a representao
grca da partcula esteja aps a representao grca da mola. Com isso,
parte do segmento de reta s1 ser sobreposto pela partcula e criar a iluso
de conexo entre estes dois objetos. Alm disso, os eixos do grco devem
possuir aproximadamente a mesma largura para evitar que as espiras mu-
dem o tamanho aparente. Este efeito est evidente na gura 1.25 (observe a
mudana aparente no comprimento do segmento s3 ).

2 2

1.5 1.5

1 1

0.5 0.5
Posio y (m)
Posio y (m)

0 0

-0.5 -0.5

-1 -1

-1.5 -1.5

-2 -2
-1.4 -1.2 -1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 -1.4 -1.2 -1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6
Posio x (m) Posio x (m)

Figura 1.25: Posies de uma partcula em um sistema massa-mola. A mola

possui uma nica espira.

Para simular uma mola com duas espiras, adicionamos uma espira no
ponto x4 da gura 1.23 e criamos mais dois pontos (veja a tabela 1.2). Con-
siderando que o comprimento da mola o mesmo e o nmero de catetos
horizontais foi duplicado, o comprimento da mola ser D + xi = 2d + 8w.
Assim, w representado por:

1
w = (D 2d + xi ) (1.13)
8

Os parmetros h e H permanecem os mesmos. Inserindo a equao 1.13

e os pontos da tabela 1.2 no modelo, obtemos a animao ilustrada na gura
1.26.
1.4. ANIMAO 25

Tabela 1.2: Coordenadas dos segmentos de uma mola com duas espiras.
Coordenada x Coordenada y
x1 = xi d y1 = 0
x 2 = x1 w y2 = h
x3 = x2 2w y3 = h
x4 = x3 2w y4 = h
x5 = x4 2w y5 = h
x6 = D + d y6 = 0
x7 = D y7 = 0

2 2

1.5 1.5

1 1

0.5 0.5
Posio y (m)

Posio y (m)

0 0

-0.5 -0.5

-1 -1

-1.5 -1.5

-2 -2
-1.4 -1.2 -1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 -1.4 -1.2 -1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6
Posio x (m) Posio x (m)

Figura 1.26: Posies de uma partcula em um sistema massa-mola. A mola

possui duas espiras.

Para criar uma mola com n espiras, a equao 1.9 deve ser generalizada
para uma mola com 2nw catetos, assim:

1
w= (D 2d + xi ) (1.14)
4n
e a mola ser representada pelo cdigo:

num=10; % Nmero de espiras da mola

plot([x(i),x(i)-d],[0,0],'k'); hold all;
% Construo da n espiras
for n=1:1:num;
plot([x(i)-d-4*(n-1)*w,x(i)-d-w-4*(n-1)*w],[0,-h],'k',...
[x(i)-d-w-4*(n-1)*w,x(i)-d-3*w-4*(n-1)*w],[-h,h],'k',...
[x(i)-d-3*w-4*(n-1)*w,x(i)-d-4*w-4*(n-1)*w],[h,0],'k');
end
plot([-D+d,-D],[0,0],'k');

Portanto, a animao, ilustrada na gura 1.27, dada por:

26 CAPTULO 1. SIMULAO NUMRICA COM MATLAB

clc
figure
for i=1:length(tout)
% Parmetros de construo da mola
D=2.5; % Comprimento total
d=0.5; % Comprimento dos conectores
C=0.3; % Fator de correo
num=10; % Nmeros de espiras da mola
w=(0.25/num)*(D-2*d+x(i)); % Cateto horizontal de um segmento
H=C*0.25*(D-2*d+max(x)); % Comprimento de um segmento
h=sqrt(H^2-w^2); % Cateto vertical de um segmento
% Representao grfica da mola
plot([x(i),x(i)-d],[0,0],'k','LineWidth',2); hold all;
% Construo de n espiras
for n=1:1:num;
plot([x(i)-d-4*(n-1)*w,x(i)-d-w-4*(n-1)*w],[0,-h],'k',...
[x(i)-d-w-4*(n-1)*w,x(i)-d-3*w-4*(n-1)*w],[-h,h],'k',...
[x(i)-d-3*w-4*(n-1)*w,x(i)-d-4*w-4*(n-1)*w],[h,0],'k',...
'LineWidth',2);
end
plot([-D+d,-D],[0,0],'k','LineWidth',2);
% Representao grfica da partcula
plot(x(i),0,'o','MarkerSize',40,'MarkerFaceColor','r'); hold off;
% Definio dos eixos horizontal e vertical
axis([-D,max(x),-2,2]);
xlabel('Posio x (m)');
ylabel('Posio y (m)');
% Captura da imagem
mov(i)=getframe(gcf);
end
% Converso das imagens no filme
movie2avi(mov,'mygreatmovie.avi','compression','None');

O fator de correo na equao 1.12 depende do nmero de espiras da

mola. Para obter uma relao direta entre estes dois parmetros, igualamos
1.12 com 1.14 quando w = H. Assim:

1
Cmin = (1.15)
n

Observe que nos cdigos anteriores, usamos C = 1,2, para a mola com
uma espira (n = C = 0,3 para a mola com dez espiras (n = 10).
1), e
Nestas duas situaes, Cmin = 1 e 0,1, respectivamente. Quando estes valores
so usados, a mola adquire a aparncia de deformao plstica. Logo, para
simular condies mais realsticas:

1
C> (1.16)
n
1.4. ANIMAO 27

2 2

1.5 1.5

1 1

0.5 0.5
Posio y (m)

Posio y (m)
0 0

-0.5 -0.5

-1 -1

-1.5 -1.5

-2 -2
-2.5 -2 -1.5 -1 -0.5 0 0.5 -2.5 -2 -1.5 -1 -0.5 0 0.5
Posio x (m) Posio x (m)

Figura 1.27: Posies de uma partcula em um sistema massa-mola. A mola

possui dez espiras.

Com o sistema massa-mola inserido, podemos adicionar mais itens e tor-

nar a animao mais informativa. Alm da representao grca do movi-
mento, podemos inserir o comportamento dos grcos de fora, posio ou
energia e os vetores de fora. O cdigo a seguir, ilustrado na gura 1.28,
apresenta o sistema massa-mola, os vetores de fora, a evoluo temporal
da fora externa, responsvel pela vibrao forada, e a fora elstica. O
leitor poder complementar o cdigo conforme desejar. Tudo depender da
imaginao!

clc
figure
for i=1:length(tout)-1
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Grficos da fora em funo do tempo %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
subplot(2,1,1)
k = 25; % Constante da mola
% Grficos
plot([tout(i),tout(i+1)],[F(i),F(i+1)],'b'); hold all;
plot([tout(i),tout(i+1)],[-k*x(i),-k*x(i+1)],'k');
% Formatao dos eixos
axis([0,10,-20,10]);
xlabel('Tempo (s)');
ylabel('Fora (N)');
legend('Externa - F_{ext}','Elstica - F_e');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Representao grfica do sistema massa-mola %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
subplot(2,1,2)
% Parmetros de construo da mola
D=2.5; % Comprimento total
28 CAPTULO 1. SIMULAO NUMRICA COM MATLAB

d=0.5; % Comprimento dos conectores

C=0.4; % Fator de correo
num=15; % Nmeros de espiras da mola
w=(0.25/num)*(D-2*d+x(i)); % Cateto horizontal de um segmento
H=C*0.25*(D-2*d+max(x)); % Comprimento de um segmento
h=sqrt(H^2-w^2); % Cateto vertical de um segmento
% Representao grfica da mola
plot([x(i),x(i)-d],[0,0],'k'); hold all;
for n=1:1:num;
plot([x(i)-d-4*(n-1)*w,x(i)-d-w-4*(n-1)*w],[0,-h],'k',...
[x(i)-d-w-4*(n-1)*w,x(i)-d-3*w-4*(n-1)*w],[-h,h],'k',...
[x(i)-d-3*w-4*(n-1)*w,x(i)-d-4*w-4*(n-1)*w],[h,0],'k');
end
plot([-D+d,-D],[0,0],'k');
% Vetores de fora
% Fora externa
plot([x(i),x(i)+F(i)/max(F)],[0,0],'b','LineWidth',2);
plot([x(i)+F(i)/max(F),x(i)+F(i)/max(F)*(1-0.1)],...
[0,F(i)/max(F)*0.2],'b','LineWidth',2);
plot([x(i)+F(i)/max(F),x(i)+F(i)/max(F)*(1-0.1)],...
[0,-F(i)/max(F)*0.2],'b','LineWidth',2);
text(x(i)+F(i)/(max(F)),0.5,'\fontsize{15} F_{ext}',...
'HorizontalAlignment','left','VerticalAlignment','middle',...
'Color','b');
% Fora elstica
plot([x(i),x(i)-k*x(i)/max(F)],[0,0],'k','LineWidth',2);
plot([x(i)-k*x(i)/max(F),x(i)-k*x(i)/max(F)*(1-0.1)],...
[0,k*x(i)/max(F)*0.2],'k','LineWidth',2);
plot([x(i)-k*x(i)/max(F),x(i)-k*x(i)/max(F)*(1-0.1)],...
[0,-k*x(i)/max(F)*0.2],'k','LineWidth',2);
text(x(i)-k*x(i)/max(F),0.5,'\fontsize{15} F_e',...
'HorizontalAlignment','right','VerticalAlignment','middle');
% Representao grfica da partcula
plot(x(i),0,'o','MarkerSize',40,'MarkerFaceColor','r'); hold off;
% Formatao dos eixos
axis([-D,max(x)+1,-2,2]);
xlabel('Posio x (m)');
ylabel('Posio y (m)');
%%%%%%%%%%%%%%%%%%%%%
% Captura da imagem %
%%%%%%%%%%%%%%%%%%%%%
mov(i)=getframe(gcf);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Converso das imagens no filme %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
movie2avi(mov,'mybestmovieofalltime.avi','compression','None');

Para inserir dois grcos na animao, usamos o comando subplot(m,n,p)

em que m o nmero de linhas, n o nmero de colunas e p a posio de
um grco na matriz mn. Os vetores de fora foram programados para
mudar o tamanho de acordo com a sua intensidade no grco fora versus
1.4. ANIMAO 29

10 10
Externa - Fext Externa - Fext
5 Elstica - Fe 5 Elstica - Fe

0 0
Fora (N)

Fora (N)
-5 -5

-10 -10

-15 -15

-20 -20
0 1 2 3 4 5 6 7 8 9 10 0 1 2 3 4 5 6 7 8 9 10
Tempo (s) Tempo (s)

2 2

1 1
Posio y (m)

Posio y (m)
Fe Fext Fe Fext
0 0

-1 -1

-2 -2
-2.5 -2 -1.5 -1 -0.5 0 0.5 1 1.5 -2.5 -2 -1.5 -1 -0.5 0 0.5 1 1.5
Posio x (m) Posio x (m)

Figura 1.28: 2/100 quadros da animao de um sistema massa-mola com

vibrao forada.

tempo; os comprimentos foram normalizados em relao fora externa e

para identic-los foi utilizado o comando:

text(x,y,'Nome','HorizontalAlignment','left','VerticalAlignment',...
'middle','Color','b').

A partir deste modelo, possvel simular e animar problemas mais ela-

borados. Utilizando os passos anteriores, podemos criar um oscilador massa-
mola acoplado com vibrao forada. Neste caso, o sistema possui trs molas
e duas partculas (gura 1.29). Os grcos da velocidade  versus  posio
mostram a periodicidade (ou no) do movimento.

3 3

2 2

F1 F2 F1 F2
posio y (m)

posio y (m)

1 1

0 0

-1 p1 p2 -1 p1 p2

-2 -2

-3 -3
-1.5 -1 -0.5 0 0.5 1 1.5 2 2.5 3 -1.5 -1 -0.5 0 0.5 1 1.5 2 2.5 3
posio x (m) posio x (m)

Partcula 1 Partcula 2 Partcula 1 Partcula 2

1.5 1 1.5 1
1 1
velocidade (m/s)

velocidade (m/s)

velocidade (m/s)

velocidade (m/s)

0.5 0.5
0.5 0.5
0 0 0 0
-0.5 -0.5
-0.5 -0.5
-1 -1
-1.5 -1 -1.5 -1
-2 -2
-0.5 0 0.5 1.2 1.4 1.6 1.8 -0.5 0 0.5 1.2 1.4 1.6 1.8
posio x (m) posio x (m) posio x (m) posio x (m)

Figura 1.29: Sistema massa-mola acoplado com vibrao forada.

Uma das aplicaes do sistema massa-mola acoplado a suspenso veicu-

1
lar passiva (gura 1.30). Neste exemplo, m = mroda + mpneu e M = Mcarro .
4
30 CAPTULO 1. SIMULAO NUMRICA COM MATLAB

1.8 1.07 1.8 1.07

1.06 1.06
1.6 1.6

Posio de M (m)

Posio de M (m)
1.05 1.05

1.4 1.04 1.4 1.04

1.03 1.03
1.2 1.2
1.02 1.02
M M
1 1
Posio (m)

Posio (m)
0 5 10 15 20 0 5 10 15 20
Tempo (s) Tempo (s)
0.8 0.8

0.6 0.6 0.6 m 0.6

m 0.58 0.58
Posio de m (m)

Posio de m (m)
0.4 0.4
0.56 0.56

0.54 0.54
0.2 0.2
0.52 0.52

0 0.5 0 0.5

0.48 0.48
-0.2 -0.2
3 3.2 3.4 3.6 3.8 0 5 10 15 20 10.2 10.4 10.6 10.8 11 0 5 10 15 20
Posio (m) Tempo (s) Posio (m) Tempo (s)

Figura 1.30: Simulao de suspenso veicular passiva com um modelo de um

quarto de carro.

Os acessrios menores (mola e amortecedor) representam o pneu do carro e

os maiores representam a suspenso. Os grcos mostram as posies dos
corpos em funo do tempo de percurso. Observe que enquanto m possui
um deslocamento mximo de aproximadamente 14 cm, M possui um des-
locamento mximo de aproximadamente 7 cm, caracterizando o efeito do
amortecimento.
Referncias Bibliogrcas
[1] B. Girod, R. Rabenstein, and A. Stenger. Sinais e sistemas. LTC, 2003.

31