![\"\"](\"http:\/\/pages.vassar.edu\/magnes\/files\/2019\/05\/eq1-2.jpg\")
<\/a><\/p>\nIn this article I will explore each term of the equation of motion, covering the various dependencies that must be accounted for in the code, as well as any assumptions or simplifications made.\u00a0 Finally, we will see the results of this simulation for different rocket models.<\/p>\n
I. Thrust Force<\/h4>\n
The majority of all rockets which have entered into orbit have been liquid-fuel rockets, as opposed to solid-fuel rockets (the notable exception being the Space Shuttle, which uses both).\u00a0 For a liquid-fuel rocket, thrust would be controlled during the ascent by a throttle, responding to rocket performance and adjusting thrust over time to control the trajectory.\u00a0 For the sake of the project, we are assuming that the thrust of the model rocket is constant for the duration of the ascent.\u00a0 A mathematical representation of the thrust force could be:<\/p>\n
In other words, the thrust force is a constant up until the time at which the rocket runs out of fuel, after which it is zero.\u00a0 The constant value of the thrust will be known for each rocket model used in the simulation.\u00a0 Therefore:\u00a0 t<\/strong>he thrust force is dependent on time only.<\/strong><\/p>\n It is important to note that this model is purely one-dimensional.\u00a0 In a realistic rocket flight, the rocket would launch perpendicular to the surface, and then rotate by a small degree in flight to begin forming a parabolic trajectory.\u00a0 Therefore the actual vertical component of the thrust would be slightly lower than it is here.\u00a0 However, we can also assume that a rocket flying at any non-zero angle would generate a small amount of aerodynamic lift that would counteract this reduction in vertical thrust.<\/p>\n The drag acting on the rocket can be separated into two components: atmospheric drag and wave drag.\u00a0 Atmospheric drag, a form of parasitic drag, is a damping force opposite to the rocket’s motion due to skin friction between the surrounding air molecules and the surface of the rocket.\u00a0 The equation for drag force due to atmospheric skin friction is:<\/p>\n where\u00a0\u03c1 is the air density (kg\/m3<\/sup>),\u00a0v<\/em> is the velocity of the rocket,\u00a0C<\/em>D<\/sub> is the drag coefficient, a\u00a0A<\/em> is the cross-sectional area of the rocket.\u00a0 The diameter of each rocket is given, so\u00a0A<\/em> will be known.\u00a0 The drag coefficient is difficult to determine without making physical tests on the rocket, but we will assume a value of\u00a0C<\/em>D<\/sub> = 0.75, which represents a flying shape somewhere between a cylinder and a cone.\u00a0 This is a relatively high drag coefficient; airfoils, for comparison, have values of\u00a0C<\/em>D<\/sub> well below 0.05.<\/p>\n Air density is more complicated to calculate, because it will decrease exponentially as the rocket ascends through the atmosphere.\u00a0 For altitudes under 86 km (the troposphere, stratosphere, and mesosphere), the air density can be given by the ideal gas law (Eq. 5). Here,\u00a0R<\/em> is the gas constant (8.314 J\u22c5K\u22121<\/sup>\u22c5mol\u22121<\/sup>) and\u00a0M<\/em>\u00a0is the molar mass of air (0.029 kg\/mol).\u00a0 Atmospheric temperature and pressure are also dependent on altitude, and are given by the barometric equations below.<\/p>\n Here\u00a0T<\/em>0<\/sub> and P<\/em>0<\/sub>\u00a0are\u00a0the baseline temperature and pressure for certain altitude regions,\u00a0L<\/em> is the temperature lapse rate for that region, and g<\/em>0<\/sub> is the acceleration due to gravity at sea level (9.81 m\/s2<\/sup>).<\/p>\n Above 86 km, in the thermosphere, temperature rises more quickly because the atmosphere no longer shields the space from the sun’s rays.\u00a0 Pressure and air density, now very small, drop off with a quartic relationship to altitude.<\/p>\n The coefficients of each term vary in different high-altitude regions.\u00a0 For a table of these coefficients, as well as the values for\u00a0T<\/em>0<\/sub>,\u00a0P<\/em>0<\/sub> and L<\/em>, please refer to this comprehensive mathematical model of the Earth’s atmosphere<\/a>.<\/p>\n To correctly calculate atmospheric drag in my rocket simulation, I needed to write a function that could accept a value for altitude and compute the pressure, temperature, and air density in each atmospheric region. Here is this function<\/a>, which stops distinguishing very small values for air density above\u00a0z<\/em> = 150 km.\u00a0 Running this function for increasing altitude produced the following graphs of temperature, pressure, and air density.<\/p>\n In summary, the code must update values of air density and drag within the for-loop of the Euler-Cromer method, because the\u00a0atmospheric drag depends on both altitude and velocity<\/strong>.<\/p>\n Earlier versions of my project failed to account for wave drag, or the drag force exerted on an aircraft by the pressure due to compressed waves of air.\u00a0 This force comes into play once the rocket exceeds its critical Mach number (typically Mach 0.8), and approaches the sound barrier.\u00a0 When an object approaches the speed of sound, it begins to move faster than the air waves in front of it, causing a rapid rise in drag.\u00a0 The peak in drag at Mach 1 causes the sound barrier, which is the primary obstacle for rockets attempting to leave the Earth’s atmosphere. When I did not account for wave drag in my simulation, rockets were able to reach space up to distances of 1,000,000 meters without difficulty.<\/p>\n Exact calculations of the relationship between drag and Mach number are difficult and rely heavily on the specific shape of the rocket.\u00a0 Furthermore, most information available on calculating wave drag is from the field of aerodynamic engineering, which deals in airplanes, not rockets.\u00a0 One equation, however, the Prandtl-Glauert Rule, can help approximate the effects of transonic drag forces in my rocket simulation.<\/p>\n The use of the Prandtl-Glauert Rule in this context requires several further caveats.\u00a0 First, it is currently accepted that the Prandtl-Glauert Rule is not valid for accurate calculations in the transonic region, only being valid for Mach numbers below 0.8 and between 1.2-5.\u00a0 I am choosing to ignore this disclaimer, because the rule in this region approximates the effect of the sound barrier well enough for the purposes of the project. Also, the Prandtl-Glauert Rule has been shown in experiments to underestimate increases in wave pressure (NACA Report No. 646, 1939).\u00a0 Two updates exist to the rule — the Karman-Tsien Rule, and Laitone’s Rule — which include higher-order terms.\u00a0 Lastly, the explicit definition of the Prandtl-Glauert Rule relates Mach number to the pressure coefficient of an aircraft, not its drag coefficient.\u00a0 Its application here is by analogy, noting the similar effects in equations for the Prandtl-Glauert corrections to coefficients of lift and moment.<\/p>\n The Prandtl-Glauert Rule (Equation 10) states that when dealing with compressible fluids, the aerodynamic coefficients (e.g. pressure, lift, drag, etc) must be divided by\u00a0\u03b2<\/em>, the Prandtl-Glauert Factor, which depends on the freestream Mach number.\u00a0 The Mach number is defined as the ratio of velocity over sound speed, and the sound speed in turn depends on air temperature.<\/p>\n In aerodynamics, in is customary to perform calculations using a frame of reference that moves along with the aircraft.\u00a0 The subscript\u00a0\u221e below the Mach number, velocity, and sound speed (a<\/em>) refers to the “freestream” value of each, meaning the speed of the moving air around the stationary vessel in this frame of reference, but far away from the vessel, where the shape of the aircraft can not yet affect the airflow.\u00a0 In Equation 13,\u00a0\u03b3<\/em> is the ratio of specific heats (about 1.4 on average), and\u00a0R<\/em>sp<\/sub> is the specific gas constant (287\u00a0J\u22c5kg\u22121<\/sup>\u22c5K\u22121<\/sup> for dry air).\u00a0 For Mach numbers higher than 1 (supersonic speeds), the terms under the square root in Equation 11 are reversed.<\/p>\n We have previously established that air temperature depends on altitude, so we can say that\u00a0wave drag depends both on velocity and altitude<\/strong> just as atmospheric drag does.\u00a0 For this reason, I also needed\u00a0a function<\/a> that could update the coefficient of drag at each iteration of the main for-loop.<\/p>\n Running the function at a fixed altitude and temperature for increasing rocket velocities produced the following results for wave drag:<\/p>\n The drag peaks at 347.18 m\/s, the speed of sound at 300 K.\u00a0 Notice that if that code had not corrected it, the drag would have been infinite at Mach 1 following Equations 10 and 11.\u00a0 This is called the Prandtl-Glauert Singularity.<\/p>\n The force due to gravity is given by the same simple equation as in other projectile physics problems:\u00a0\u00a0FG<\/sub> = mg<\/em>.\u00a0 However, in this case, calculating the gravitational acceleration and the mass of the rocket introduces new problems.<\/p>\n The mass of the rocket will always be decreasing as it expels fuel.\u00a0 Continuing to assume constant thrust, and therefore a constant mass flow rate out of the rocket, we can model the mass of the rocket with a linear time dependence as in Equation 15.\u00a0 The gravitational acceleration at an altitude\u00a0z<\/em> from the Earth’s surface is given in Equation 16, where\u00a0RE<\/sub><\/em> is the radius of the Earth (6,371,000 m) and\u00a0ME<\/sub><\/em> is the mass of the Earth (5.972 x 1024<\/sup> kg).<\/p>\n The mass flow rate can be calculated from the Tsiolkovsky Rocket Equation, given the thrust force and the specific impulse of the rocket.\u00a0 The specific impulse, in seconds, is a measure of the efficiency of the rocket’s engines, and is known for most rockets.<\/p>\n Because mass is decreasing over time, and gravity is decreasing with altitude,\u00a0gravitational force depends on time and altitude<\/strong>.<\/p>\n We have now reduced the equation of motion to depend only on time, altitude as a function of time, and velocity (the derivative of altitude) as a function of time.\u00a0 We can now apply the Euler-Cromer method to solve for altitude and velocity and determine the rocket’s one-dimensional trajectory.<\/p>\n My code for solving the equation of motion is attached\u00a0here<\/a> in full.\u00a0 A pseudo-code summary of the program is given below.<\/p>\n I would like to call attention to the three “if conditions” in the for-loop, which were inserted as the result of trial and error.\u00a0 First, the “drag vector condition” is necessary because the code cannot determine the direction component of the velocity vector, only its magnitude.\u00a0 If the rocket fails, and begins to fall, the other forces will still point upwards, but the drag force will need another line of code to tell it to change sign when the velocity is negative.<\/p>\n The “crash condition” simply says to stop the for-loop and end the simulation when the rocket hits the ground, when z = 0 or z < 0.\u00a0 Otherwise, the rocket would fall through the Earth.\u00a0 This is also why the rocket’s altitude was initialized at z = 1, so as not to prematurely trigger this condition.<\/p>\n Finally, the “empty condition” says to set the thrust and mass flow rate equal to zero when the rocket runs out of fuel. The empty mass (“dry mass”) of the rocket was initialized, so the condition simply checks to see when\u00a0m<\/em> =\u00a0m<\/em>dry<\/sub>.\u00a0 When I had not yet written this condition, the rocket’s mass would decrease to negative infinity.<\/p>\n The test rocket used in the simulations was the Falcon 1e, a theoretical rocket designed by SpaceX in 2011 as an update to the earlier Falcon 1.\u00a0 The Faclon 1e was never launched, which gave me the inspiration to test its theoretical viability in my simulation.\u00a0 The Falcon 1e specifications, compared to other rockets which I tried in the simulation, are given below.<\/p>\n The next two figures show the rising altitude and velocity of the rocket, as well as the relative strength of each component force over time.<\/p>\n We can see that the Falcon 1e will fly for about 169 seconds before it runs out of fuel and the thrust force drops to zero. By that time it has reached an altitude of about 125 km, still short of the 200 km required for a comfortable low-earth orbit.\u00a0 This is consistent with our expectations: no single-stage rocket in history has every made it to orbit from Earth.\u00a0 The topic of a single-stage-to-orbit spacecraft (SSTO) is an area of current research, mostly among defense contractors.<\/p>\n The flatline in the velocity curve around 65-95 seconds represents the rocket approaching the sound barrier.\u00a0 The corresponding sudden increase in drag represents effects of compressibility and wave drag taking hold.\u00a0 Note how the graph of drag directly across the sound barrier around 95 seconds seems to exhibit more erratic, possibly chaotic fluctuation.\u00a0 This occurs because thrust and drag at this point in time are competing for dominance over the rocket’s net force.\u00a0 Once the sound barrier is breached, drag quickly becomes negligible, as the air waves move behind the rocket, and air density becomes very thin.\u00a0 Making the jump to supersonic speeds, therefore, is critical in order to leave the Earth’s atmosphere; although the jump requires a high amount of thrust and a sturdy spacecraft, choosing to remain at subsonic velocities would be less effective overall due to the higher drag on the rocket.<\/p>\n In order to get to space, our simulated rocket will need two stages.\u00a0 Specifications for both stages of the Falcon 1e are public on Wikipedia and the SpaceX website.\u00a0 In order to achieve this in the code, I added the following lines inside the for-loop:<\/p>\n In other words, when the rocket runs out of fuel, the program will continue for another 15 iterations of the loop, and then will break (nextstage was initialized as zero).\u00a0 After that, I reinitialized the specifications of the rocket for its second stage and then ran a second for-loop identical to the first, but without reinitializing the holding matrices for the variables.<\/p>\n The results for the two-stage Falcon 1e’s altitude trajectory, velocity, and force diagrams are shown below.\u00a0 The program was instructed to stop when the rocket reached an altitude of 250 km.<\/p>\n<\/a><\/p>\nII.1 Atmospheric Drag<\/h4>\n
<\/a><\/p>\n<\/a><\/p>\n<\/a><\/p>\n<\/a>\u00a0\u00a0<\/a><\/a><\/p>\nII.2 Wave Drag<\/h4>\n
<\/a><\/p>\n<\/a><\/p>\nIII. Gravitational Force<\/h4>\n
<\/a><\/p>\n<\/a><\/p>\nIV. Computational Method<\/h4>\n
\u2022 Initialize universal constants (G, M_earth, mol, R, P0, T0, rho0, g0)\r\n\u2022 Initialize rocket specifications (thrust, drag coef., diameter, area, wet and dry mass,\r\n specific impulse, mass flow rate)\r\n\u2022 Euler-Cromer method\r\n \u2023 initialize altitude (1 m), velocity (0 m\/s), mass (m0), air density (rho0), and\r\n gravity (g0)\r\n for<\/strong> t from 0 to tmax\r\n \u2023 update gravity, g(z)\r\n \u2023 update mass, m(t)\r\n \u2023 call density function, [rho,T,P] = f(z)\r\n \u2023 call wave drag function, Cd = f(v,T,C0)\r\n \u2023 calculate thrust force, FT(m)\r\n \u2023 calculate drag force, FD(rho,v,Cd,m)\r\n if<\/em> [drag vector condition]\r\n \u2023 Euler-Cromer: v = v + (thrust - drag - g)*dt\r\n z = z + v*dt\r\n \u2023 save matrix values for v, z, m, g, thrust, drag, and rho\r\n \u2023 save end time\r\n if <\/em>[crash condition]\r\n if <\/em>[empty condition]\r\n end\r\n<\/strong>\r\n\u2022 Plot z vs. t, v vs. t\r\n\u2022 Plot FT, FD, and g vs. t\r\n\u2022 Define density function, [rho,T,P] = f(z)\r\n\u2022 Define wave drag function, Cd = f(v,T,C0)<\/pre>\n
V. Results, Single-Stage Rocket<\/h4>\n
<\/a><\/p>\n<\/a><\/p>\n<\/a><\/p>\nVI. Results, Two-Stage Rocket<\/h4>\n
if z < 0 % rocket crashes or fails to launch\r\n break\r\n elseif m < empty % rocket runs out of fuel, mass becomes stable\r\n FT = 0;\r\n dm = 0;\r\n nextstage = nextstage + 1<\/strong>;\r\n end\r\n \r\n if nextstage == 15<\/strong>\r\n nextstage = 0;<\/strong>\r\n break<\/strong>\r\n end<\/strong><\/pre>\n