# 3D animation of rocket taking off from a flat Earth
# Version 6.0 9-25-09
# Plots the trajectory of the rocket and then its distance and speed vs. time
# Includes wind resistance in Version 2.0
# Includes variable mass in Version 3.0
# Includes a fairly accurate model of time dependent thrust in Version 4.0
# Also includes the fall-off in air density with altitude in Version 5.0
# Finally Version 6.0, includes dependence of the force of gravity upon altitude
# and the final fall-off in thrust in the last few seconds of the solid booster's burn
# Lots of useful info from this hobbyist rocketeer site: http://www.braeunig.us/space/propuls.htm
from visual import *
from visual.graph import * # import graphing features
# initialization section
# Set up the values I will need
# space shuttle website with specs at: http://www.braeunig.us/space/specs/shuttle.htm
dt = 1 # timestep in seconds
m = 2.0e6 # mass in kilograms
RE = 6378.e3 # radius of the Earth in meters
Ntimesteps = 124 # total number of timesteps for the simulation
air_density = 1.0 # air density kg/m^3 < --- correct for air pressure falling off with altitude!
drag_coefficient = 0.5 # values of 0.5 to 0.75 quoted in misc. rocketry sources as the appropriate drag for this geometry (0.5 for shuttle)
area = 400 # in m^2
solid_mass=2.*5.e5 # solid propellant mass in kg per booster * 2-- boosters separate after 2 minutes
solid_burn_time = 124. # solid burn time in s -- combined thrust about 30 MN: maximum thrust at launch
dm_dt_solid = solid_mass/solid_burn_time # rate of ejected solid mass
liquid_mass = 7.3e5 # liquid propellant mass in kg
liquid_burn_time = 530. # liquid burn time in s
F_thrust_solid = 29.36e6 # solid rocket booster thrust at maximum
v_exhaust = F_thrust_solid / dm_dt_solid # exhaust speed
v_exhaust_SSME = 3600 # in m/s according to wikipedia
F_thrust_SMME = 1.8e6*3.*0.65# each of 3 SMME's used can generate 1.8 MN of thrust, at 65% after 35 s
dm_dt_SMME = F_thrust_SMME/v_exhaust_SSME
m_SMME = liquid_mass
# data plotted on Wikipedia show a relativley constant thrust profile vs. time; this depends on the fuel geometry and the surface
# area burning at each time. Some solid fuel designs result in relatively constant thrust, others increasing, decreasing. See main reference.
#Launch thrust: 34,677 kN sea level (three SSMEs firing at 104%/1,754 kN each, two SRBs at 14,680 kN each; SSMEs throttle down to 65% by 35 s and up to 104% by 65 s)
#SSME = space shuttle main engines; SRB = solid rocket booster
#SRB Thrust: combined thrust 29.36 MN SL (maximum thrust at launch reducing by 1/3 after 50 s)
#Pressure: 62.1 atm max, 45.0 atm average
oxa=70000 # For plotting, the wireframe dimensions along Cartesian axes
oya=70000
oza=70000
g = 9.8 # acceleration of gravity--constant at first for a flat Earth
v_yo = 0. # initial speed along y axis--zero at launch
v_xo = 0. # later maybe we will make it 2D!
v_zo = 0. # or even 3D!
x_o = 0 # initial value of x, y and z--starts at the origin
y_o = 0
z_o = 0
# create wireframe with dimensions of the pixel sizes along x,y,z for now
square1 = curve(pos=[(oxa,oya,oza),(-oxa,oya,oza),(-oxa,-oya,oza),(oxa,-oya,oza),(oxa,oya,oza)])
square2 = curve(pos=[(oxa,oya,-oza),(-oxa,oya,-oza),(-oxa,-oya,-oza),(oxa,-oya,-oza),(oxa,oya,-oza)])
square3 = curve(pos=[(oxa,oya,oza),(oxa,oya,-oza),(oxa,-oya,-oza),(oxa,-oya,oza),(oxa,oya,oza)])
square4 = curve(pos=[(-oxa,oya,oza),(-oxa,oya,-oza),(-oxa,-oya,-oza),(-oxa,-oya,oza),(-oxa,oya,oza)])
# end of wireframe
rocket = cone (pos=(x_o,y_o,z_o), radius=1000, axis=(0,5000,0), color=color.white) # initialize rocket's position at (xo,yo,zo), pointing along +y
rocket.velocity = vector(v_xo,v_yo,v_zo) # initialize rocket's velocity too
rocket_position=[]
rocket_speed=[]
c = curve(pos=(x_o,y_o,z_o),color=color.red) # plots the rocket's trajectory
n=0
while n < Ntimesteps-10:
rate (20) # screen refresh rate
air_density = exp(-mag(rocket.pos)*0.00011856) # correction for altitude
F_drag = -0.5*air_density*drag_coefficient*area*mag(rocket.velocity)*mag(rocket.velocity)
F_thrust = F_thrust_solid*(1-(F_thrust_solid-14.e6)*(n/(Ntimesteps-10))) + (F_thrust_SMME* m_SMME/liquid_mass)
# thrust varies due to the actual solid booster thrust vs. time variation: 32MN to 14 MN linear decrease over 110s, then sharp falloff to zero to end.
m_SMME = m_SMME - dm_dt_SMME*dt # reduce the liquid propellant available to the SMME's
m = m-(dm_dt_solid+dm_dt_SMME)*dt # reduce the rocket mass by the ejected solid fuel
F_net_y = F_thrust - m*g*(RE/(RE+mag(rocket.pos)))**2 + F_drag
rocket_force = vector(0.,F_net_y,0.)
rocket.velocity = rocket.velocity + (rocket_force/m)*dt # these update rules conform to Euler-Cromer method; they
# conserve energy over a periodic orbits: important point:
rocket.pos = rocket.pos + rocket.velocity*dt # first update velocity, then update position second!
# If switched in order, it looses energy and orbits don't close!
rocket_position.append(mag(rocket.pos))
rocket_speed.append(mag(rocket.velocity))
c.append(pos=rocket.pos)
n=n+1
while n < Ntimesteps:
rate (20) # screen refresh rate
air_density = exp(-mag(rocket.pos)*0.00011856) # correction for altitude
F_drag = -0.5*air_density*drag_coefficient*area*mag(rocket.velocity)*mag(rocket.velocity)
F_thrust = 14.e6*(1-(Ntimesteps -10 - n)/Ntimesteps) + (F_thrust_SMME* m_SMME/liquid_mass)
# thrust varies due to the actual solid booster thrust vs. time variation: 32MN to 14 MN linear decrease over 110s, then sharp falloff to zero to end.
m_SMME = m_SMME - dm_dt_SMME*dt # reduce the liquid propellant available to the SMME's
m = m-(dm_dt_solid+dm_dt_SMME)*dt # reduce the rocket mass by the ejected solid fuel
F_net_y = F_thrust - m*g*(RE/(RE+mag(rocket.pos)))**2 + F_drag
rocket_force = vector(0.,F_net_y,0.)
rocket.velocity = rocket.velocity + (rocket_force/m)*dt # these update rules conform to Euler-Cromer method; they
# conserve energy over a periodic orbits: important point:
rocket.pos = rocket.pos + rocket.velocity*dt # first update velocity, then update position second!
# If switched in order, it looses energy and orbits don't close!
rocket_position.append(mag(rocket.pos))
rocket_speed.append(mag(rocket.velocity))
c.append(pos=rocket.pos)
n=n+1
print "final speed at solid booster separation = ",mag(rocket.velocity)," m/s (measured 1390 m/s)",100*(mag(rocket.velocity)-1390)/1390,"percent error"
# SRB's take the shuttle to 45km and 1390 m/s before separation
# a canned functional plotting program from the vpython manual--use to plot distance and speed vs. time
scene2 = display(title='Graph of position and speed')
#scene2=display() # new screen with the plots
funct1 = gcurve(color=color.cyan)
funct2 = gcurve(delta=0.05, color=color.yellow)
print "cyan (light blue-green) curve = rocket position, yellow curve = 40.*rocket speed vs. time (s)"
n=0
while n < Ntimesteps:
funct1.plot(pos=(n*dt,rocket_position[n])) # rocket altitude vs. time
funct2.plot(pos=(n*dt,40.*rocket_speed[n])) # rocket speed vs. time--scaled to be on same plot as position
n=n+1