Feeds:
Posts

## Python: 1D mass and spring system with friction

In a set of posts, I show the solution of a few simple 1D and 2D problems in physics. The numerical solutions and plots all are
done using Python. As we go on, the problems will be more and more difficult. Although, I did all that just for fun!

As the simplest problem, I discuss a mass and spring system in one dimension with friction. As you see in the below plot, the system oscillates and the amplitude of the oscillation decreases gradually. That is due to loss of the energy (work of friction force). The rate of decrease of the amplitude is determined by the friction coefficient. I would like to note that although this problem has analytical solution (Haliday, Resnik, Walker, chapters 7 and 8), I used a numerical scheme to solve it (see the Python code below). The next plot shows the variation of the potential, kinetic, and total energies. As expected, when e.g. the kinetic energy is at zero, the potential energy has a maxima. Note that I only used the simplest form of the drag force (-kv). In reality, one may consider terms of higher order. Alternatively, one might take the friction coefficient, k, as a function of velocity. And here is the code. I didn’t manage to put it in a way it should look like. In case you have an idea, let me know.

“””
Program to solve 1D mass and spring system in presence of friction

Equations

F = dp/dt = – kx

import pylab
import numpy
import scipy

from pylab import *
from numpy import *
from scipy import *

#———————–
# initial condition, constants
#———————–
h0 = .50     # initial deviation from equilibrium in m
v0 = 1.0     # m/s
t0 = 0.0     # s

#———————–
# constants
#———————–
time = 10.        # duration in sec
dt = 1.0e-3       # time step in sec
nstep = int(time/dt) # number of time steps

k = 10.0         # spring constant in N/m
m = 1.0          # mass in kg
k1_drag = 0.5    # kg/s

eta = 0.8        # v_in/v_out in an inellastic encounter
Energy = 0.5 * m * pow(v0,2) + 0.5 * k * pow(h0,2)

#———————–
# array definitions
#———————–
h = zeros((nstep,), dtype=float32)  # height in meters
v = zeros((nstep,), dtype=float32)  # velocity in m/sec
a = zeros((nstep,), dtype=float32)  # acceleration in m/sec^2
t = zeros((nstep,), dtype=float32)  # time in sec
K = zeros((nstep,), dtype=float32)  # kinetic energy
V = zeros((nstep,), dtype=float32)  # potential energy
En = zeros((nstep,), dtype=float32) # total energy

#———– initial conditions ———
i = 0
v = v0
h = h0
K = 0.5 * m * pow(v0,2)
V = 0.5 * k * pow(h0,2)
t = t0
En = K + V

for i in arange(1,nstep):#———– main loop ———

a[i-1] = -(k/m) * h[i-1] – (k1_drag/m) * v[i-1]
v[i] = v[i-1] + a[i-1] * dt
dE = – k1_drag * abs(v[i-1]) * abs(v[i-1]) * dt # drag
Energy = Energy + dE
h[i] = h[i-1] + v[i-1]*dt
t[i] = t[i-1] + dt
K[i] = 0.5 * m * pow(v[i],2)
V[i] = 0.5 * k * pow(h[i],2)
En[i] = K[i] + V[i]

tz = [0, time]
hz = [0, 0]
#——————————————————
pylab.figure()
plot(t, h, tz, hz,’k’)
xlabel(‘Time (s)’, fontsize=18)
ylabel(‘Amplitude (m)’, fontsize=18)
tex = r’\$F\,=\,-kx\, -\, k_dv\$’
text(5, .4, tex, fontsize=24)
savefig(‘spring_oscil.png’, dpi=100)
#——————————————————
pylab.figure()
plot(t, K, t, V, t, En, ‘r’, tz, hz,’k’)
xlabel(‘Time (s)’, fontsize=18)
ylabel(‘Energy (J)’, fontsize=18)