Space Shuttle Simulation Program

Turbo Pascal Program Files
shuttle.pas Shuttle simulation main program
earth.pas Earth constants
rocket.pas General purpose trajectory simulation procedures
const.pas Shuttle constants
stage1.pas Solid Rocket Booster (SRB) procedure
stage2.pas Main Propulsion System (MPS) procedure
stage3.pas Orbital Maneuvering System (OMS) procedure
srb.dat SRB thrust verus time table
shuttle.exe DOS executable

Gnuplot Files
speed.plt Gnuplot file to generate speed v time curve
speed.ps Speed v time postscript file from "shuttle"
height.plt Gnuplot file to generate height v time curve
height.ps Height v time postscript file from "shuttle"

Similar files for the Space Shuttle with liquid rocket boosters are in lrb.zip.

I have written a general purpose simulation program that simulates the Space Shuttle. The heart of the simulation program is rocket.pas which contains procedures used to determine the trajectory of any rocket using the Runga-Kutta fourth order method. I first wrote this for a Saturn V trajectory simulation program. rocket.pas is separate from the stage*.pas procedures which generate the appropriate values (such as thrust F, propellant rate Rp, time increment dt, etc.) for input to rocket.pas. const.pas contains the constants used in the stage*.pas procedures.

This program was compiled using Free Pascal.

When you execute "shuttle" you will get the following response

Enter output filename (return is standard output): 
gs = 9.79930 m/s^2, mu = 398645.0 km^3/s^2. 
At Ht Pt = 22683 Pa, vst = 295.18 m/s

   t      a     vi     h0       r0      alpha  beta  theta     Pq     m0+Me
  sec   m/s^2  m/s   metres   metres     deg    deg   deg      Pa       kg
-----------------------------------------------------------------------------
   0.00   9.8   409       0         0  90.00   0.00  90.00      0.0 2054857.0
Turn time (s)? 

You should then enter the duration in which you want the rocket to pitch over. The pitch angle has been set to -1.0 degrees, but you can adjust this angle in const.pas with name angle1. The rocket then takes off vertically before pitching over for the time specified. This time does not include the time to move the rocket to and from angle1, so even if you set the angle to 0, the rocket will slightly pitch over. The output from the program has in each line

t	time, seconds
a	acceleration (excluding any external forces such as gravity and
        air drag); metres per second per second
vi	inertial speed; metres per second
h0	altitude above planet's surface; metres
r0	range; metres
alpha	thrust angle relative to inertial velocity vector or angle of attack; 
        degrees
beta	velocity angle relative to motionless planet; degrees
theta   velocity angle relative to rotating planet; degrees
Pq	Dynamic pressure; Pascals
m0+Me   the mass of the rocket

The first stage performs a gravity turn keeping the thrust vector and the air velocity vector the same. Thus the air angle of attack is zero. You will also see events such as SSME throttling down and then back up during the burn. At the end of the first stage burn you are asked for the maximum angle of attack. This is relative to the inertial velocity vector. When this angle is greater than the air angle of attack, the angle of attack will gradually increase. This is due to the trajectory algorithm trying to maintain the rate of altitude increase. We have h₀ as the height, h₁ = dh₀/dt is the rate of increase of height, and h₂ = dh₁/dt. Our orbit algorithm has h₂ proportional to sign(h₁)|h₁|^pow. pow is a constant and is set to 2.0. This seems to be not the most optimal algorithm, but I have found that it does a reasonable good job, getting you where you want to go.

The angle of attack will increase until the maximum angle is reached and be maintained there until centrifugal acceleration becomes strong enough. The angle of attack will then naturally decrease. The reason that "shuttle" needs this algorithm is that the MPS doesn't have enough thrust to allow for a zero angle of attack. Well, this is what I found anyway. If anyone has a better way that uses a smaller angle of attack, I would be glad to hear from them.

At the end of stage 2 you should be in an elliptical orbit and at an altitude of between 80 and 90 km. At this altitude and speed, drag is still significant and so the apogee needs to be higher than the nominal 203.7 km (for a 24.9 t payload, 1 t = 1000 kg). After transferring to apogee the apogee should have decreased to 203.7 km due to drag. You then perform the OMS2 burn in a direct ascent trajectory to circularise your orbit. You are then told how many kilograms (or seconds) of propellant is left in the Orbiter. I have been able to get into a 203.7 x 203.9 km orbit with 3200.6 kg left. I leave it to you to determine the turn time and angle of attack. Have fun!

One area I would like to improve is determining the coefficient of drag (cd) versus speed. The values I used were from the Mars Project by Werner von Braun for a 20 m diameter rocket. If anyone can help me out on this, I would greatly appreciate it.

A concern I have is the need for my altitude to be between 80 and 90 km at the end of stage 2. The only trajectory curves I have seen for the shuttle are from a 1979 book where the poorer performance dual OMS burn is illustrated. There, the altitude is around 100 km at burnout.

I would very much like to thank Linder Metts who provided me with vacuum thrust versus time tables and the vacuum exhaust speed of the SRB.

23 Nov 1999 First release
17 Aug 2006 Free Pascal compilation