Sunday, March 14, 2010

AP Physics: Graphing in Visual Python

We're going to spend some time this week studying air resistance. Mathematically, the concept can be a little tricky, since the standard Newtonian principle that changes in velocity depend on net forces is complicated by the fact that the net force now itself depends on velocity!

A ball in free-fall obeys the 2nd Law equations

F_net = ma = m(dv/dt)

F_net = F_gravity + F_air resistance F_gravity = -mg ("-" meaning downward in this case)
F_air_resistance = -bv^n ("-" meaning opposite to the direction of velocity in this case)

So we get what's called a differential equation:

m(dv/dt) = -mg -bv^n

We're going to spend some time solving this analytically but it is also useful to simulate it in Visual Python. The following graphs were made using the Visual Python graphing module:




(These are for b = 0.4 and n = 1.0 -- I chose pretty steep b so we could see the curve in this limited time span.) Note that if the ball had been dropped from a higher altitude, the acceleration curve would have continued to approach zero before the ball hit the ground -- this is known as terminal velocity.

Here's the Visual Python code I used -- feel free to use as you like:

from __future__ import division # makes sure is decimal and not integer
from visual import *
from visual.graph import * # graphing capability

# cast of characters
floor = box(length=30, height=0.5, width=30, color=color.blue) # the apparent ground
ball = sphere(pos=(-15,20,0), radius = 0.5) # distances in meters

# initial conditions
ball.velocity = vector(0,0,0) # m/s
F_fric = vector(0,0,0) # N

# constants of simulation - note upward and rightward are + (positive) in sign
g = 9.8 # magnitude of gravitational field in m/s^2
m = 1 # mass kilograms
dt = 0.01 # time step in seconds
b = 0.4 # drag force constant, where F_drag = -bv^n
n = 1.0 # drag force power, where (again) F_drag = -bv^n
coeff_restitution = 0.70 # fractional decrease in speed at each collision
coeff_friction = 0.15 # coefficient of rolling friction between ball and ground
time = 0.00 # time elapsed since start

# set up velocity graph
graph1 = gdisplay(x=0, y=0, width=600, height=400,
title='Velocity vs. Time', xtitle='time (s)', ytitle='v (m/s)',
xmax=7., xmin=0., ymax=20, ymin=-20,
foreground=color.black, background=color.white)
ballv = gcurve(gdisplay = graph1, color = color.black)

# set up acceleration graph
graph2 = gdisplay(x=0, y=0, width=600, height=400,
title='Acceleration vs. Time', xtitle='time (s)', ytitle='a (m/s^2)',
xmax=7., xmin=0., ymax=20, ymin=-20,
foreground=color.black, background=color.white)
balla = gcurve(gdisplay = graph2, color = color.black)

# gather a few seconds of data
while (time <>

# move time forward
rate(100)
time += dt

# forces acting on the object
F_grav = vector(0,-m*g,0) # gravitational force
F_drag = - norm(ball.velocity) * b * math.pow(ball.velocity.mag, n)
F_net = F_grav + F_drag + F_fric
# note friction force is handled below because we need to know if we've reached the floor

# kinematics
a = F_net / m
ball.pos = ball.pos + ball.velocity*dt
if ball.y - ball.radius <= floor.pos.y + (floor.height/2.0): # if ball reaches floor
F_fric = - norm(ball.velocity) * coeff_friction * m * g # only friction if rolling on ground
ball.velocity.y = -ball.velocity.y * coeff_restitution # loses KE on collision; flip y-velocity
ball.pos = ball.pos + ball.velocity*dt # avoids a problem I had with "sticking" -- kluge
else: # if ball is still in the air
F_fric = vector(0,0,0) # no friction with ground
ball.velocity = ball.velocity + a*dt # do ordinary kinematics

# graph the most recent points
ballv.plot(pos = (time, ball.velocity.y))
balla.plot(pos = (time, a.y))




No comments:

Post a Comment