Saturn V Simulation Program

Turbo Pascal Program Files
apollo.pas Saturn V simulation main program
earth.pas Earth constants
rocket.pas General purpose trajectory simulation procedures
const.pas Saturn V constants
stage1.pas S-IC procedure
stage2.pas S-II procedure
stage3.pas S-IVB procedure
apollo.exe DOS executable

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

I have written a general purpose simulation program that simulates the Saturn V launch vehicle from the Apollo program. 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. This 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 "apollo" 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   0.0     0       0         0  90.00   0.00  90.00      0.0 2913754.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 -0.1 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 Centre Engine Cutoff 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 initially, changed to 1.5 on S-II Centre Engine Cut-Off, and then set to 1.0 at the beginning of the first S-IVB burn. 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 "apollo" needs this algorithm is that the S-II 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.

Once you have reached orbit (nominally a 185.2 km circular orbit) the S-IVB is fired again with a zero angle of attack. You are then told how many kilograms (or seconds) of propellant is left in the S-IVB. I have been able to get into a 185.2 x 185.3 km orbit with 2915 kg (13.7 seconds) left after the second S-IVB burn. This compares to 17.1 seconds using the data from the Apollo 14 Flight Evaluation Report. So, my software seems to give a pessimistic result. 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 Space Shuttle simulation program is also available. If you have any questions about this software, please let me know.

24 Jul 1995 Update

I have made a number of improvements to the software since the 12 Jun 1995 release.

rocket.h (rocket.pas)
The speed of sound, vs, is now calculated as a function of altitude so that the correct Mach number can be determined. Function cd now inputs the Mach number. Procedure atmosphere now calculates vs as well as the kinematic coefficent of viscosity, Vk. Vk is not used yet, but could be used in cd. The speed and height files are now against seconds.

earth.h (earth.pas)
G and M changed to reflect values given in sci.space FAQ. Rs and Vas added to allow determination of vs and Vk.

saturnV.h, saturnVs2.h, saturnVs3.h (stage1.pas, stage2.pas, stage3.pas)
Many values changed and added to reflect expected flight of Apollo 14 (for which I have the Flight Evaluation Report, FER). Biggest change is the thrust and v_e of the F-1. I previously used the nominal thrust (for standard air pressure and pump inlet conditions) of the F-1 from Apollo 15 on (6770 kN SL). The actual thrust is much higher (as given by the FER). The v_e is also higher using the thrust and propellant rate in the FER. I have also put the expected values of the J-2 (slighty worse than before). For Apollo 14, we have Centre Engine Cut-Off (CECO) and then a change in mixture ratio to 4.8:1. The SIV-B also has a change in the mixture ratio, changing to 4.5:1 at the beginning of the SIVB second burn and then changing back to 5.0:1. The S-IC/S-II interstage is held on by the S-II for 28.3 s (I previously thought that the interstage was released on S-II mainstage). LES_h slightly increased to what I measured in FER. There are probably some other small changes that I have forgotten about. All these changes resulted in a 2.1 s increase in the remaining propellant mass from the last version.

I used a turn time of 8.3 s and a maximum angle of attack of 16.49 degrees to obtain a 185.2x185.3 km orbit.

9 Jan 1996 Update

rocket.h (rocket.pas)
I found a bug which caused the rocket to diverge when the vertical speed (h1) became negative, thus making the program difficult to use. The fix allows you to go into reasonable orbits much easier. The angle of attack (alpha), acceleration (a), and dynamic pressure (Pq) are now properly calculated at the end of the trajectory program. The inertial speed (vi) is now a global variable. I have swapped the alpha and angle variables that I used in the previous version so that alpha becomes the global variable.

saturnVs3.h (stage3.pas)
The variable vit is removed and replaced with the global variable vi. Since trajectory calculates vi, all previous vit calculations are removed.

29 Oct 1996 Update

rocket.h (rocket.pas)
I have modified the program so that you can specify dF (change of thrust per second, N/s) and dRp (change of propellant rate per second, kg/s^2). This allows you to accurately model thrust buildup, decay, and change of thrust as a straight line continuous change, instead of an unrealistic discrete change. I have extensively modified saturnVs[1-3].h so that thrust changes use this feature. The total mass of the rocket is now given in the output.

saturnVs2.h (stage2.pas)
The burn time of the ullage rocket motor has been increased so that burnout now occurs after mainstage.

saturnVs3.h (stage3.pas)
The jettison of the ullage motor rocket casing is now included in the first burn, as well as the thrust decays of the first and second burns. The total mass of the stage has also been corrected.

9 Jun 1997 Update

Converted files from Sun pascal to TMT pascal. Added ability to send output to a file.

24 Dec 1997 Update

Corrected bug in rocket.pas which produced extra lines in an output file. Changed names of output files from "height" and "speed" to "height.dat" and "speed.dat", respectively. Corrected bug in height.gnu which produced an empty postscript file.

17 Aug 2006 Update

Compiled program with Free Pascal.