workshops Sistema Massa-Mola¶
Vamos agora considerar um sistema constituído por uma mola presa por uma das pontas a uma parede e a outra ponta presa a um bloco que se pode mover livremente na direção horizontal.
O nosso objetivo é resolver para a posição do bloco em função do tempo $x(t)$, conhecendo a sua posição e velocidade iniciais.
Para aplicar o método de Euler temos apenas que especificar as forças aplicadas no corpo.
Esperamos que a força exercida pela mola no bloco seja zero numa certa posição de equilíbrio $x_0$, para o qual a mola está no seu comprimento natural. Já quando a mola se encontra esticada ou comprimida esperamos ver uma força que aponte no sentido oposto à perturbação, procurando levar o bloco para a posição de equilíbrio.
Podemos obter facilmente uma força com estas propriedades escrevendo:
$F = -k \Delta x = -k (x-x_0)$
Esta expressão é chamada de Lei de Hooke e assume uma relação linear entre a força e o deslocamento em relação ao equilíbrio, com constante negativa de modo a que a mola puxe o bloco para a posição de equilíbrio.
Usando esta expressão podemos passar à aplicação do método de Euler
#definição de constantes do sistema necessários aos cálculos
M = 1. #massa do bloco
K = 3. #constante da mola
L = 1. #comprimento natural da mola
dt = 0.002 #intervalo infinitesimal de tempo
ttotal = 10 #tempo que desejamos simular
X0 = 1.1 #condições iniciais
V0 = 0.1
lista_t = [0] #listas onde vamos guardar os dados do sistema a cada iteração
lista_x = [X0]
lista_v = [V0]
#método de Euler
for i in range(int(ttotal/dt)):
forca = - K * (lista_x[-1] - L) #o último elemento da lista_x é a posição atual da mola
aceleracao = forca/M
novo_v = lista_v[-1] + aceleracao*dt
novo_x = lista_x[-1] + lista_v[-1]*dt #o último elemento da lista_v é a velocidade atual da mola
lista_x.append(novo_x) #acrescentamos novos valores no fim das listas
lista_v.append(novo_v)
lista_t.append(lista_t[-1]+dt)
#representação gráfica do movimento calculado
import matplotlib.pyplot as plt
plt.figure(figsize=(13,5))
plt.subplot(1,2,1)
plt.plot(lista_t, lista_x,lw=2)
plt.xlabel("t")
plt.ylabel("x")
plt.subplot(1,2,2)
plt.plot(lista_t, lista_v,lw=2)
plt.xlabel("t")
plt.ylabel("v")
plt.show()
Obtemos soluções oscilatórias $-$ a posição oscila em torno do comprimento natural da mola e a velocidade oscila em torno de zero.
Solução Exata¶
A equação do movimento considerado é:
$F = ma = -k(x-x_0) \;\Longleftrightarrow\;\; a = -\frac{k}{m}(x-x_0)$
A solução geral desta equação pode ser obtida matematicamente e é:
$x(t) - x_0 = \Delta x_0 \cdot \cos{(\omega t)} + \frac{v_0}{\omega} \cdot \sin{(\omega t)} $
$v(t) = v_0 \cdot \cos{(\omega t)} - \Delta x_0 \cdot \omega\cdot \sin{(\omega t)}$
onde as constantes $\Delta x_0$ e $v_0$ são o deslocamento inicial em relação à posição de equilíbrio, e a velocidade inicial. A frequência de oscilação é dada por $\omega = \sqrt{\frac{k}{m}}$.
Podemos comparar estas soluções com as obtidos pelo método de Euler.
#Comparação do método de Euler com a solução exata
import matplotlib.pyplot as plt
import math
plt.figure(figsize=(13,5))
#definição de funções para cáculo da solução exata
def x_exato(t):
w = math.sqrt(K/M)
return L + (X0-L)*math.cos(w*t) + V0/w * math.sin(w*t)
def v_exato(t):
w = math.sqrt(K/M)
return -(X0-L)*w*math.sin(w*t) + V0 * math.cos(w*t)
plt.subplot(1,2,1)
plt.plot(lista_t[1::40], lista_x[1::40],lw=0, marker="o",ms=3,label="Euler")
plt.plot(lista_t, [x_exato(t) for t in lista_t],label="Sol. Exata")
plt.legend(loc="upper right")
plt.xlabel("t")
plt.ylabel("x")
plt.subplot(1,2,2)
plt.plot(lista_t[1::40], lista_v[1::40],lw=0, marker="o",ms=3,label="Euler")
plt.plot(lista_t, [v_exato(t) for t in lista_t],label="Sol. Exata")
plt.legend(loc="upper right")
plt.xlabel("t")
plt.ylabel("v")
plt.show()
Verificamos que a solução numérica da equação (com método de Euler) nos permitiu visualizar o movimento do sistema sem conhecimento prévio da solução exata, e que coincide bastante bem com esta durante o tempo simulado.
Forças de Atrito¶
Vamos adicionar um pouco mais de física ao nosso sistema. Considerámos que a única força aplicada no bloco era a força de restituição da mola, dada pela lei de Hooke, e obtivémos soluções que oscilam para sempre no tempo.
Contudo, sabemos que se tentássemos reproduzir este sistema na vida real a massa acabaria por oscilar cada vez menos até parar, uma vez que perde energia no contacto com a mesa. Na descrição Newtoniana do sistema, dizemos que existe uma força de atrito entre a mesa e o bloco que se opõe ao movimento.
Esta força tem então que ser oposta à direção do movimento (oposta à velocidade do bloco) e deve ser mais forte quando o bloco se move mais rapidamente. Então, podemos contruir uma força proporcional à velocidade mas com direção oposta:
$F_{atrito} = -U v$
onde $U$ é um coeficiente de atrito que deve ser positivo.
Vamos voltar a aplicar o método com a adição desta força e avaliar o movimento resultante.
#definição de constantes do sistema necessários aos cálculos
M = 1.
K = 3.
L = 1.
U = 0.3 #<---- acrescentamos este coeficiente
dt = 0.002
ttotal = 30
X0 = 1.1
V0 = 0.1
lista_t = [0]
lista_x = [X0]
lista_v = [V0]
#método de Euler
for i in range(int(ttotal/dt)):
forca = - K * (lista_x[-1] - L) - U*lista_v[-1] #<---- acrescentamos força de atrito
aceleracao = forca/M
novo_v = lista_v[-1] + aceleracao*dt
novo_x = lista_x[-1] + lista_v[-1]*dt
lista_x.append(novo_x)
lista_v.append(novo_v)
lista_t.append(lista_t[-1]+dt)
#representação gráfica do movimento calculado
import matplotlib.pyplot as plt
plt.figure(figsize=(13,5))
plt.subplot(1,2,1)
plt.plot(lista_t, lista_x,lw=2)
plt.xlabel("t")
plt.ylabel("x")
plt.subplot(1,2,2)
plt.plot(lista_t, lista_v,lw=2)
plt.xlabel("t")
plt.ylabel("v")
plt.show()
Como esperado, a introdução da força de atrito no sitema levou a uma decaimento da amplitude da solução oscilatória.