Spacecraft Flight Dynamics is the study of spacecraft performance, stability, and control. If you ever listen to mission audio from the American Space Program, the flight controller responsible for this particular aspect of the mission has the call-sign “FIDO.” In terms of popular culture, “FIDO” was the guy in the film Apollo 13 who determined that the mission could still make orbit after the second stage center engine cut out prematurely. In that case, Apollo 13 was able to make orbit by extending the burn time of the four remaining J-2 engines.
Every rocket, ranging from simple fireworks to the mighty Saturn V have the same four basic forces acting on them: weight, thrust, lift, and drag. Simulating the behavior of the rocket mathematically requires us to quantify these forces and establish mathematical relationships between them.
Weight – is the gravitational force acting on the rocket. This is distinct from the mass of the rocket and is dependent upon the local gravitational field.
Thrust – is the reaction force generated by the engine and is determined by the mass of propellant burned in the engine and the velocity with which the exhaust gas is expelled through the engine’s nozzle.
Lift – is an aerodynamic force generated by air flowing past the surface of the rocket. Lift is the component of this force that is perpendicular to the oncoming flow direction. This becomes confusing when the rocket’s flight-path is more vertical than horizontal and we begin treating lift as a lateral force rather that the usual vertical force seen in airplanes and helicopters.
Drag – is a similar aerodynamic force, this time exerted in a direction parallel to that of the airflow.
Ultimately, the behavior of our rocket – and its simulation on the Apple ][ – come back to Sir Isaac Newton and his laws of motion. In particular, the second and third law find application here since our rocket’s thrust is dependent on his third law (“For every action, there is an equal and opposite reaction.”) and the motion of the rocket is dependent on the second law:
F = ma
Where F is the sum of all forces acting on the rocket, m is the current mass of the rocket (constantly changing as it expends fuel and rocket stages), and a is the acceleration vector, the instantaneous rate of change of velocity (v), which in turn is the instantaneous rate of change of displacement. Solving for a, acceleration equals the force sum divided by mass. Acceleration is integrated over time to get velocity, and velocity is in turn integrated to get position.
The weight of the rocket is computer by multiplying its mass (in kg) by the acceleration imparted by the local gravitational field. Although gravity varies as the distance from the center of mass changes, according to the inverse square law, we can assume that for objects near the surface of the Earth, the gravitational acceleration is approximately 9.8 ms2.
As we are using metric units to simplify our arithmetic, we should note that while mass is expressed in kilograms, weight (like thrust) will be specified in a derived unit, Newtons.
As mentioned above, thrust is the reaction force generated by the rocket engine and is determined by the mass of propellant burned in the engine and the velocity with which the exhaust gas is expelled through the engine nozzle. The basic equation for finding the thrust of a rocket is:
Where Fn is the thrust (in Newtons), m is the mass of the fuel and oxidizer burned in the engine and ve is the velocity of the exhaust gas.
Without spending too much time on the mathematics, our simulation has to integrate all four of these forces acting on the rocket. The BASIC source code here is a little rough and is still in a very preliminary stage, but is an attempt at simulating the flight dynamics of a rocket attempting to reach orbit. Again, we’re sticking with the early space program… In this case, the weight and general characteristics of our rocket is a Mercury-Atlas similar to MA-6 that put John Glenn into orbit on February 20, 1962.
Launch Date: 02/20/1962 14:47 UTC
Pilot: John Glenn
Launch Vehicle Mass: 120,000 kg
Payload (Spacecraft) Mass: 1,224.7 kg
Approximate Fuel Mass: 114.420 kg
Main Engine (Sustainer) Thrust: 300 kN
Booster Engines - Thrust: 1,300 kN
Main Engine Burn Time: 300 seconds
Booster Engines Burn Time: 134 seconds
Velocity at Orbital Insertion: 7,844 m/s
Altitude at Orbital Insertion: 248 km
Running this simulation on the Apple IIe Emulator (Source code below), it appears we get the time, range, altitude, and velocity approximately correct, but our translation on the X, Y, and Z axes is way off. Apparently, I’m not resolving the forces correctly in all three dimensions…
I’ll review the code and post more later… I need to check my arithmetic!
1000 REM PROGRAM MAIN
1010 GOSUB 1390: REM INITIALIZE STATE VECTORS
1020 GOSUB 1630: REM INITIALIZE LINEAR INCREMENTAL CHANGE
1030 GOSUB 1230: REM DRAW THE PLOT AXIS
1040 GOSUB 1280: REM PLOT THE COURSE
1050 GET A$: IF A$ = "" THEN 1050
1070 REM RESOLVE FORCE IN THREE DIMENSIONS
1090 LET TZ = 90 - T: REM T IS ANGLE THETA FROM 6000
1100 LET TY = A2: REM LAUNCH AZIMUTH FROM 7000
1110 LET TX = T: REM CURRENT PITCH FROM 5000
1120 LET R9 = (22 / 7) / 180: REM PI OVER 180
1130 LET FX = F * COS (TX * R9): REM FORCE ALONG X AXIS
1140 LET FY = F * COS (TY * R9): REM FORCE ALONG Y AXIS
1150 LET FZ = F * COS (TZ * R9): REM FORCE ALONG Z AXIS
1160 LET F1 = SQR (FX ^ 2 + FY ^ 2 + FZ ^ 2): REM CHECK YOUR ANSWERS!
1190 REM CALCULATE ACCELERATION GIVEN THRUST AND WEIGHT
1210 LET A9 = T9 / VM
1230 REM DRAW SCALE
1240 HGR : HOME : HCOLOR= 3
1260 HPLOT 0,159 TO 279,159
1280 REM PLOT INITIAL BOOST PHASE
1290 VTAB 21: HTAB 1: PRINT "INITIAL BOOST PHASE:"
1300 LET PX = 0:PY = 0:R = 159:PI = 22 / 7
1310 FOR T = 0 TO 90
1320 LET X = R * COS (T * (PI / 180))
1330 LET Y = R * SIN (T * (PI / 180))
1340 LET PX = 279 - X:PY = 159 - Y
1350 HPLOT PX,PY
1360 GOSUB 1680
1370 NEXT T
1390 REM INITIALIZE STATE VECTORS
1400 SX = 0: REM POSITION X METERS
1410 SY = 0: REM POSITION Y METERS
1420 SZ = 0: REM POSITION Z METERS
1430 VX = 0: REM X VELOCITY MPS
1440 VY = 0: REM Y VELOCITY MPS
1450 VZ = 0: REM Z VELOCITY MPS
1460 MT = 0: REM MISSION TIME SECONDS
1470 A1 = 0: REM ALTITUDE MSL
1480 R1 = 0: REM DOWNRANGE METERS
1490 V1 = 0: REM NET VELOCITY
1500 FM = 114420: REM FUEL MASS KG
1510 AX = 0: REM ATTITUDE X AXIS
1520 AY = 0: REM ATTITUDE Y AXIS
1530 AZ = 0: REM ATTITUDE Z AXIS
1540 RX = 0: REM RATE CHANGE X AXIS
1550 RY = 0: REM RATE CHANGE Y AXIS
1560 RZ = 0: REM RATE CHANGE Z AXIS
1570 VM = 120000: REM VEHICLE MASS KG
1580 A2 = 57.5: REM LAUNCH AZIMUTH (90 DEG DUE EAST)
1590 LET G = 9.8: REM ONE G OF ACCELERATION
1600 LET FR = 384.1: REM FUEL CONSUMPTION PER SECOND
1610 LET T9 = 160000: REM THRUST IN NEWTONS
1630 REM LINEAR INCREMENTAL CHANGE TO STATE VECTORS
1640 DIM D(4): RESTORE
1650 DATA 320,188000,248000,7844
1660 FOR I = 1 TO 4: READ D(I): LET D(I) = D(I) / 90: NEXT I
1680 REM UPDATE STATE VECTORS
1690 LET MT = MT + D(1)
1700 LET A1 = A1 + D(1) * V1 * COS (T * (PI / 180))
1710 LET R1 = R1 + D(1) * V1 * SIN (T * (PI / 180))
1720 LET V1 = V1 + D(4)
1730 GOSUB 1190
1740 LET F = T9
1750 GOSUB 1070
1760 LET VX = VX + (FX / VM) * D(1)
1770 LET VY = VY + (FY / VM) * D(1)
1780 LET VZ = VZ + (FZ / VM) * D(1)
1790 LET SX = SX + VX * D(1)
1800 LET SY = SY + VY * D(1)
1810 LET SZ = SZ + VZ * D(1)
1820 LET VM = VM - (FM / 90)
1840 VTAB 21: HTAB 1: PRINT "TIME: "; INT (MT)
1850 VTAB 22: HTAB 1: PRINT "RANGE: "; INT (R1)
1860 VTAB 23: HTAB 1: PRINT "ALTITUDE: "; INT (A1)
1870 VTAB 24: HTAB 1: PRINT "VELOCITY: "; INT (V1);
1880 VTAB 21: HTAB 20: PRINT "POS X: "; INT (SX);
1890 VTAB 22: HTAB 20: PRINT "POS Y: "; INT (SY);
1900 VTAB 23: HTAB 20: PRINT "POS Z: "; INT (SZ);
1920 REM DELAY LOOP
1930 FOR N = 0 TO 2000
1940 NEXT N