UFDC Home  myUFDC Home  Help 



Full Text  
PAGE 1 1 A SIXDEGREEOFFREEDOM LAUNCH VE HICLE SIMULATOR FOR RANGE SAFETY ANALYSIS By SHARATH CHANDRA PRODDUTURI A THESIS PRESENTED TO THE GRADUATE SCHOOL OF THE UNIVERSITY OF FLOR IDA IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF MASTER OF SCIENCE UNIVERSITY OF FLORIDA 2007 PAGE 2 2 Sharath Chandra Prodduturi PAGE 3 3 To my parents. PAGE 4 4 ACKNOWLEDGMENTS I would like to express my sincere gratitude to my supervisory committee chair (Dr. Norman G. FitzCoy) for his con tinuous guidance, support, and hel p. I am really thankful to him. I would also like to express my grat itude to my supervisory committee members (Dr. Warren E. Dixon and Dr. Gloria J. Wien s) for their suppor t and guidance. I would like to express my grat itude to my parents for all th eir moral and financial support, without which this task could not have been acc omplished. I would be nowhere without them. I would like to acknowledge my sisters (Shiri sha and Swetha) for their help and support throughout my life. I would like to thank my friends and colle agues from AMAS (Frederick Leve, Shawn Allgeier, Sharan Asundi, Takashi Hiramatsu, Ja ime Jos Bestard, Andrew Tatsch, Andrew Waldrum, AiAi Cojuangco, Dante Buckley, Ni ck Martinson, Josue Mu noz, Jessica Bronson and Gustavo Roman) for their advice, help and support. PAGE 5 5 TABLE OF CONTENTS page ACKNOWLEDGMENTS...............................................................................................................4 LIST OF FIGURES................................................................................................................ .........7 ABSTRACT....................................................................................................................... ..............9 CHAPTER 1 INTRODUCTION AND BACKGROUND...........................................................................11 2 EQUATIONS OF MOTION FORMULATION....................................................................19 Coordinate Frames.............................................................................................................. ....19 Kinematic Equation of Motion...............................................................................................24 Dynamical Equations............................................................................................................ ..27 Generalized External Forces...................................................................................................30 External Forces................................................................................................................30 Thrust force..............................................................................................................30 Aerodynamic forces (drag and lift)..........................................................................32 Gravitational force....................................................................................................33 External Moments...........................................................................................................34 Aerodynamic moments............................................................................................34 Gravitational moment...............................................................................................35 Thrust moment.........................................................................................................36 3 DESCRIPTION OF MODELS USED...................................................................................38 Gravity Model.................................................................................................................. .......38 Inertia Model.................................................................................................................. ........49 Strapon booster...............................................................................................................50 Cylindrical segment..................................................................................................50 Parabolic nose cone..................................................................................................52 Fins...........................................................................................................................54 Liquid Engine..................................................................................................................57 Solid Motor.................................................................................................................... ..59 Payload........................................................................................................................ ....61 Drag Coefficient Model......................................................................................................... .63 Center of Pressure Model.......................................................................................................64 Nose........................................................................................................................... ......66 Cylindrical Body.............................................................................................................67 Conical Shoulder.............................................................................................................67 Conical Boattail...............................................................................................................68 Fins (Tail Section)...........................................................................................................68 The WGS84 Ellipsoid Model.................................................................................................69 PAGE 6 6 4 SIMULATION RESULTS AND DISC USSION...................................................................73 Simulation..................................................................................................................... ..........73 Validation..................................................................................................................... ..........87 5 CONCLUSION AND FUTURE WORK...............................................................................91 Conclusions.................................................................................................................... .........91 Future work.................................................................................................................... .........92 APPENDIX A MATLAB FUNCTIONS AND SCRIPT................................................................................93 B SIMULATION CONFIGURATION....................................................................................116 LIST OF REFERENCES.............................................................................................................132 BIOGRAPHICAL SKETCH.......................................................................................................135 PAGE 7 7 LIST OF FIGURES Figure page 11 Spacebased range and range safety, today and future......................................................13 21 Relative orientation of the various frames.........................................................................21 22 Euler angles and the relative orientati on between the vehicle frame and the vehiclecentered horizontal frame..................................................................................................24 23 Geometry of the launch vehi cle and various position vectors...........................................28 24 External forces acting on a launch vehicle during its flight...............................................31 31 Representation of a position vector in Cartesian and Spherical coordinates.....................42 32 Cylindrical segment of the strapon booster......................................................................51 33 Parabolic nose cone........................................................................................................ ....53 34 Fin........................................................................................................................ ..............55 35 Liquid engine.............................................................................................................. .......58 36 Solid motor................................................................................................................ .........60 37 Payload.................................................................................................................... ...........62 38 Conical shoulder........................................................................................................... .....67 39 Conical Boattail........................................................................................................... ......68 310 Fin and Tail section...................................................................................................... ......69 311 Geodetic Ellipsoid and Geodetic coor dinates of an arbitrary point P............................70 41 Various parameters of the launc h vehicle as a function of time........................................77 42 Velocity of the launch vehi cle in the inertial frame...........................................................79 43 Position of the launch vehi cle in the inertial frame...........................................................80 44 Launch vehicle during the time of launch as seen from the J2000 inertial frame.............80 45 Moments of inertia of the launch vehicl e about its instantane ous center of mass.............81 46 Moments of inertia of the strapon booste r about the instantaneous center of mass of the launch vehicle and about its instantaneous center of mass..........................................82 PAGE 8 8 47 Moment of inertia of the first stage ab out the instantaneous center of mass of the launch vehicle and about its in stantaneous center of mass................................................84 48 Moment of inertia of the second stage a bout the instantaneous center of mass of the launch vehicle and about its in stantaneous center of mass................................................85 49 Moment of inertia of the third stage a bout the instantaneous center of mass of the launch vehicle and about its in stantaneous center of mass................................................86 410 The need for instrumental data or thrust vector in the vehicle frame................................90 B1 The DELTA II Launch vehicle geometry........................................................................116 B2 Strapon booster geometry...............................................................................................116 B3 Elements of DELTA II Launch vehicle and Strapon Booster........................................120 B4 Cylindrical shell.......................................................................................................... .....121 B5 Propellant shell........................................................................................................... ......122 B6 Parabolic nose cone........................................................................................................ ..123 B7 Fins....................................................................................................................... ............124 B8 First stage................................................................................................................ .........125 B9 Second stage............................................................................................................... ......127 B10 Third stage............................................................................................................... ........128 B11 Payload................................................................................................................... ..........130 B12 Strapon boosters around the Rocket...............................................................................130 PAGE 9 9 Abstract of Thesis Presen ted to the Graduate School of the University of Florida in Partial Fulfillment of the Requirements for the Degree of Master of Science A SIXDEGREEOFFREEDOM L AUNCH VEHICLE SIMULATOR, FOR RANGE SAFETY ANALYSIS By Sharath Chandra Prodduturi August 2007 Chair: Norman G. FitzCoy Major: Mechanical Engineering Failure of a launch vehicle during its launch or flight might pose a hazard to the general public. The United States Air Force Space Co mmand (USAFSC) operates the United States launch facilities and ensures safety to the gene ral public, launch area an d personnel, and foreign land masses in case of such a failure. To ensu re safety, USAFSC currently uses extensive groundbased systems, which are expensive to ma intain and operate and are limited to the geographical area. To overcome these draw backs, NASA proposed a concept called SpaceBased Telemetry and Range Safety (STARS) which uses spacebased assets to ensure safety. The STARS concept requires support tools in the form of simula tion softwares that provide the ability to quickly analyze new (or changes in) concept and ideas, an option not easily accomplished with hardware only. Trajectory and link margin analysis tool is one of these crucial support tools required by STARS. My study focused on modeling the full dynamics of a launch vehicle and development of a MATLAB based sixdegreeoffreedom simulato r for generating nominal and offnominal trajectories as part of the trajec tory and link margin analysis. In my study, the J2000 coordinate frame and the vehiclecentered horizontal frame were used as the reference frames to define the position and orientation of a la unch vehicle, respectively. Or ientation and the kinematic PAGE 10 10 equation of a launch vehicle are ex pressed in terms of quaternions instead of Euler angles, to avoid intensive computations and singularities. The equations of motions of a launch vehicle are developed by accounting for the variability in its mass and geometry. Various models are developed for calculation of quantities such as gravity, inertia, center of pressure and drag coefficient required for solving the equations of motion. The developed gravity model uses the spherical harmonic representation of the gravitati onal potential to account for the variability in Earths mass distribution and uses EGM96 (360 X 360) spherical harmonic coefficients and WGS84 Earth ellipsoid model. Th e gravity model is singularityfre e and numerically efficient. A novel way of calculating the variable mass/in ertial properties of a launch vehicle was developed. This inertia model is a simple and approximate model and considers general geometries to develop the inertia characteristics of a launch vehicle. The drag coefficient model from the Missile Datcom database is used in this research. The ki nematic equations, dynamic equations, gravity, inertia, center of pressure, dr ag coefficient and other models are implemented in MATLAB to form a sixdegreeoffreedom launch vehicle simulator. The results and discussions of a simulation perfor med using the developed simulator are presented in this thesis. A validation of the developed simulator was atte mpted with flight data available from NASA Kennedy Space Center; however, critical data need ed for the validation could not be provided due to ITAR restrictions. PAGE 11 11 CHAPTER 1 INTRODUCTION AND BACKGROUND Ensuring safe, reliable and affordable access to space is the fundamental goal of the U.S. range safety program [23]. The Public Law 60 established the national range system based on two primary concerns/factors: location and public sa fety. Thus, Range Safety, in the context of national range activities, is rooted in PL 60 [14]. Range is defined to be the volume through whic h the launch vehicle must pass in order to reach its destination from the la unch point, and its projection on ea rth (in case of a space vehicle, the destination can be outer space or a location on earth) [26]. A range includes space, facilities, equipment and systems necessary for testing and monitoring launches, landing and recovery operations of launch vehicles and other technica l and scientific program s and activities [18]. The United States launch facili ties are divided into Eastern and Western Ranges. The Air Force Space Command operates the launch f acilities of the United States. The 30th and the 45th space wings manage and operate the Western Ra nge and Eastern Range respectively. The Eastern Range comprises of Cape Canaveral Air Station and its owned or leased facilities and encompasses the Atlantic Ocean, including all surrounding land, sea, and air space within the reach of any launch vehicle extending eastward into the Indian and Pacific Oceans. The Western Range comprises of Vandenberg Air Force Base (VAFB) and its owned or leased facilities and encompasses the Pacific Ocean, including all su rrounding land, sea, and air space within the reach of any launch vehicle extending westwa rd through the Pacific and Indian Oceans[14]. The Eastern and Western Ranges, using a Range Safety Program provide safety to the public by ensuring that the risk to the general public from launch and flight of launch vehicles and payloads is no greater than that imposed by the over flight of conventional aircraft. Apart from public protection, the national range safety includes launch area safety, launch complex PAGE 12 12 safety, and the protection of national resources [ 14]. The objective of the Range Safety Program as stated in Eastern and Western Range 1271, Range Safety Requirements [14] is The objective of the Range Safety Program is to ensu re that the general public, launch area personnel, foreign land masses, and launch area resources ar e provided an acceptable level of safety and that all aspects of prelaunch a nd launch operations adhere to pub lic laws and national needs. The mutual goal of the Ranges and Range Users sh all be to launch launch vehicles and payloads safely and effectively with commit ment to public safety (14, p. 15). Range safety personnel evaluate vehicle de sign, manufacture and in stallation prior to launch and monitor vehicle and environmental conditions during countdown. Range safety personnel also monitor the performa nce of launch vehicles in flight and are responsible for their remote destruction/termination if it should be judged that they pose a hazard. For all vehicle termination cases, propulsion is terminated and based on the vehicle type, stage of flight, and other circumstances of failure, the method of te rmination might vary. Depending on factors like geographic location and population, the vehicle may be destroyed to disperse the propellants before surface impact, or it may be kept intact to minimize the debris footprint. The launch vehicle is also equipped with a br eakwire or lanyard pull to initiate a flight termination in case of a premature stage separation [23]. Extensive groundbased systems are utilized by the current United States Eastern and Western Ranges for realtime tracking, communicat ions, and command and control of the launch vehicles. These groundbased assets are very ex pensive to maintain and operate and are limited to the geographical area [31]. Therefore the current range sy stems need to be upgraded or replaced. According to Whiteman et al. [31], NASA Dryden Flight Res earch Center, Future spaceports will require new technologies to pr ovide greater launch and landing opportunities, PAGE 13 13 support simultaneous missions, and offer enhan ced decision support models and simulation capabilities. These ranges must also have lowe r costs and reduced complexity, while continuing to provide unsurpassed safety to the public, flig ht crew, personnel, vehi cles, and facilities. Commercial and government spacebased assets for tracking and commu nications offer many attractive possibilities to help achieve thes e goals (31, p. 2). Figure 11 shows the current primary Easter n and Western Ranges instrumentation sites (solid lines) and a possible futu re spacebased configuration with fewer groundbased assets (dashed lines). From Fig 11, it should be noted that the future spaceba sed configuration might still include some launchhead gr oundbased assets for visibility and rapid response times shortly after liftoff [31]. Figure 11. Spacebased range and range safety, today and future Reprinted with permission from D. E. Whiteman, L. M. Valencia, and J. C. Simpson, SpaceBased Range Safety and Future Space Range Appli cations, NASA Dryden Flight Research Center, Edwards, Californ ia. Rep. H2616, NASA TM2005213662, 2005. PAGE 14 14 SpaceBased Telemetry and Range Safety (STARS) SpaceBased Telemetry and Range Safety (S TARS) is a multifaceted and multicenter project to determine the feasibility of using sp acebased assets, includi ng the Tracking and Data Relay Satellite System (TDRSS) and Global Positioning System (GPS), to reduce operational costs and increase reliability. The STARS study was established by the National Aeronautics and Space Administration (NASA) to demonstrate th e capability of spacebased assets to provide communications for Range Safety (lowrate, ultrahi gh reliability metric tracking data, and flight termination commands) and Range User (video, voi ce, and vehicle telemetry) [31]. To support the envisioned future space range, new and improved systems with Range Safety and Ranger User capabilities are under tes ting and development. A brief description of the planned and completed phases of the STARS project is given below [31], [30], [10], [21]. Phase 1 Developed and tested a new Sband Range Safety system. During JuneJuly 2003, seven test flights were performed on a F15B aircraft at Dryden Flight Research Center using a Range User system representative of those on the current launch vehicles. Successfully demonstrated the basic ability of the STARS to establish and maintain satellite links with TDRSS and the GPS. Phase 2 The objective is to increase the Range User data rates by an order of magnitude by enhancing the Sband Range Safety system a nd a new telemetry system which utilizes a Kuband phasedarray antenna. TDRSS is the spacebased communication link (i.e., TDRSS provides th e tracking and data acquisition services between the launch vehicle/low earth orbiting spacecraft and NASA/customer control and data processing facilities [22]). PAGE 15 15 Phase 3 Phase 3 uses a small, lightweight hardware co mpatible with a fully operational system and demonstrates the ability to maintain a Kaband TDRSS communications link during a hypersonic flight. Develop smaller, lighter version of the Range Safety Unit for the Range Safety system in fiscal year 2006. TDRSS is the spacebased communication link. Test flights are planned fo r late fiscal year 2007. Space Based Range Safety system will be complete by the completion of Phase 3 development. Phase 4 Develop Kaband transmitter (NASA) and phased array antenna (AFRL) for Range User system in fiscal year 20062007. Perform flight test on aircraft (Flight Demo 3a) to test pe rformance of Glenn Research Centers (GRC) Kaband active phased array antenna in fiscal year 2007. Perform flight test of Kaband sy stem on F15B in fiscal year 2008. Refly phase 3 Range Safety Unit design with enhancements. Certification Phase Perform Certification of Range Safety and Range Users systems in fiscal year 2009. The STARS program was renamed to Space Base d Range Demonstration and Certification (SBRDC) program [20]. From the available information on the World Wide Web/internet, Phases 1, 2 and 3 are completed and the curren t status of the STARS/ SBRDC program is as stated in Phase 4 above [19]. The STARS concept requires support tools in the form of simulation softwares which provide the ability to quickly analyze new (or changes in) concepts a nd ideas, an option not easily accomplished with hardware only. Trajectory and link margin analysis tool is one of these crucial support tools required by STARS. The trajectory portion of the trajectory and link PAGE 16 16 margin analysis involves genera ting trajectories (and orie ntation) of a launch vehicle. The link margin portion involves calcula ting the telemetry link margin during the flight of a launch vehicle. Link margin is defined as the differe nce in dB, between the ma gnitude of the received signal at the receiver input and the receiver sensitivity (i.e., the minimum level of signal required for reliable operation). The higher the link marg in, the more reliable the communications link [24]. The trajectory and orientation of the launch vehicle calculated using the trajectory portion and dynamic parameters such as vehicl e antenna patterns, locations of ground stations and others are taken into account in order to compute the link margin. Trajectory and link margin analysis is frequently required to ensure adequate link closure for range safety [15], [8]. Trajectory and link margin analysis involve s simulating the launch vehicle for various failure scenarios and checking if the command upli nk can be closed with sufficient margin under the worst possible conditions and from any intend ed ground site(s). The worst possible failure scenario includes trajectories that result due to to tal loss of control of the launch vehicle. These trajectories might include a sudden heading change to a populated area or may consist of series of tumbles [8]. The Space Systems Group (University of Flor ida) and UCF collaborated to develop a MATLAB based tool for trajectory and link marg in analysis. The Space Systems Group is responsible for modeling the dynamics of a laun ch vehicle while UCF is responsible for the communications link model. The thrust of this thesis is to develop a MATLAB based launch vehicle model/simulator which is capable of simulating a launch vehicle in flight. Nominal launch vehicle trajecto ry simulation models have been done by many researchers. Researchers have also attempted to simulate th e offnominal trajectory of launch vehicles; e.g., Chen et al. [8] has suggested a threestep algor ithm to estimate the deviation of the launch PAGE 17 17 vehicle from the nominal trajectory. In their ap proach, the sudden accelerations are treated as the artificial maneuver controls focusing on kinematics instead of dynamics. This research intends to model the full dynamics of the launch vehicle. In Chapter 2, the equations of motions of the launch vehicle are developed. The definitions of the various coordinate frames used in the development of the simulator and the transformations between them are discussed in detail. Following the above, the development of the kinematic and dynamic equations of motion of the launch vehicle is presented. Finally, the various external forces and exte rnal moments (acting on the launch vehicle) to be included in the external force and moment terms in the equations of motion of the launch vehicle are discussed. In Chapter 3, the models used in the deve lopment of the simulator are presented. The development of the gravity, inertia, drag coefficient, center of pressure and WGS84 ellipsoid models is presented in detail in this chapter. The gravity model computes the acceleration due to gravity of Earth at a point of interest using its ECEF coordinates; the inertia model computes the mass properties of a launch vehicle as a function of time; the drag coefficient model computes the coefficient of drag for the launch vehicle as a function of position and ve locity; the center of pressure model computes the center of pressure for a specific geometry of the launch vehicle; the WGS84 ellipsoid model defines a reference Earth ellipsoid and is used to compute the altitude of a point of interest using its ECEF coordinates. In Chapter 4, the simulation results are pres ented and discussed. Simulation results of a DELTA II launch vehicle model for a fictitious thrust profile (constant axial thrust) are discussed. Following the above, the details of th e attempt to validate the developed simulator with flight data available from NASA Kennedy Sp ace Center were presented. It is shown that PAGE 18 18 the simulator cannot be validated due to the lack of availabi lity of critical data (an ITAR1 issue). Finally, in Chapter 5, the conclusions of this res earch and the possible future work are discussed. 1 ITAR International Traffic in Arms Regulations PAGE 19 19 CHAPTER 2 EQUATIONS OF MOTION FORMULATION This chapter discusses the equations of moti ons (i.e., the dynamic and kinematic equations) of a launch vehicle. First the background is presented and then the derivations of the equations of motions of an expendable laun ch vehicle are presented. Finally the generalized forces acting on a launch vehicle during its flight are discussed. The following assumptions are made in this research [9]. The launch vehicle (with the strapon boosters) is assumed to be rigid. The center of mass of the launch ve hicle lies on the longitudinal axis. The longitudinal axis is the principal axis of inertia. Coordinate Frames In order to derive the equations of motion of a launch vehicle that describe its position and orientation as a function of time, various coordi nate frames are considered. These frames are discussed below. Inertial frame (XiYiZi): For studying the launch vehicle motion in the vicinity of Earth and at an interplanetary level, the J2000 frame is consider ed as an inertial frame. This frame has the origin at the Earths center of ma ss; its positive Zaxis points towards the Earths North Pole and coincides with the Earths rotational axis. The pos itive Xaxis lies in the Earths equatorial plane and points towards the vernal equi nox in J2000 epoch. The Yaxis lies in the equatorial plane and completes a righthanded Cartesian frame [9], [28]. Rotating geocentric frame (XgYgZg): This frame rotates with the rotating Earth. This frame has its positive Zaxis pointed towards the Eart hs North Pole and coincides with the Earths rotational axis. The posi tive Xaxis lies in the equatorial pl ane, crossing the upper branch of the PAGE 20 20 Greenwich meridian. The Yaxis lies in the equatorial plane and co mpletes a righthanded Cartesian frame [9]. Vehiclecentered horizontal frame (XvYvZv): The orientation of the launch vehicle and its velocity vector relative to the Ea rths surface can be described using this frame. The origin of this frame coincides with the initi al center of mass of the launch ve hicle. The orientation of the frame remains fixed through out the flight of th e launch vehicle. The XY plane of this frame coincides with the initial local ho rizontal plane (the local horizont al plane is the plane normal to the radius vector from the center of mass of the Earth to the center of mass of the launch vehicle). The positive Xaxis po ints north and lies along the nort hsouth direction. The positive Yaxis points east and lies along the eastwest direction. The Zvaxis is along the radius vector from the center of Earth and is positive downwards [9]. Vehicle frame (XrYrZr): The origin of this frame coincides with the initial center of mass of the launch vehicle. The Xaxis coincides with the longitudinal axis of the launch vehicle and is positive forwards (i.e., towards the nose of the la unch vehicle). The Yaxis and Zaxis lie along the other two principal axes of inertia of the vehicle such that they complete a righthanded Cartesian frame [9]. Relative Orientations Figure 21 shows the relative orientations of th e various frames. The details of the relative orientations and transformations between th e above described frames are given below. Inertial frame/Rotating geocentric frame [9] : The Earth and therefore the rotating geocentric frame rotate about the Zaxis of the inertial frame with an angular velocity of Earth (e ). Thus, the relative orientation of these frames is determined by a rotation about the Zaxis PAGE 21 21 Figure 21. Relative orientat ion of the various frames PAGE 22 22 through an angle that is equa l to the angle between the Xiaxis and Xgaxis. This angle is equal to the Greenwich hour angle of the vernal equinox HG. If both frames coincide at t = t0, the angle HG at any time is given in Eq. 21. 0 e GHtt (21) Since the inertial frame in our case is the J2000 frame, the term 0tt is equal to the time elapsed in seconds from January 1, 2000, 12: 00 UTC until the time t of interest. The transformation between the frames is given in Eq. 22. The vectors GE and IE in Eq. 23 represent an arbitrary vector E coordinatized in the rotating ge ocentric frame and the inertial frame respectively. The transformation matrix is given in Eq. 23. GI GIECE (22) cossin0 sincos0 001GG GIGGHH CHH (23) Rotating geocentric frame/Vehiclecentered horizontal frame [9] : The relative orientation of these two frames can be determined by mean s of two successive rotations. The rotating geocentric frame (XgYgZg frame) is first rotated about its Zaxis (i.e., Zg axis) by an angle the geographic longitude of the launch vehicle. This new frame is then rotated about its new Yaxis by an angle 2 where is the geocentric latitude of the launch vehicle. The resulting frame has the same orientation as the vehiclecentere d horizontal frame. The transformation between the frames is given in Eq. 24. The vectors VE and GE in Eq. 24 represent an arbitrary vector E coordinatized in the vehiclecentere d horizontal frame and the rotating geocentric frame respectively. The tran sformation matrix is given in Eq. 25. VG VGECE (24) PAGE 23 23 sincossinsincos sincos0 coscoscossinsinVGC (25) Vehiclecentered horizontal frame/Vehicle frame [9] : The relative orientation of these two frames can be determined by means three successi ve rotations as shown in Fig. 22. The three angles through which these three successive rota tions are performed are called Euler angles. The vehiclecentered horizontal frame is fi rst rotated about its Zaxis (i.e., v Z axis) by an angle to obtain a new frame 111vvvXYZ. is called the yaw angle, the a ngles between the vertical plane through the longitudinal axis of the launch vehicle and the vXaxis. Then the new frame 111vvvXYZ is rotated about its Yaxis (i.e., 1vY axis) by an angle to obtain another new frame 222vvvXYZ. is called the pitch angle, the angle be tween the longitudinal axis of the launch vehicle and the local horizontal plan e. Finally, the newest frame, 222vvvXYZ, is rotated about its Xaxis (i.e., 2vXaxis) by an angle to obtain the vehicle frame rrrXYZ. is called the bank angle, the angle between the r Z axis and the vertical plan e through the longitudinal axis of the launch vehicle. The tran sformation between the frames is given in Eq.26. The vectors RE and VE in Eq. 26 represent an arbitrary vector E coordinatized in the vehicle frame and the vehiclecentered horizontal fr ame respectively. The transforma tion matrix is given in Eq. 27. In Eq. 27, C and S are used to represent the cosine and sine of an angle RV RVECE (26) RVCCCSS CCSSSCCCSSSSC SSCSCSCCSSCC (27) PAGE 24 24 Figure 22. Euler angles and the relative orient ation between the vehicle frame and the vehiclecentered horizontal frame Inertial frame/Vehicle frame [9] : The transformation from the in ertial frame to the vehicle frame can be obtained by successively applying the transformations GIC, VGC and R VC to the inertial frame. The transformation between th e frames is given in Eq. 28. The vectors RE and IE in Eq. 28 represent an arbitrary vector E coordinatized in the vehi cle frame and the inertial frame respectively. The transformation matrix is given in Eq. 29. R I RIECE (28) RIRV VGGICCCC (29) Kinematic Equation of Motion The rotational kinematic equation of motion rela tes the orientation and the angular velocity of a launch vehicle. The derivation of th e kinematic equation is presented below. Let 1 2 3 be the angular velocity of the vehi cle frame with respect to the vehiclecentered horizontal frame expressed in the vehicle frame. Since the vehiclecentered horizontal frame is an inertial frame, is the absolute angular velocity of the launch vehicle. Let PAGE 25 25 and be the Euler angle rates for the 321 Eule r rotation sequence from the vehiclecentered horizontal frame to the vehicle frame. The angular velocity of the launch vehicle can be expressed in terms of the Euler rates as given in Eq. 210. The rotation matrices 1 RVC andRVC in Eq. 210 are given in Eqs. 211 and 212. 11 2 30 0 00 00RVRVCC (210) 1100cos0sin 0cossin010 0sincossin0cosRVC (211) 100cos0sincossin0 0cossin010sincos0 0sincossin0cos001RVC (212) Equation 210 can be rewritten as Eq. 213 where the matrix X in Eq. 213 is given in Eq. 214. Equation 213 can be rewritten as Eq. 215. The matrix X in Eq. 214 is inverted and substituted into Eq. 215 to obtain Eq. 216. 1 2 3X (213) 21 21 211,11,21,3 2,12,22,3 3,13,23,3RV RVRV RV RVRV RV RVRVXCCC CCC CCC (214) 11 2 3X (215) 1 2 3cossinsincossin 1 0coscossincos cos 0sincos (216) PAGE 26 26 Eq. 216 is the kinematic equation of motion of the launch vehicle. This Euler angle representation of the relative orientation of the vehiclecentered horizontal frame and the vehicle frame has the following disadva ntages (i) singularity at 2 and (ii) solving the kinematic equation of motion Eq. 216 is com putationally intensive as it invo lves trigonometric quantities. To avoid these problems, quaternions are used to represent the relative orie ntation of the vehiclecentered horizontal frame and the vehicl e frame. The transformation matrix RVC can also be expressed in terms of quaternions as shown in Eq. 217. The quantites 0123and,, qqqq in Eq. 217 are calculated using the e xpressions in Eqs. 21821. 22 0112031302 22 1203022301 22 12032301032212222 2222122 2222221RVqqqqqqqqqq qqqqqqqqqq qqqqqqqqqqC (217) 0coscoscossinsinsin 222222q (218) 1coscossinsinsincos 222222q (219) 2cossincossincossin 222222q (220) 3sincoscoscossinsin 222222q (221) The kinematic equation of motion in terms of quaternion rates is give n in Eq. 222. The quantites 123and, in Eq. 222 are the components of the angular velocity vector, of the vehicle frame with respect to the vehiclecentered horizontal frame expressed in the vehicle frame 1 2 3i.e., Since the vehiclecentered horizontal frame is an inertial frame, is the absolute angular veloc ity of the launch vehicle. PAGE 27 27 123 132 231 32100 11 22 330 0 1 0 2 0qq qq qq qq (222) Dynamical Equations In this section, the derivation of the dynamic equations of a launch vehicle is presented. Figure 23 shows the geometry of the launch vehi cle and the various position vectors considered in the derivation of the dynamical equations. The instantaneous center of mass of the launch vehicle is represented by C a nd the initial center of mass of th e launch vehicle is represented by C1. The position vector of the mass element dm relative to the origin of the inertial frame (coordinatized in the inertial frame) is represented by IdmR. The position vector of the initial center of mass of the launch vehicle C1 relative to the origin of the inertial frame (coordinatized in the inertial frame) is represented by 1I CR. The position vector of the instantaneous center of ma ss of the launch vehicle C with re spect to the initial center of mass of the launch vehicle C1 (coordinatized in the vehi cle frame) is represented by Rcr. The position vector of the mass element dm with resp ect to the instantaneous center of mass of the launch vehicle C (coordinatized in the vehicle frame) is represented by R r. From Fig. 23, the position vector of the mass element dm can be written as expressed in Eq. 223. The matrix R IC in Eq. 223 is the transformation matrix from the inertial frame If to the vehicle frame Rf obtained from Eq. 29. The acceleration of the mass elemen t dm is obtained by differentiating Eq 223 twice with respect to tim e t as shown in Eq. 224. Th e resulting expression for the acceleration of the mass element dm is given in Eq. 225. PAGE 28 28 Figure 23. Geometry of the launch vehicle and various position vectors 1IITRR RI C dmcrRRCr (223) 12 2IITRR RI C dmcrd RRCr dt (224) 12RRRRR c IIT RI C dm RRRRRRRc ccr rrrrr RRC rr (225) Applying Newtons second law, we obt ain Eq. 226. The expression for I dm R from Eq. 225 is substituted into Eq. 226 and then inte grated over the entire launch vehicle mass as shown in Eq. 227. The resulting expressi on after integration is given in Eq. 228. I dmdFRdm (226) 12RRRRR IextIIT RI C dm RRRRRRR MMcc ccrr rrrr F RdmRCdm rr (227) PAGE 29 29 12IextITRRRRRRRR RI CccccrrrrFMRCMMMM (228) ( 0and0Mrrrdm for a rigid body) Equation 228 is the translati onal equation of motion of a la unch vehicle. The term ext F in Eq. 228 represents the result ant of the external forces ac ting on the launch vehicle. The external forces acting on the laun ch vehicle during its flight are described in the next section. The term cr and its time derivatives in Eq. 228 ar e obtained by computing the instantaneous center of mass of the launch vehi cle with respect to the initia l center of mass of the launch vehicle as a function of time using the mass prope rties of the launch vehicle, the mass/inertia properties of a launch vehicle are discussed in the inertia model in Chapter 3. A brief description of the procedure used to compute cr and its time derivatives in the simulator is presented in Chapter 4. Taking moments of all the forces ab out the instantaneous center of mass C, we obtain Eq. 229. Substituting the expression for I dm R from Eq. 225 into Eq. 226 and then substituting the resultant expression for dF into Eq. 229, and then integrating the resultant expression over the entire launch vehicle mass, we obtain Eq. 230. cdMrdF (229) RextRRRRRRR c MM M rrdmrrdm (230) The terms andRRRRRRR MMrrdmrrdm in Eq. 230 are represented in terms of moment of inertia tensor of the launch vehicle, I as shown in Eqs. 231 and 232. Substituting Eqs. (231) and (232) into Eq.230, we obtain Eq. 233. RRRR M Irrdm (231) PAGE 30 30 RRRRRR M Irrdm (232) RRRRext cIIM (233) Equation 233 is the rotationa l equation of motion of a launch vehicle. The term ext cM in Eq. 233 represents the result ant moment (about the instantaneous center of mass of the launch vehicle) of the external forces acting on the launch vehicle. The moments of the external forces acting on the launch vehicle du ring its flight are descri bed in the next section. Generalized External Forces To solve the translational and rotational equa tions of motions of th e launch vehicle given by Eqs. 228 and 233, the external forces and th e moments (of the external forces) acting on the launch vehicle need to be calcula ted. This section discusses th e external forces and moments (due to external forces) acting on a launch vehicle. First, the external forces acting on a launch vehicle are discussed and then the moments due to the external forces acting on a launch vehicle are discussed. External Forces Figure 24. depicts the external forces acti ng on a launch vehicle during its flight. The external forces acting on the laun ch vehicle during its flight are (i) thrust force, (ii) aerodynamic forces (lift and drag) and (iii) gravity (weight) Thrust force The thrust force acting on a launch vehicle is eeaeTmVppA where eVis the exhaust speed of the gases rela tive to the launch vehicle; m is the propellant mass flow rate; ep is the pressure at nozzle exit; apis the ambient pressure; e A is the nozzle exit (exhaust) area. PAGE 31 31 Figure 24. External forces acting on a launch vehicle during its flight The thrust force of a launch vehicle is genera lly expressed in the vehicle frame as given in Eq. 234. The simulator requires the thrust force to be coordinatized in the inertial frame, this can be obtained by using the rotation matrix from th e vehicle frame to the inertial frame from Eq. 29, T IRRICC. The expression for the thrust force coordinatized in the inertial frame is given in Eq. 235 R Thrust x y zT F T T (234) ITR RI ThrustThrustFCF (235) For a launch vehicle composed of strapon boos ters, solid motors and liquid engines, the thrust acting on a launch vehicle at any instant is equal to the v ector sum of the thrusts provided by all of the strapon boosters, solid motors and liquid engines at that in stant. The simulator requires the user to i nput the thrust profile of a launch vehicle. PAGE 32 32 Aerodynamic forces (drag and lift) The aerodynamic forces acting on a launch vehi cle can be neglected at altitudes greater than or equal to 600 km [29]. However, th e aerodynamic forces acting on a launch vehicle below 600 km cannot be neglected and the expr essions for these forces are shown below. Drag The drag force acting on a launc h vehicle is expressed as 1 2DragD F CAvv where DC is the drag coefficient; A is the crosssectional area perpendicular to the flow; is the density of the medium, andvv are the speed and velocity of the launch vehicle relative to the medium. The simulator requires the drag force to be coordinatized in the inertial frame. The expression for the drag force coordinatized in the inertial frame is given in Eq. 236. 1 2 I I DragD F CAvv (236) For a launch vehicle composed of strapon booste rs and a main section (i.e., the section consisting of the different stages of the launch vehicle and payload), th e resultant drag force acting on a launch vehicle is equal to the vector sum of the drag forces acting on the strapon boosters and main section of the launch ve hicle. The density of the medium/air, depends on the altitude of the launch vehicl e. The altitude of the launch vehicle is computed using the WGS84 ellipsoid model discussed in Chapter 3. The procedure to calculate the drag coefficient is presented in the drag coefficient model in Chapter 3 Lift For low angles of attack, the lift force can be neglected. However, for high angles of attack, the lift force cannot be neglected. The lift force acting on a launch vehicle is expressed as 21 2L Lift F CAvv where L C is the lift coefficient; A is the surface area of the lifting PAGE 33 33 surface; is the density of the medium; v is the speed of the launch vehicle relative to the medium and v is a unit vector normal to th e velocity of the launch vehicle. It should be noted that the lift coefficient, L C, is a function of the angle of attack. The simulator requires the lift force to be coordinatized in the inertial frame. The expression for the lift force coordinatized in inertial frame is given in Eq. 237. 21 2I L LiftI F CAvv (237) For a launch vehicle composed of strapon booste rs and a main section (consisting of the different stages of the launch vehicle and payl oad), the resultant lift force acting on a launch vehicle is equal to the vector sum of lift forces acting on the st rapon boosters and main section of the launch vehicle. The density of the medium/air, depends on the altitude of the launch vehicle. The altitude of the launch vehicl e is computed using th e WGS84 ellipsoid model discussed in Chapter 3. In the current simula tor, the lift force acting on a launch vehicle is neglected. Gravitational force The gravitational force ac ting on a launch vehicle is WMg where M is the mass of the launch vehicle; g is the acceleration due to Earths gravit ational field. The gravitational force can be best expressed in the inertial frame, the gravitational force acts approximately in the negative direction along the radius vector from the center of the ea rth to the center of mass of the launch vehicle. The expression for the gravitationa l force coordinatized in the inertial frame is given in Eq. 238. I g x y zW F W W (238) PAGE 34 34 For a launch vehicle composed of strapon boosters, solid motors, liquid engines and payloads, the gravitational for ce acting on the launch vehicle is equal to the product of the instantaneous mass of the launch ve hicle (i.e., sum of instantaneous masses of all the elements of the launch vehicle) and the accele ration due to gravity vector ac ting at the instantaneous center of mass of the launch vehicle. The gravitat ional force is assumed to be acting at the instantaneous center of mass of th e launch vehicle. The procedur e to calculate the acceleration due to gravity vector is presente d in the gravity model in Chapte r 3. The instantaneous mass and the instantaneous center of mass of a launch vehi cle can be calculated using the mass properties of the launch vehicle, the mass/in ertia properties of a launch vehicl e are discussed in the inertia model in Chapter 3. A brief description of th e procedure used to compute the instantaneous mass and center of mass of a launch vehicle in the simulator is presented in Chapter 4. The resultant external force acting on a launc h vehicle is the vector sum of the all the forces acting on the launch vehicle as shown in Eq. 239. Substituting the external forces from Eqs. 23538 into Eq. 239, we obtain the expre ssion for the resultant external force acting on the launch vehicle given in Eq. 240 IextIIII g Drag ThrustLift F FFFF (239) 211 22IextTI RIDL x x I yy zzTW F CTCAvvCAvvW TW (240) External Moments The moments due to the external forces (i) thru st, (ii) drag, (iii) lift and (iv) gravity acting on a launch vehicle are discussed below. Aerodynamic moments Aerodynamic moments due to the separation of the center of pressure and center of mass are typically nonzero. The moments due to th e aerodynamic forces acting on a launch vehicle PAGE 35 35 can be neglected at altitudes gr eater than or equal to 600 km [29]. However, the aerodynamic moments acting on a launch vehi cle below 600 km cannot be negl ected and the expressions for these moments are shown below. The resultant aerodynamic for ces (i.e., lift and drag) acti ng on a launch vehicle are assumed to be acting at the center of pressure of the launch vehicle. The moments due to drag force and lift force about the instantaneous center of mass of the launch vehicle are given in Eqs. 241 and 242 respectively, where Pr is the position vector of th e center of pressure of the launch vehicle with respect to the instantaneous cent er of mass of the launch vehicle. The vector Pr can be computed by computing the center of pressure and cen ter of mass locations of a launch vehicle. The procedure to compute the center of pressure of a launch vehicle is presented in the center of pressure model in Chapter 3. The center of mass of the launch vehicle can be calculated using the mass propertie s of the launch vehicle, the ma ss/inertia properties of a launch vehicle are discussed in the inertia model in Chapter 3. DragDrag C pIIIMFr (241) L iftLift C pIII M Fr (242) In the current simulator, the lift force and its moment acting on the launch vehicle are neglected. Gravitational moment The gravitational force acts at the instantane ous center of mass of the launch vehicle. Therefore the gravitational moment about the same point (instantaneous center of mass of the launch vehicle) is a zero vector. The expressi on for the gravitational moment is given in Eq. 243. 0g CIM (243) PAGE 36 36 Thrust moment If thrust vectoring is consid ered, the moment due to thrust force of the launch vehicle about the instantaneous center of mass of the la unch vehicle will be nonzero and will be the major factor affecting the attitude of the vehicle. If thrust force is considered to act always along the longitudinal axis (i.e., no thrust vectoring) the moment due to thrust force of the launch vehicle about the instanta neous center of mass of the launch ve hicle will be zero. The expression for the moment due to thrust force of the laun ch vehicle about the inst antaneous center of mass of the launch vehicle is given in Eq. 244, where inr is the position vector of the center of nozzle of the ith thrusting element (from which the burnt fuel/gases exits from the launch vehicle) with respect to the instantaneous center of mass of the launch vehicle. The position vectors inr of all the thrusting elements can be calculated from the knowledge of the ge ometry of the launch vehicle and the location of the cen ter of mass of the launch vehicl e. The center of mass of a launch vehicle can be calculated using the mass pr operties of the launch ve hicle, the mass/inertia properties of a launch vehicle are di scussed in the inertia model in Chapter 3. A brief description of the procedure used to compute the position vectors inr in the simulator is presented in Chapter 4. ThrustThrust C ni iIIIiMFr (244) Therefore the resultant external moment acting on a launch vehicle is the vector sum of all the moments due to the external forces as shown in Eq. 245. Drag LiftThrust C C CIextIII CMMMM (245) Substituting Eq. 240 into Eq. 228, yields Eq. .246, which is the composite translational equation of motion for a launch vehicle. Similarl y, substituting Eqs. 29 and 245 into Eq. 233 yields Eq. 247, which is governing equation fo r the rotational motion of a launch vehicle. PAGE 37 37 2111 22 2TI RIDL ITRRRRRRRR cccc RI Cxx I yy zzTW CTCAvvCAvvW TW MRCMrMrMrMr (246) TRRR DragRI LiftThrust C C CIIICMMMII (247) This chapter presented the development of the kinematic equation of motion and the dynamic equations of motion of a launch vehicle. To solve the dynamic equations given in Eqs. 246 and 247, quantites such as gr avity, inertia, center of pressu re, drag coefficient and others need to be calculated. Chapter 3 presents so me of the models developed to calculate these quantities. PAGE 38 38 CHAPTER 3 DESCRIPTION OF MODELS USED In the previous chapter, the derivation of th e equations of motion and the discussion of the generalized forces were presented. In or der to implement these in the MATLAB based simulator, quantities such as inerti a, gravity, drag coefficient and al titude need to be calculated as time dependent functions. In this chapter the mo dels employed to calculate the above quantities in the simulator are described. The models desc ribed in this chapter are (i) gravity model (ii) inertia model (iii) drag coefficient model (iv) center of pressure model and (v) WGS84 ellipsoid model. The gravity model calculates the gravity v ector at a point of interest. The inertia model calculates the mass properties of a launch vehicle; the center of pressure model calculates the location of center of pressure necessary fo r computing the aerodynami c moments about the center of mass of a launch vehicl e; the drag coefficient model calculates the drag coefficient necessary for the computation of drag for ces and moments; the WGS84 ellipsoid model calculates the altitude of a launch vehicle. Th e altitude of a launch ve hicle is required to calculate the density of air and the Mach number of the launch vehi cle which in turn are required for the computation of the aerodynamic forces and th e drag coefficient respectively. The details of the models are given below. Gravity Model The gravity model presented below calculates the acceleration due to gravity acting at a point of interest due to the gravitational field of the Earth. The acceleration vector obtained from this model is coordinatized in the ECEF frame. The details of the model are presented below. The Earth is not a spherically symmetric mass body. The Earth bulges at the equator as a consequence of its rotation. The density of the Ea rth also varies from location to location. This variability of the Earths mass is modeled using a spherical harm onic expansion of the PAGE 39 39 gravitational potential [17]. The spherical harmon ic representation of the gravitational potential function (V) is given in Eq. 31 [27]. max201sincossinn n n m nnmnm nmGMa VPCmSm rr (31) In Eq 31, ()GM is the Earths gravitational constant; ,,r are the distance from the Earths center of mass to the point of interest, geocentric latitude and geocentric/geodetic long itude of the point of interest, respectively; ais the semimajor axis of the WGS84 ellipsoid (discussed later); nm are the degree and order of the spherical harmonic function respectively; ,nmnmCS are the spherical harmonic coefficients of degree n and order m; m nP is the associated Lege ndre function of degree n and order m. The associated Legendre function m nP is defined as follows [27]. sinm nP cossin sinm m n md P d sinnPLegendre polynomial 21 sin1 2! sinnn n nnd d A gravitational model is defined by the set of constants a and the spherical harmonic coefficients ,nmnmCS The gravitational model used in this research is the WGS84EGM96 model. The spherical harmonic representation of the gravitational potential function in Eq. 31 has numerical computation problems in the fo rm of the (unnormalized ) spherical harmonic coefficients, ,nmnmCS and the associated Legendre functions, sinm nP The (unnormalized) spherical harmonic coefficients, ,nmnmCS, tend to very small values and the associated Legendre functions, sinm nP tend to very large values as the degree increases. These problems can be circumvented by normalizi ng the associated Legendre function and the PAGE 40 40 spherical harmonic coefficients. In general, this normalization is achieved by multiplying the spherical harmonic coefficients (a nd dividing the associat ed Legendre function) by a scale factor which depends on the degree and order of the f unction. Denoting the normalized quantities by an overbar, the normalized spherical harmonic coefficients and the normalized associated Legendre polynomial are defined in Eqs. 32 and 33 respectively [27]. 1 2! !21nm nm nm nmC nm C S nmnk S for 0,1; 0,2 mk mk (32) 1 221 sinsin! !mm nnnmnk PP nm for 0,1; 0,2mk mk (33) In Eqs. 32 and 33, ,nmnmCS are the normalized spherical harmonic coefficients; m nP is the normalized associated Legendre function. The advantages of the normalization of the associated Legendre polynomial and spherical harm onic coefficients can be retained by rewriting the spherical harmonic representation of the gravitational potential function (V) in Eq. 31 in terms of normalized spherical harmonic coefficients and normalized associated Legendre function as given in Eq. 34 [27]. max201sincossinn n n m nnmnm nmGMa VPCmSm rr (34) The above normalized spherical harmonic repr esentation of the gravitational potential Eq. 34 is used to compute the acceleration due to gravity. The acceleration due to gravity can be computed by taking the gradie nt of the gravitational potential V as given in Eq. 35. The gradient of the potential V is shown in Eq. 36. It should be noted that the ,, x yz in Eq. 36 represents the Cartesian coor dinates in the ECEF frame. gV (35) PAGE 41 41 VVV Vijk x yz (36) Let ,, x yzggg be defined as xV g x yV g y and zV g z Since V is a function of ,,r ,and x yzggg must be expressed in terms of the partial derivatives of V with respect to ,and r as in Eqs. 379. xVVrVV g x rxxx (37) yVVrVV g y ryyy (38) zVVrVV g zrzzz (39) The expressions for the terms ,andVVV r in Eqs. 379 are given in Eqs. 310313, where the expression for the term sinm nP in Eq. 311 is given in Eq. 312. The term 1C in Eq. 311 is defined as 1 2 1!21 !nmnk C nm max 2201 1sincossinn m nnmnm nn n nmna VGM PCmSm rrr (310) max1 20sincossinn n n m nnmnm nmVGMa PCmSm rrC (311) 1sinsintansinmmm nnnPPPm (312) max20sinsincosn n n m nnmnm nmVGMa PmCmmSm rr (313) To calculate all other partial derivatives in Eq s. 37, 38 and 39, the spherical coordinates ,,r must be represented in term s of the Cartes ian coordinates ,, x yz. To do this, let us consider the geometry shown in Fig. 31. PAGE 42 42 Figure 31. Representation of a position vect or in Cartesian and Spherical coordinates From Fig. 31, the Cartesian position vector, ,, R xyz can be represented in terms of spherical coordinates, ,,r as given in Eqs. 314316. The partial derivatives of the spherical coordinates ,,r with respect to Cartesian coordinates ,, x yzare calculated using Eqs. 31416 and are expressed in Eqs. 317319. 222rxyz (314) 1tan y x (315) 1 22tan z xy (316) 22 222,,rxzxy x rxxxy xyr (317) 22 222,,ryzyx y ryyxy xyr (318) PAGE 43 43 22 2,0,xy rz zrzrz (319) Substituting Eqs. 31019 into Eqs. 37, 38 and 39, we obtain the X, Y and Z components of the acceleration vector in the ECEF frame. The computation of the acceleration vector using the above procedure becomes cu mbersome at or near poles (i.e., at 90 ) because Eq. 312 contains tan which becomes infinite when approaches 90. Furthermore, Eqs. 317 and 318 also become infinite as and/or x y approach zero. A change of coordinates is car ried out which avoids the above difficulty yet retains similar recursive and orthogonality prop erties. The coordinate chan ge and the derivation of the acceleration vector in term s of these new coordinates (as a re sult of the coordinate change) is presented below [25]. Let the new coordinates be and,,urst where each coordinate is defined as 1 222 2rxyz ,, x yz stu rrr(i.e., r is the scalar magnitude of the vector R ,,ust are the three components of the unit vector R ). From the above definition of the new coordinates, it can be seen that 2221 stu and s R rt u To obtain the acceleration vector g in terms of these new coordinates, the gravitational potential V in Eq. 34 should be first expressed in terms of these new coordinates. In place of the normalized associated Legendre polynomial .m nP, a new polynomial .m nA is defined as given in Eq. 320. In place of sinandcos mm in the gravitational potential V expression Eq. 34, define ,m R st and PAGE 44 44 ,m I st as given in Eqs. 32122. Using the definitions in Eqs. 32022, the gravitational potential V in Eq. 34 can be expressed in terms of these definitions as given in Eq. 323. 21 221 1 1 2! !!nm n m n nnmuunmnk d A nmndu (320) ,coscosmmRstmReal part of msit (321) ,sincosmmIstmImaginary part of msit (322) max201,,n n n m nnmmnmm nmuGMa VACRstSIst rr (323) The above spherical harmonic representation of the gravitational poten tial Eq. 323 is now expressed in terms of the new coordinates and,,urst and is used to compute the acceleration vector. The acceleration vector can be computed by taking the gradient of the gravitational potential V as given in Eq. 324. The gradient of the potential V is defined in Eq. 325. It should be noted that the ,, x yz in Eq. 325 represents the Cart esian coordinates in the ECEF frame. gV (324) VVV Vijk x yz (325) Let ,, x yzggg be defined as xV g x yV g y and zV g z Since, V is a function of ,,, rstu, ,and x yzggg must be expressed in terms of the partial derivatives of V with respect to ,,andrstu as given in Eqs. 32628. xVVrVsVtVu g x rxsxtxux (326) yVVrVsVtVu g y rysytyuy (327) zVVrVsVtVu g zrzsztzuz (328) PAGE 45 45 Calculating and substi tuting the quantities ,,,,,,rrrsss x yzxyz ,,,ttt x yz ,anduuu x yz in Eqs. 32628 and then substituting the resulting expressions of Eqs. 326 328 into Eqs. 324 and 325, we obtain the expression for the acceleration vector g as given in Eq. 329. 111 VVVVVVV gVRijk rstursrtru s tu rrr (329) The acceleration vector can be calculated by s ubstituting the expressions for the quantites ,,andVVVV rstu into Eq. 329. Before computing and substituting the quantities ,,andVVVV rstu into Eq. 329, the definitions in Eqs. 33035 are used to express the gravitational potential V in a compact form. a r (330) 0GM r (331) 10 (332) ,,,mmm nnmnmDstCRstSIst (333) 11,,,mmm nnn mmEstCRstSIst (334) 11,,,mmm nnn mmFstSRstCIst (335) Using Eqs. 33032, the recursion equations in Eqs. 33637 are formed. Using Eqs. 331, 333 and 336, the spherical harmonic re presentation of the gravitational potential V in Eq. 323 can be compactly represented as given in Eq. 338. 1for1nnn (336) 1nnra 11n nn ra (337) PAGE 46 46 max020,mm nnn n n nmVAuDst (338) The spherical harmonic representation of the gr avitational potential V in Eq. 338 is used with Eq. 329 to comput e the acceleration vector g The recursion relationships (Eqs. 339 and 340) and partials for ,m R st and ,m I st with respect to and st (Eqs. 341 and 342) are useful in computing the acceleration vector g 11,,,mmm R stsRsttIst (339) 11,,,mmm I stsIsttRst (340) 1,, ,mm mRstIst mRst st (341) 1,, ,mm mIstRst mIst st (342) Define 1234,,and aaaa as given in Eqs. 34346. Fr om Eqs. 329 and 34346, the acceleration vector g can be expressed in terms of 1234,,and aaaa as given in Eq. 347. Using Eqs. 33042, the quantities ,,and VVVV rstu are computed and substituted into Eqs. 34347 to obtain the expressions for 1234,,and aaaa as given in Eqs. 34851. 11 V a rs (343) 21 V a rt (344) 31 V a ru (345) 4VVVV a rstu s tu rrr (346) 4123 gaRaiajak (347) max11 20,mm nnn n n nmaAumEst a (348) max21 20,mm nnn n n nmaAumFst a (349) PAGE 47 47 max1 31 20,mm nnn n n nmaAuDst a (350) max0 4 1231 201 ,mm nnn n n nmn aAuDstsaaa ratu (351) The expressions for 1234,,and aaaa from Eqs. 34851 are substituted into Eq. 347 to obtain the acceleration vector g The acceleration vector obtained by the above procedure is coordinatized in ECEF frame. The above pro cedure for calculating th e acceleration vector g Eqs. 32051, however, has a numerical comput ation problem in the form of the polynomial m n A u which is addressed below. From Eq. 320, the polynomial m n A u is defined as 21 221 1 1 2! !!nm n m n nnmuunmnk d A nmndu (352) From Eq. 352, it can be seen that the computation of m nAu for high values of degree nand order m becomes cumbersome specifically due to the presence of terms nm and 21nm n nmud du This problem can be circumvented by using the recursion relationships for the polynomial m n A u given in Eqs. 35356 [16], [6]. 1 1 2 000 2 12121 21211 23nnnn uunnunu AAA nn (353) 1 11 221 2nn nnn uu AA n (354) 1 21 2 1 2 21 21 1 21 1 1mm nn m numuu AA nmnm nmnm uu A nmnm (355) 0,m nAumn (356) PAGE 48 48 It should be noted that the recursion relationships Eqs. 35355 cannot be used if the values of the polynomial foratleast 23,0m nAunmn are not provided. Therefore the values of polynomials ,23,0m n A unmn need to be calculated manually. Once the values of ,23,0m nAunmn are calculated, they can be used along with the recursion relationships Eqs. 35356 to generate ,3,0m nAunmn by following the procedure below [6]. 1. Equation 353 is used to generate 0,3nun A 2. Equation 354 is used to generate the diagonal elements ,3n nun A 3. Equations 355 and 356 are used to ge nerate all the remaining elements. The recursion relationships Eqs. 35356 have been proved to be the most numerically stable of all the recursion relati onships [16], [6]. The polynomials ,2360,m nAun 0 mn generated using the recurs ion relationships in Eqs. 35356 have a maximum percentage error of about 0.007983. This ma ximum percentage error was determined by computing the polynomial ,2360,0m nAunmn (i) using the recursion relationships Eqs. 35356 and (ii) using a MATLAB inbuilt function m nNu related to m n A u by the expression 2 2,where1 2 1m n mm nuukN uA 1for0 2for0 km km and then taking the difference between the co rresponding values of m nAu obtained by the above two methods and dividing it by the corresponding value of m n A u obtained by using th e second method (i.e., using m nNu). It should be noted that the maxi mum percentage error was computed for PAGE 49 49 0.990.99u because m nAu cannot be calculated at 1u using the second method (i.e., using m nNu). It can be concluded that the gravit y model described above i.e., Eqs. 32056 is singularityfree and numerically efficient. Inertia Model The inertia model presented here is used to calculate the mass pr operties of the launch vehicle (i.e., mass and mass mome nt of inertia of the launch vehicle). The mass moment of inertia of the launch vehicle is calculated about a Cartesian axes through the instantaneous center of mass of the launch vehicle. A simple but reasonable inertia model was developed and is being used in this research. The details of the developed inertia model are presented below. The launch vehicle is assumed to have (i) soli d strapon boosters, (ii) solid motors, (iii) liquid motors and (iv) a payload. To develop th e inertia characteristics, general geometries were considered. The simulator require s the inertia tensor of the launc h vehicle to be computed about the instantaneous centroidal axes of the launch ve hicle. Therefore the inertia tensors of the individual elements (solid strapon boosters, solid motors, liquid motors and payload) of the launch vehicle need to be computed about the instantaneous centroidal axes of the launch vehicle. This is calculated by first computing th e inertia tensor of the individual elements about their own centroidal axes and then using the pa rallel axis theorem, the inertia tensors of the individual elements about the in stantaneous centroidal axes of th e launch vehicle are computed. The sum of the instantaneous iner tia tensors of the individual el ements about the instantaneous centroidal axes of the launch vehicle is the inst antaneous inertia tensor of the entire launch vehicle about its instantaneous centroidal axes. These details are presented in the following subsections. PAGE 50 50 Strapon booster The solid strapon booster is divided into followi ng segments (i) a cylindrical segment, (ii) a parabolic nose cone and (iii) f our fins. To calculate the iner tia tensor of a strapon booster about the instantaneous centroidal axes of the laun ch vehicle, the inertia te nsors of the individual segments (of the strapon booster) about their ow n centroidal axes are calculated and then using the parallel axis theorem, the in ertia tensors of the individual segments about the instantaneous centroidal axes of the strapon boos ter are calculated. The sum of these inertia tensors of the individual segments about the instantaneous centroidal axes of the strapon booster is the instantaneous inertia tens or of the entire strapon booster abou t its instantaneous centroidal axes. This instantaneous inertia tensor of the strapon booster is used with the parallel axis theorem to calculate the instantaneous inertia tensor of the strapon booster about the instantaneous centroidal axes of the launch vehicle. The iner tia properties of the se gments of the strapon booster are presented below. Cylindrical segment The cylindrical segment of the strapon booster de picted in Fig. 32 is assumed to be made of a cylindrical shell structur e and solid propellant (the solid propellant is assumed to be distributed as a cylindrical shell). The cy lindrical shell structure of outer radius, os R inner radius, is R and height, c L is modeled as a solid cylinder of radius, os R and height, c L from which another cylinder of radius, is R and height, c L is removed (the bases and the centroidal axes of both cylinders coincide). The propellant shell is modeled similar to the cylindrical shell structure with a time varying radius of the inne r cylinder as the propellan t burns. Therefore the propellant shell of outer radius, op R inner radius, ip R and height, c L is modeled as a solid cylinder of radius, op R and height, c L from which another cylinder of time varying radius, ip R PAGE 51 51 and height, c L is removed (the bases and the centroidal ax es of both cylinders coincide). It is assumed that as the propellant burns the radius of the inner cylinder, ip R increases. When all the propellant is burnt, op ip R R. The densities of the materials of the cylindrical shell structure and the solid propellant are assumed to be and s p respectively. The mass of the cylindrical segment, c M is the sum of masses of the cylindrical shell structure and cylindrical propellant shell as show n in Eq. 357. The mass moment of inertias, and,CCCXXYYZZ I IIof the cylindrical segment about its centroidal axis is the sum of mass moment of inertias of the cylindr ical shell structure and the prope llant shell about the centroidal axis of the cylindrical segment as shown in Eqs. 358, 359 and 360, respectively. Figure 32. Cylindrical se gment of the strapon booster 2222 cscospcop isip M LRRLRR (357) PAGE 52 52 444411 22Cscospcop isip XX I LRRLRR (358) 222222 22222211 33 1212 11 33 1212Cscososcscc isis YY p copopcpcc ipipILRRLLRRL L RRLLRRL (359) 222222 22222211 33 1212 11 33 1212Cscososcscc isis ZZ p copopcpcc ipipILRRLLRRL L RRLLRRL (360) Therefore the inertia tensor of the assumed homogeneous cylindrical segment about its centroidal axis is 00 00 00C C CXX YY ZZI II I The instantaneous inertia tensor of the cylindrical segment about the instantaneous centro idal axes of the solid booster, C I is computed using the parallel axis theorem as shown in Eq.361. CM1TT c CCCCIIrrrr (361) The vector, Cr, in Eq. 361 is the position vector of the instantaneous center of mass of the solid strapon booster with respect to the ce nter of mass of the cy lindrical segment and 1 is the 3X3 identity matrix. Parabolic nose cone The nose cone of the solid strapon booster de picted in Fig. 33 is modeled as a solid parabolic nose cone of radius, N C R and length, N C L The density of material of the nose cone is assumed to be N C .The details of the mode l are presented below. The mass of the parabolic nose cone, N C M is given by the expression as shown in Eq. 362. The mass moment of inertias, and,NCNCNCXXYYZZ I II, of the parabolic nose cone about PAGE 53 53 its centroidal axis are given by the expressi ons as shown in Eqs. 363, 364 and 365 respectively. Figure 33. Parabolic nose cone 21 2 N CNCNCNC M RL (362) 21 3NCXX N CNC I MR (363) 2211 63NCYY N CNCNCIMRL (364) 2211 63NCZZ N CNCNCIMRL (365) Therefore the inertia te nsor of the parabolic nose cone about its centroidal axis is 00 00 00NC NC NCXX YY ZZI II I The instantaneous inertia tensor of the parabolic nose cone about the instantaneous centroidal axes of the solid booster, N C I is computed using the parallel axis theorem as shown in Eq. 366. NCM1TT NCNCNCNCNCIIrrrr (366) PAGE 54 54 The vector, N Cr, in Eq. 366 is the position vector of the instantaneous center of mass of the solid strapon booster with respect to the center of mass of th e parabolic nose cone. Fins The geometry of the fin is modeled by removing two triangular slabs 2 f and 3 f from a rectangular slab 1 f as shown in Fig. 34. Each strapon booster is assumed to have four fins of thickness f t towards the aft of the solid strapon booster. The masses of the three slabs 123,and f ff are given by the expressions as shown in Eqs. 367, 368 and 369. The mass moment of inertias of the slabs 123,and f ff about their own centroidal axes are given in Eqs. 37078. The expressions for 12and L L in Eqs. 37080 are given by 1TR L XC and 2TTR L XCC. The Xand Ycomponents of the location of center of mass of the fin, andCMCMXY, are given by the expressions as shown in Eqs. 379 and 380 respectively. The mass of the fin, f M ,is obtained by subtracting the mass of the slabs 2 f and 3 f from the mass of the slab 1 f as shown in Eq. 381. The moment of iner tia of the fin about its centroidal axes is obtained by subtracting the mome nt of inertias of slabs 2 f and 3 f from the corresponding moment of inertias of the slab 1 f as shown in Eqs 382, 383 and 384 respectively. 11 f ff M LtS (367) 21 2T f ff M XtS (368) 321 2 f ff M LtS (369) 11 12 21 122cm XX fff fS IMStMY (370) PAGE 55 55 Figure 34. Fin 111 12 2 11 122cm YY fff fL IMLtMX (371) 111 12 2 2 11 1222cmcm ZZ ff fLS IMLSMXY (372) 22 22 22112 18123cm XX fff fIMStMSY (373) 22 22 22111 18123cm YYTT fff fIMXtMXX (374) 22 222 22112 1833cmcm ZZTT ff fIMSXMXXSY (375) 33 32 22111 18123cm XX fff fIMStMSY (376) 233 32 22 21111 18123cm YY fff fIMLtMLLX (377) 233 322 22 21111 1833cmcm ZZ ff fIMSLMSYLLX (378) 22 1212 1211 33 2T CM T L XLLL X LXL (379) PAGE 56 56 12 1221 33 2T CM T L SXSLS Y LXL (380) 123 f fff M MMM (381) 123fXXXXXXXX f ffIIII (382) 123fYYYYYYYY f ffIIII (383) 123 fZZZZZZZZ f ffIIII (384) Therefore the inertia te nsor of the fin about its centroidal axis is 00 00 00f f f XX YY ZZI II I The instantaneous inertia tensors of the four fins about the inst antaneous centroidal axes of the solid booster, 1234and,, f fff I III, are computed using the paralle l axis theorem as shown in Eqs. 38588. 11111fM1TT fffffIIrrrr (385) 22222fM1TT fffffIIrrrr (386) 33333fM1TT fffffIIrrrr (387) 44444fM1TT fffffIIrrrr (388) The vectors, 1234,,and f fffrrrr, in Eqs. 38588 are th e position vectors of the instantaneous center of ma ss of the solid strapon booster with respect to the center of mass of the four fins respectively. The instantaneous inertia tensor of the entire solid strapon booster about its instantaneous centroidal axes, S I ; is obtained by summing the inertia tensors of the individual segments of the strapon booster (i .e., cylindrical segment, parabo lic nose cone and fins) about the instantaneous centroidal axis of the strapon booster as shown in Eq. 389. The instantaneous PAGE 57 57 mass of the solid booster is obtained by summing the masses of the individual segments of the strapon boosters as shown in Eq. 390. 1234SCNC f fff I IIIIII (389) 4SCNC f M M MM (390) As mentioned at the beginning of this mass moment of inertia section, the inertia tensor of the solid strapon booster need s to be computed about the inst antaneous centroidal axes of the launch vehicle. Therefore the instantaneous iner tia tensor of the entire solid strapon booster about the instantaneous centroidal axes of the launch vehicle C S I is computed using the parallel axis theorem as shown in Eq. 391. SM1TT CSSSSS SIIrrrr (391) The vector Sr, in Eq. 391 is the position vector of the instantaneous center of mass of the launch vehicle with respect to the instantaneous center of mass of the solid strapon booster. Liquid Engine The liquid engine depicted in Fig. 35 is assume d to be made of a cyli ndrical shell structure and liquid propellant. The cylindrical shell structure of outer radius, os R inner radius, is R and height, L is modeled as a solid cylinder of radius, os R and height, L from which another cylinder of radius, is R and height, L is removed (the bases and the centroidal axes of both cylinders coincide). The liquid propellant is modeled as a solid cylinder of radius, p R and height, L whose density (and mass) decreases with tim e (as the propellant burns). The densities of the materials of the cylindrical shell struct ure and the liquid propellant are assumed to be and s p respectively. PAGE 58 58 The mass of the liquid engine, L M is the sum of masses of the cylindrical shell structure and propellant cylinder as shown in Eq. 392. The mass moment of inertias, and,LLLXXYYZZ I II, of the liquid engine about its centroida l axis is the sum of mass moment of inertias of the cylindrical shell structure and the propellant cy linder about the centroidal axis of the liquid engine as shown in Eqs. 393, 394 and 395 respectively. Figure 35. Liquid engine 222 s ospp is L M LRRLR (392) 44411 22L s ospp is XX I LRRLR (393) 222222 22211 33 1212 1 3 12Csososs isis YY pppILRRLLRRL LRRL (394) PAGE 59 59 222222 22211 33 1212 1 3 12Csososs isis ZZ pppILRRLLRRL LRRL (395) Therefore the inertia tensor of the liquid engine about its centroidal axis is 00 00 00L L LXX YY ZZI II I As mentioned at the beginning of this inertia model section, the inertia tensor of the liquid engine needs to be computed about the instantaneous centroidal axes of the launch vehicle. Therefor e the instantaneou s inertia tensor of the liquid engine about the instantaneous centroidal axes of the launch vehicle, L I is computed using the parallel axis theorem as shown in Eq.396. LM1TT LLLLLIIrrrr (396) The vector, L r, in Eq. 396 is the position vector of the instantaneous center of mass of the launch vehicle with respect to the center of mass of the liquid engine. Solid Motor The solid motor model is similar to the cylindrical segment model of the solid strapon booster model discussed above in the strapon boos ter section. However, the details of model are presented again for the sake of completeness The solid motor depicted in Fig. 36 is assumed to be made of a cylindrical shell struct ure and solid propellant (t he solid propellant is assumed to be distributed as a cy lindrical shell). The cylindrical shell structure of outer radius, os R inner radius, is R and height, L is modeled as a solid cylinder of radius, os R and height, L from which another cylinder of radius, is R and height, L is removed (the bases and the centroidal axes of both cylinders coincide). The propellant shell is modeled similar to the PAGE 60 60 cylindrical shell structure with a time varying radi us of the inner cylinder as the propellant burns. Therefore the propellant shell of outer radius, op R inner radius, ip R and height, L is modeled as a solid cylinder of radius, op R and height, c L from which another cylinder of time varying radius, ip R and height, c L is removed (the bases and the centroidal axes of both cylinders coincide). It is assumed that as the propell ant burns, the radius of the inner cylinder, ip R increases. When all the propellant is burnt, op ip R R The densities of the materials of the cylindrical shell structure and the solid propellant are assumed to be and s p respectively. 2222 sospop isip SM M LRRLRR (397) 444411 22SMsospop isip XX I LRRLRR (398) Figure 36. Solid motor PAGE 61 61 222222 22222211 33 1212 11 33 1212SMsososs isis YY popopp ipipILRRLLRRL L RRLLRRL (399) 222222 22222211 33 1212 11 33 1212SMsososs isis ZZ popopp ipipILRRLLRRL L RRLLRRL (3100) The mass of the solid motor, SM M is the sum of masses of th e cylindrical shell structure and cylindrical propellant sh ell as shown in Eq. 397. The mass moment of inertias, and,SMSMSMXXYYZZ I II, of the solid motor about its centroidal axis is the sum of mass moment of inertias of the cylindrical sh ell structure and the pr opellant shell about the centroidal axis of the solid motor as shown in Eqs. 398, 399 and 3100 respectively. Therefore the inertia tensor of the so lid motor about its centroidal axis is 00 00 00SM SM SMXX YY ZZI II I The instantaneous inertia tensor of the solid motor about the instantaneous centroidal axes of the launch vehicle, SM I is computed using the parallel axis theorem as shown in Eq.3101. CM1TT SMSMSMSMSMIIrrrr (3101) The vector, SMr, in Eq. 3101 is the position vector of the instantaneous center of mass of the launch vehicle with respect to the inst antaneous center of mass of the solid motor. Payload The payload model presented below is simila r to the parabolic nos e cone model of the strapon booster model discussed ab ove in the strapon booster secti on. However, the details of model are presented again for the sake of completeness. The payl oad depicted in Fig. 37 is PAGE 62 62 modeled as a solid parabolic nose cone of radius, P R and length, P L The density of material of the nose cone is assumed to be P .The details of the mode l are presented below. Figure 37. Payload The mass of the parabolic nose cone, P M is given by the expression as shown in Eq. 3102. The mass moment of inertias, and, P PPXXYYZZ I II, of the parabolic nose cone about its centroidal axis are given by the expressions as shown in Eqs. 3103, 3104 and 3105 respectively. 21 2PPPP M RL (3102) 21 3PXXPP I MR (3103) 2211 63PYYPPP I MRL (3104) 2211 63P Z ZPPP I MRL (3105) PAGE 63 63 Therefore the inertia te nsor of the parabolic nose cone about its centroidal axis is 00 00 00P P P XX YY ZZI II I The instantaneous inertia te nsor of the payload about the instantaneous centroidal axes of the launch vehicle, P I is computed using the parallel axis theorem as shown in Eq. 3106. PM1TT PPPPPIIrrrr (3106) The vector, Pr, in Eq. 3106 is the position vector of the instantaneous center of mass of the solid strapon booster with respect to the center of mass of th e parabolic nose cone. Drag Coefficient Model The drag force acting on the launch vehicle depe nds on the drag coefficient, crosssectional area perpendicular to the flow, de nsity of the medium and the ve locity of the launch vehicle relative to the medium. The drag force act ing on a launch vehicle can be expressed as 1 2DragD F CAvv where DC is the drag coefficient; A is the crosssectional area perpendicular to the flow; is the density of the medium; andvv are the speed and velocity of the launch vehicle relative to the medium. The simulator requires the drag force to be coordinatized in the inertial frame. The expression for the drag force coordinatized in inertial frame is given as 1 2ii DragD F CAvv. In this research, the drag coefficient is calcula ted using a drag coeffici ent model obtained from the Missile Datcom database [13]. The descripti on of the drag coefficient model is given below. The total drag coefficient acting on the launc h vehicle can be expressed as the sum of body zerolift wave drag coeffi cient, body base drag coeffici ent and body skin friction drag PAGE 64 64 coefficient i.e., ,,ooooDDD B odyBodyWaveBaseBodyFrictionCCDCC where oD B odyCis the body zerolift drag coefficient or simply the drag coefficient; ,oD B odyWaveC is the body zerolift wave drag coefficient; oD B aseC is the body base drag coefficient and ,oD B odyFrictionC is the body skin friction drag coefficient. The expressions for the body zer olift wave drag coefficient, body base drag coefficient and body skin fric tion drag coefficient are given below. 1.69 1 2 ,1.830.5 1.59tanoD BodyWave NC L M d for M>1 (Based on Bonney reference, 1tan in radians) ,0.25oD BaseCoastC M if M>1 and 2,0.120.13oD BaseCoast M C, if M<1 Ref0.25 1oe D BasePoweredA C SM if M>1 2 Ref10.120.13oe D BasePoweredA CM S if M<1 0.2 ,0.053oD BodyFrictionLM C dqL (Based on Jerger reference, turbulent boundary layer, q in psf, L in feet) where NL is the nose length; d is the rocket diameter; L is the rocket body length; e A is the nozzle exit area; RefSis the reference area; q is the dynamic pressure and M is the mach number. Center of Pressure Model The aerodynamic forces acting on the launch vehicl e are assumed to be concentrated at the center of pressure of th e launch vehicle. In order to co mpute the moments of the aerodynamic forces (i.e., drag and lift) about the instantaneous center of mass, the center of pressure needs to be calculated. The center of pressure for a launch vehicle can be computed using Barrowman PAGE 65 65 equations [2]. The assumptions considered in the calculation of the center of pressure using Barrowman equations are given below [3]. The angle of attack is very near zero. The flow is steady state and subsonic. Flow over launch vehicl e is potential flow. Point of the nose is sharp. Fins are thin plates with no cant. The launch vehicle is a rigid body. The launch vehicle is axially symmetric. The first two assumptions severely restrict the applicability of the Barrowman equations for calculation of the center of pressure for a real time launch vehicle. De spite these restrictions, the Barrowman equations are used for the calcul ation of the center of pressure of the launch vehicle as it gives a rough estimate of the location of the center of pressure of the launch vehicle. However, to be more conservative, the location of the center of pressure can be chosen so that the launch vehicle is statically stable throughout its flight. Fo r static stability the center of pressure must be located aft of the center of ma ss and the static margin must be at least one caliber (i.e., distance between the center of pressu re and the center of mass of the launch vehicle is equal to the largest diameter of the launch vehi cle) [1]. It should be noted that a modern full scale launch vehicles do not rely on aerodynamics for stability. Active control systems like thrust vector control provides st ability and control to the modern full scale launch vehicles [5]. Specifically, the above assumption of one caliber st atic margin was made in order to compute the aerodynamic moments acting on the launch vehicle. However, the procedure to compute the center of pressure of a launch vehicle using Ba rrowman equations is presented below for the sake of completeness. The procedure for calculating the center of pr essure of a launch ve hicle using Barrowman equations involves dividing the vehicle into separa te sections, analyzing e ach section separately, PAGE 66 66 analyzing the interference effects between sect ions and then recombining the results of the separate analyses to obtain the final answer i.e., the normal force coefficients and center of pressure locations of the separate sections w ith respect to a common poi nt are calculated and then the center of pressure of the launch vehicl e is calculated by using the formula in Eq. 3107 [2], [4]. The term cp x in Eq. 3107 is the location of the cen ter of pressure of the launch vehicle with respect to an arbitrary point O; N iC is the normal force coefficient of the ith section of the launch vehicle; i x is the location of the center of pressure of the ith section of the launch vehicle with respect to the same arbitrary point O. i N i i cp N i iCx x C (3107) A launch vehicle is assumed to be composed of the following sections (i) nose, (ii) cylindrical body, (iii) conical s houlder, (iv) conical boattail and (v) fins (tail section). The normal force coefficients and locations of the cente r of pressure for the above sections (i)(v) are given below. Nose The normal force coefficient for all shapes of the nose is always equal, nose NC = 2. However, the location of the center of pressure is different for different shapes of the nose. The locations of the center of pressure measured from fore of the nose for various shapes are given below. (i) conical nose: 2 3cpLx; (ii) ogive nose: 0.466cpLx and (iii) parabolic nose: 1 2cpLx PAGE 67 67 Cylindrical Body The normal force coefficient for a cylindrical bod y at low angles of attack is zero, i.e., cylindrical_bodyNC = 0. The location of center of pressure for a cylindrical body at low angle of attack is undefined or not necessary. Conical Shoulder The normal force coefficient for a conical s houlder is given in Eq. 3108. The term d jn Eq. 3108 is the diameter of the nose of the la unch vehicle. The loca tion of the center of pressure of the conical s houlder as shown in Fig. 38 is given in Eq. 3109. 22 21 conical_shoulder2Ndd C dd (3108) 1 2 2 1 21 1 3 1cpd d L x d d (3109) Figure 38. Conical shoulder PAGE 68 68 Conical Boattail The normal force coefficient for a conical bo attail is given in Eq. 3110. The term d in Eq. 3110 is the diameter of the nose of the la unch vehicle. The loca tion of the center of pressure of the conical s houlder as shown in Fig. 39 is given in Eq. 3111. 22 21 conical_boattail2Ndd C dd (3110) 1 2 2 1 21 1 3 1cpd d L x d d (3111) Figure 39. Conical Boattail Fins (Tail Section) The normal force coefficient for a fin geometry as shown in Fig. 310 is given in Eq. 3112. The term R in Eq. 3112 is the radius of the tail section of the launch vehicle; d is the diameter of the nose of the launch vehicle; n is the number of fins andinterK is the fin interference factor to be consider ed in order to account for the in terference between the fins and the tail section (i.e., laun ch vehicle body section to which the fi ns are attached) and is calculated PAGE 69 69 using Eq. 3113. The location of the center of pre ssure of the tail secti on of the launch vehicle (i.e., the section of the launch vehicle body to whic h the fins are attached) as shown in Fig. 310 is independent of the number of fins and is give n by Eq. 3114 (where all the fins attached to the tail section are of same size and shape). 2 inter 2fin4 2 11Ns n d K L abC (3112) interK = 1 R sR for 3 or 4 fins per stage (3113) = 0.5 1 R sR for 6 fins per stage 2 1 36RT RRT cpRT RTRTCC XCC xCC CCCC (3114) Figure 310. Fin and Tail section The WGS84 Ellipsoid Model WGS (World Geodetic System) defines a reference model for the earth, for use in geodesy and navigation. WGS84 is used as the reference ellipsoid in this research. WGS84 was PAGE 70 70 last revised in 2004 and will be valid up to about 2010 [11]. The model presented below computes the geodetic coordinates (i.e., geodetic la titude, longitude and altitude) from the Earth centered Earth fixed (ECEF) coordinates. The ge odetic coordinates and th e ECEF coordinates of an arbitrary point P are depicted in Fig. 311. The computed al titude is used to calculate the density of air at that altitude. The calculated de nsity of air in turn is necessary to compute the aerodynamic forces acting on the launch vehicle. Th e details of the model are given below [12]. The geodetic coordinates ,, h are defined as follows. Geodetic latitude, is the angle between the ellipsoidal normal N and the equatorial plane. Longitude, is the angle measured in the equato rial plane between the prime meridian (i.e., the Xaxis of the ECEF frame) and the pr ojection of the point P onto the equatorial plane. Altitude, h, is the distance between the surface of the ellipsoid and the point P measured along the ellipsoidal normal. Figure 311. Geodetic Ellipsoid and Geodetic coordinates of an arbitrary point P The WGS84 ellipsoid is define d by the following parameters. Semimajor axis length of the ellipsoid, a = 6378137.0 meters Semiminor axis length of the ellipsoid, b = 6356752.3142 meters PAGE 71 71 Flatness of the ellipsoid, ab f a = 0.0034 Eccentricity of the ellipsoid, 2 eff = 0.0818 The relationship between ECEF coordinates ,, x yz and geodetic coordinates ,, h of an arbitrary point of interest P can be expressed by Eqs. 3115117. coscos xNh (3115) cossin yNh (3116) 21sinzNeh (3117) Given the ECEF coordinates of an arbitrary point P, geodetic coordinates of that point can be computed using the following algorithm. 1. From Eqs. 3115 and 3116, the longitude of the point P can be calculated using Eq.3118. 1tan y x (3118) 2. Initialize h = 0 and Naand define 22 p xy. 3. Iterate the following expressions until converges as given below. 2sin 1z Neh (3119) 2 1sin tan zeN p (3120) 221sina N e (3121) cosp hN (3122) Using the above algorithm (i.e., Eqs. 311912 2), geodetic coordinate s of an arbitrary point P can be computed from the ECEF coordinates of P. This chapter presented the models develope d for calculating quanti ties such as gravity, inertia, mass, center of pressure and drag coeffici ent. These models along with the equations of motion of the launch vehicle developed in Chapter 2 were implemented in MATLAB to form a PAGE 72 72 sixdegreeoffreedom launch vehi cle simulator. The results a nd discussions of a simulation performed by this simulator are presented in the next chapter. PAGE 73 73 CHAPTER 4 SIMULATION RESULTS AND DISCUSSION The dynamic equations (translational and ro tational), kinematic equations, gravity, inertia, drag coefficient, center of pressure and WGS84 ellipsoid models developed in Chapters 2 and 3 were implemented in MATLAB, these (alo ng with the other MATL AB scripts/functions) collectively form the sixdegreeoffreedom la unch vehicle simulator. Using the developed simulator, a DELTA II launch vehicle model was si mulated for a fictitious constant axial thrust profile. This chapter presents the results and discussions of that simulation. Following the results and discussions, the details of the attempt to reconstruct the thrust profile of an ATLAS IIA launch vehicle from the flight data provide d by NASA Kennedy Space Ce nter are presented. This reconstructed thrust profile is required for the valid ation of the simulator. It was shown that the simulator cannot be validated due to the lack of availability of sufficient realtime data necessary to reconstruct the thrust profile as required by the simulator. Simulation In this section, the results and discussions of a DELTA II la unch vehicle model simulation are presented. The various parameters used for this simulation are listed below. 1. Launch vehicle DELTA II [7] o Four solid strapon boosters o First stage solid motor o Second stage liquid engine o Third stage liquid engine o Payload enclosed in a parabolic nose cone 2. Thrust: Constant magnitude and acting along the longitudinal axis of the launch vehicle. The thrust profile for this simulati on can be seen in Fig. 41 (Tx vs time). 3. Place of launch: Kennedy Sp ace Center, Merrit Island, Florida. 4. Date of launch: 15 June 2004. PAGE 74 74 5. Time of launch: 1640:30 UTC. The developed simulator requires the following data to be placed in individual text files (a) thrust profile of the main section of the launc h vehicle (the section consisting of solid motors, liquid engines and payloads), (b) thrust profil es of the individual strapon boosters, (c) the instantaneous mass of the entire launch vehicle, (d) the position of the instantaneous center of mass of the launch vehicle with respect to its in itial center of mass and its time derivatives, (e) instantaneous inertia tens or of the launch vehicle about its instantaneous center of mass and (f) the positions of the nozzles of all the thrusting elements (solid motors, liquid engines and strapon boosters) with respect to the instantaneous center of mass of the launch vehicle. The procedure used to populate the text files stated in (a)(f) for this simulation is explained in brief in the following (1)(4) steps. (1) The assumed fictitious thrust profile was used to populate the text files stated in (a) and (b). (2) The text files populated in the step (1) i.e., the thrust profiles text files were used with Eqs. 41 and 42 to compute the instantaneous mass of the launch vehicle and populate the text file stated in (c). i ioi s p ttdm TgI dt (41) i it ttittdm mmt dt (42) The term itT in Eq. 41 represents the magnitude of the thrust provided by the ith thrusting element at time t, itT is calculated from the text files stated in (a) and (b); itdm dt in Eqs. 41 and 42 is the rate of change of mass of the ith thrusting element at time t; the term i s pI in Eq. 41 is the specific impulse of the ith thrusting element; tim in Eq. 42 is the mass of the ith thrusting PAGE 75 75 element at time t; ittm in Eq. 42 is the mass of the ith thrusting element at time tt ; og in Eq. 42 is the acceleration due to gravity at the surface of the Earth. (3) The launch vehicle instantaneous mass text file popula ted in step (2) is used along with the inertia model discussed in Chapter 3 and DELTA II launch vehicle geomet ry (obtained from DELTA II payload planners guide [7]) to compute the instantaneous mass, th e initial and the instantaneous center of mass and the instantaneous inertia tensor of the launch vehi cle and populate the text files stated in (d) and (e). (4) The text file stated in (d) i.e., the instantaneous center of mass position text file populated in step (3) is used with the DELTA II launch vehi cle geometry (obtained from DELTA II payload planners guide [7]) to compute positions of the nozzles of all the thrusting elements with respect to the instantaneous cente r of mass of the launch vehicle and populate the text file stated in (f). (For detailed DELTA II launch vehicle parameters and calculations, please refer to APPENDIX B.) The simulation resu lts and discussions are presented below. Figure 41 shows the following parameters of the launch vehicle as a function of time. Ms is the instantaneous mass of the booster. M1 is the instantaneous mass of the first stage. M2 is the instantaneous mass of the second stage. M3 is the instantaneous mass of the third stage. rc is the location of the instantaneous cen ter of mass of the launch vehicle with respect to its initial center of mass. Tx is the magnitude of the thrust acting on the launch vehicle as a function of time and is user defined. The strapon booster provides thrust to th e launch vehicle by burning and expelling its fuel/propellant. Therefor e the instantaneous mass of the strapon booster Ms decreases as a function of time which can be seen in Fig 41. The mass of the strapon booster after 64th second PAGE 76 76 is zero because the strapon boos ter is jettisoned at the 64th second, and therefore its mass contribution to the launch vehicle is zero. Similar kind of argume nts as the above explains the instantaneous mass plots of the firs t, second and third stages of the launch vehicle. The plot of rc in Fig. 41 is the Xcomponent of the position ve ctor of the instantaneous center of mass with respect to the initial center of mass in the vehi cle frame. The vehicle frame was described in Chapter 2 and its origin coincides with the initial center of mass of the la unch vehicle. The Xaxis of the vehicle frame coincides with the long itudinal axis of the launc h vehicle and is positive forwards (i.e., towards the nose of the launch vehicle) As the fuel is burnt and expelled from the bottom stages/strapon boosters, the instantaneous center of mass of the launch vehicle shifts up (i.e., towards the nose of the launch vehicle) comp ared to the instantaneous center of mass of the launch vehicle at the previous instant because re latively more mass is concentrated towards the top than towards the bottom when compared to the previous instant in time. Therefore rc increases with time as fuel/propellant is burnt and expelled. The Y and Z components of the position vector of the instantaneous center of mass with respect to the initial center of mass are zero always because the center of mass of the launch vehicle lies always on its longitudinal axis. The Tx in Fig. 41 shows the X component of the th rust vector acting on the launch vehicle as a function of time. The time history of the thrust vector is called the thrust profile. Generally, thrust profiles are (i) coordina tized in the vehicle frame and ( ii) designed based on the trajectory of mission (i.e., user defined). In this present case, constant axial thrust is considered for all the stages and strapon boosters of th e launch vehicle. Therefore th e X component is equal to the magnitude of the thrust vector (s ince thrust is axial, Y and Z components of th e thrust vector in the vehicle frame are zero), and the magnitude of the thrust remains constant between two successive jettisons as seen in Fig. 41. Specifically, between the liftoff and the 64th second, PAGE 77 77 Figure 41. Various parameters of th e launch vehicle as a function of time 2673736 N of thrust is provided by the four strapon boosters and the first stage liquid engine; between the 64th second and 269th second, 889644 N of thrust is pr ovided by the first stage liquid engine; between the 269th second and the 700th second, 43657 N of thrust is provided by the second stage liquid engi ne and between the 700th second and the 787th second, 66367 N of thrust is provided by the third stage solid motor. After the 787th second, no thrust is provided to the launch vehicle i.e., Tx = 0. Figure 42 shows the x, y, z components and magnitude of the velo city vector of the launch vehicle in the J2000 inertial frame as a function of time. From the velocity magnitude plot, it can be seen that the magnitude of veloc ity of the launch vehicle in creases slowly and after the 64th second from the time of launch, the magnitude decreases by a small amount and then PAGE 78 78 increases until the 289th second from the time of launch and then suddenly decreases, a similar decrease and increase in the ve locity magnitude can be observe d again at the 880th second from the time of launch. This response can be analyzed as follows. The strapon boosters and the first stage of the launch vehicle provide 2673736 N of th rust to the launch vehicle from the beginning of the launch until the 64th second at which the strapon boosters are jettisoned. During this time, the fuel/propellant of the strapon boosters and the first stag e of the launch vehicle is burnt and expelled which decreases the mass of the launch vehicle. Since, the thrusts of the strapon boosters and the first stage of the launch vehicl e are constant, it can be seen that the same magnitude of the thrust (2673736 N) is propell ing a launch vehicle of decreasing mass and therefore it results in accelerating the launch vehi cle and hence the increase in the magnitude of velocity of the launch vehicle is seen. At the 64th second, the strapon boosters are jettisoned while the first stage of the launch vehicle c ontinues to provide a thrust of 889644 N until the 289th second at which the first stage is jettisoned. Just before the 64th second, both the strapon booster and the first stage of the launch vehicle provide a total thrust of 2673736 N to the launch vehicle and just after the 64th second (and until the 289th second), only the first stage provides a thrust of 889644 N. Therefore it can be seen th at after the 64th sec ond, lesser thrust (when compared to the thrust just before the 64th second ) is propelling the launch vehicle of almost the same mass (since, the mass of th e jettisoned boosters frames are very small compared to the mass of the entire launch vehicle) and therefore it results in decelerating the launch vehicle and hence the slight decrease in the magnitude of the velocity can be seen after the 64th second. The magnitude of velocity of the launch vehicle k eeps on decreasing until a point where the constant thrust of 889644 N from the first stage of the la unch vehicle is sufficient to accelerate the launch vehicle of decreasing mass and therefore the ma gnitude of velocity increases until the 289th PAGE 79 79 second. At 289th second, the first stage is jet tisoned and the magnitude of velocity decreases again and then the cycle repeats. This respons e can be explained by means of similar arguments as stated above. Figure 42. Velocity of the launc h vehicle in the inertial frame Figure 43 shows the x, y, z components of th e position vector of the launch vehicle and the position of the launch vehicle in 3D in the J2000 in ertial frame (discussed in Chapter 2) as a function of time. At the time of launc h, the launch vehicle points in the +XI, YI, +ZI direction as seen from the J2000 inertial frame which is depicted in Fig. 44. Since the thrust is assumed to be constant and axial, the thrust force always acts along the longitudi nal axis of the launch vehicle, and since the longitudinal axis of the launch vehicle is al most parallel to a vector from the center of Earth to the center of launch vehicle, the thrust force continuously accelerates the PAGE 80 80 Figure 43. Position of the launc h vehicle in the inertial frame Figure 44. Launch vehicle during the time of launch as seen from the J2000 inertial frame PAGE 81 81 launch vehicle approximately in the +XI, YI, +ZI direction and therefor e the launch vehicle position changes accordingly (i.e., position becomes more positive in X direction, more negative in Y direction and more positive in the Z direction) which can be seen in Fig. 43. The 3D plot in Fig. 43 shows the launch vehicle starti ng at point A and reaching point B after 1000 seconds from the time of launch.. The rates at which the X, Yand Zcomponents of the position vector change in Fig. 43 are dictated by the corresponding X, Yand Zcomponents of the velocity vector in Fig 42. Figure 45 shows the moments of inertia of the launch vehicle about its instantaneous center of mass as a function of time. xxyyzz,andIII are the (instantaneous) moments of inertia of the launch vehicle about the X, Y, and Zaxis respectively of a coordinate frame whose orientation is the same as that of the vehicle frame and passes through the instantaneous center of mass of the launch vehicle. Since the launc h vehicle is symmetric about the Xaxis, Figure 45. Moments of inertia of the launch vehicle about its instan taneous center of mass PAGE 82 82 yyzzandII are equal. The change in th e inertia of the launc h vehicle is due to the change in (i) mass of the launch vehicle and (ii) location of the instantaneous center of mass of the launch vehicle. These changes in turn are due to the consumption of propellant and jettisoning of the consumed stages of the launch vehicle. Figure 46 shows the moments of inertia of a strapon booster of th e launch vehicle about (i) a centroidal axis passing through the instanta neous center of mass of the launch vehicle and (ii) a centroidal axis passing th rough its instantaneous center of mass as a function of time. Figure 46. Moments of inertia of the strapon booster about the instanta neous center of mass of the launch vehicle and about its instantaneous center of mass xxyyzz ss s,andIII are the (instantaneous) moments of inertia of the strapon booster about the X, Y, and Zaxis respectively of a coordinate frame whose orientation is the same as that of PAGE 83 83 the vehicle frame and passes through its instantaneous center of mass. ccccsxxyy s,II cczz sandI are the (instantaneous) moments of inertia of the strapon booster about the X, Y, and Zaxis respectively of a coordinate frame whose orientation is the same as that of the vehicle frame and passes through th e instantaneous cente r of mass of the launch vehicle. The change in the inertia of the st rapon booster is due to the chan ges in (i) mass of the strapon booster and (ii) locations of the instantaneous cen ter of mass of both the launch vehicle and the strapon booster. These changes in turn are due to the consumption of propellant and jettisoning of the consumed stages of the launch vehicle. Figure 47 shows the moment of inertia of the first stage of the launc h vehicle about (i) a centroidal axis passing through th e instantaneous center of mass of the launch vehicle and (ii) a centroidal axis passing through its instantaneous center of ma ss as a function of time. xxyy 1 1,II zz 1andI are the (instantaneous) moments of inertia of th e first stage about the X, Y, and Zaxis respectively of a coordinate fr ame whose orientation is the same as that of the vehicle frame and passes through its instantaneous center of mass. cccc1xxyy 1,II cczz 1andI are the (instantaneous) moments of inertia of the first stage about the X, Y, and Zaxis respectively of a coordinate frame whose orientation is the same as that of the vehicle frame and passes through the instantaneous cent er of mass of the launch vehicle. From Fig. 47, it can be seen that there is a slight incr ease (as opposed to decrease) in cc1yyI cczz 1andI around 64th second. This is due to jettisoning of the strapon boosters. When the strapon boosters are jettisoned, there is a sudden increase in the di stance between the instantaneous center of mass of the first stage and the instantaneous center of ma ss of the launch vehicle. Therefore the increase PAGE 84 84 in the inertia of the first stage about the cente r of mass of the launch vehicle (computed using parallelaxis theorem)is greater than the decreas e in the inertia of the first stage due to the consumption of its propellant which results in a slight increase in cc1yyI cczz 1andI when the strapon boosters are jettisoned. After this slight increase in inertia of the first stage, its Figure 47. Moment of inertia of the first stage about the inst antaneous center of mass of the launch vehicle and about its instantaneous center of mass inertia decreases again until it is jettisoned. The ch ange in the inertia of the first stage is due to the changes in (i) mass of the firs t stage and (ii) locations of th e instantaneous center of mass of both the launch vehicle and the first stage. Thes e changes in turn are due to the consumption of propellant and jettisoning of the consumed stages of the launch vehicle. PAGE 85 85 Figure 48 shows the moments of inertia of th e second stage of the la unch vehicle about (i) a centroidal axis passing through th e instantaneous center of mass of the launch vehicle and (ii) a centroidal axis passing through its instantaneous center of ma ss as a function of time. Figure 48. Moment of inertia of the second st age about the instantaneous center of mass of the launch vehicle and about its instantaneous center of mass xx 2I yy 2,Izz 2andI are the (instantaneous) moments of inertia of the second stage about the X, Y, and Zaxis respectively of a coordinate frame whose orientation is the same as that of the vehicle frame and passes through its instantaneous center of mass. cccc2xxyy 2,II cczz 2andI are the (instantaneous) mome nts of inertia of the second stage about the X, Y, and Zaxis respectively of a coordinate frame whose orientation is the same as that of the vehicle frame and passes through the instantaneous center of mass of the launch vehicle. The change in PAGE 86 86 the inertia of the second stage is due to the changes in (i) mass of both the second stage and (ii) locations of the instantaneous center of mass of bo th the launch vehicle and the second stage. These changes in turn are due to the consumpti on of propellant and jett isoning of the consumed stages of the launch vehicle. Figure 49 shows the moments of inertia of the third stage of the launc h vehicle about (i) a centroidal axis passing through th e instantaneous center of mass of the launch vehicle and (ii) a centroidal axis passing through its instantaneous center of ma ss as a function of time. Figure 49. Moment of inertia of the third stage about the inst antaneous center of mass of the launch vehicle and about its instantaneous center of mass xx 3I yy 3,Izz 3andI are the (instantaneous) moments of inertia of the third stage about the X, Y, and Zaxis respectively of a coordinate frame whose orientation is the same as that of PAGE 87 87 the vehicle frame and passes through its instantaneous center of mass. cccc3xxyy 3,II cczz 3andI are the (instantaneous) moment s of inertia of the third stage about the X, Y, and Zaxis respectively of a co ordinate frame whose orient ation is the same as that of the vehicle frame and passes through the instantaneou s center of mass of the launch ve hicle. The change in the inertia of the third stage is due to the changes in (i) mass of both the third stage and (ii) locations of the instantaneous center of mass of both the launch vehicle and the third stage. These changes in turn are due to the consumption of propellant and jettisoning of the consumed stages of the launch vehicle. Validation Validation of the simulator requires comparing the simulator output to the realtime flight data. This can be done by provi ding realtime input of a launch vehicle for a particular mission (date, place and time of launch and thrust profile data coordinatized in vehicle frame) to the simulator and comparing the output of the simula tor (attitude and position vectors) to the realtime flight data (attitude and position vectors) of the same launch vehicle for the same mission. As the validation requires all re altime input and flight data and since the realtime thrust profile of a launch vehicle is not available to publi c (an ITAR issue), one of the alternatives is to extract the approximate thrust profile from (i) place, date and time of launch, (ii) position, velocity, and acceleration vectors of the launch vehi cle as a function of time, (iii) attitude of the launch vehicle as a function of time and (iv) launch vehicle configura tion and ascent profile Assuming that all the above data is available, approximate thrust profile of the launch vehicle can be calculated by following the procedure presented below. The translational equation of motion of a launch vehicle disc ussed in Chapter 2 is given in Eq. 43. In Eq. 43, I a, is the acceleration vector of the launch vehicle coordinatized in the PAGE 88 88 inertial frame. The forces and,,IIII g Drag ThrustLift F FFF in Eq. 43 can be computed using Eqs. 23739 developed in Chapter 2 I extIIIII g Drag ThrustLiftmaFFFFF (43) For simplicity, the frame of re ference notation is dropped i.e., a instead of I a. Also, the notation t x is used to represent a vector quantity x at a time t. Using Eq. 43 and the above notation, the expression for tThrustF is given in Eq. 44. The relation between thrust and rate of change of mass is given in Eq. 45. ttt t ttg Drag ThrustLiftmaFFFF (44) ot tt Thrust s pdm dt gIF (45) In Eq. 43, og is the acceleration due to gravit y at the surface of the Earth and t s p I is the specific impulse. The algorithm (steps 1) presented be low is used to compute the thrust profile. 1. Calculate the mass at time t using 1 1 tt tdm mmdt dt 2. Calculate the thrust vector at time t using ttt t ttg Drag ThrustLiftmaFFFF 3. Calculate the rate of change of mass at time t using ot tThrust t s pdm gI dtF 4. Repeat steps 1,2 and 3 to calculate the thru st vectors at time intervals t+1, t+2, ,until the final time interval t+n (i.e., step 1 to calculate the mass at time t+1 and step 2 to calculate the thrust vector at time t +1 and step 3 to calculate the rate of change of mass at time t+1). 5. The thrust vectors at time inte rvals t = 0, 1, 2, t+n represents the approximate thrust profile of the launch vehicle (the la unch vehicles configuration and ascent profile play a critical role in the computation of thrust profile). The thrust profile obtained by the above procedure is coordinatized in the inertial frame, but, the simulator requires the input thrust profile to be c oordinatized in the vehicle frame. The attitude PAGE 89 89 of the launch vehicle can be used to transform the thrust profile from the inertial frame to the vehicle frame using the transfor mation matrix from the inertial frame to the vehicle frame R IC (Eqs. 28 and 29) discussed in Chapter 2. Validation Attempt Partial data of an ATLASIIA launch vehi cle for a realtime mission is available. An attempt was made to validate the simulator as described at the beginning of this (validation) section using the available data. The attempt proved to be futile; the details of the attempt are given below. The following realtime flight data of an ATLASIIA launch vehicle available for the validation of the simulator (i) pla ce, date and time of launch, (ii) position and velocity vectors of the launch vehicle as a function of time and (iii) launch vehicle configuration and ascent profile. The realtime attitude of the ATLASIIA launch vehicle is not available. As mentioned previously, the simulator requires the thrust profil e to be coordinatized in the vehicle frame. Since the realtime attitude of th e launch vehicle is unavailable, it is not possible to coordinatize the thrust profile in the vehicle frame as requi red by the simulator. Therefore the simulator cannot be validated with the available data. Therefore in addition to the available data, inst rumented flight data (i.e., accelerations and gyro on board the vehicle) or realtime thrust pr ofile of the launch vehicle is necessary to validate the simulator. Figure 410 shows the graphical depiction of the need for instrumental data or thrust vector in the vehicle frame. PAGE 90 90 Figure 410. The need for instrumental data or thrust vector in the vehicle frame PAGE 91 91 CHAPTER 5 CONCLUSION AND FUTURE WORK Conclusions This research presents the development of a sixdegreeoffreedom launch vehicle simulator which along with the communications link to ol being developed by UCF forms the trajectory and link margin analysis tool. This tr ajectory and link margin an alysis tool is one of the crucial support tools for STARS which provi de the ability to quickly analyze new (or changes in) concepts and ideas, an option not easily accomplished w ith hardware only. Chapter 1 provided an overview of the United States range safety, STARS concept and the trajectory and link margin analysis which motivat ed the development of a sixdegreeoffreedom launch vehicle simulator. Chapter 2 modeled th e full dynamics of an expendable launch vehicle by developing its kinematic and the dynamic equa tions. It was shown that the kinematic equation expressed in terms of Euler angles had problems in the form of a singularity and numerical computation. Thes e problems were circumvented by expressing the kinematic equation in terms of quaternions. Following th e development of the kinematic equation, the dynamic equations of an expendable launch ve hicle were developed by accounting for the variability in its mass and geometry. Chapter 3 developed the models for the calculation of quantities such as gravity, inertia, center of pressure and drag coefficient required for solving the equations of motion presented in Chapter 2. The gravity model u tilized the spherical harmonic re presentation of the gravitational potential to account for the variability in Eart hs mass distribution and used EGM96 (360 X 360) spherical harmonic coefficients and WGS84 Earth ellipsoid m odel. The gravity model was shown to be singularityfree and num erically efficient. A novel way of calculating the variable mass/inertial properties of a launch vehicle was de veloped. This inertia model was a simple and PAGE 92 92 approximate model and considered general geometri es to develop the inertia characteristics of a launch vehicle. The drag coefficient model used in this research was obtained from the Missile Datcom database. The models developed in Chapters 2 and 3 were implemented in MATLAB to form a sixdegreeoffreedom launch vehicle simulator. The results and discussions of a simulation performed using this simulator were presented in Chapter 4. These resu lts with the support of the discussions were proven to be logically tr ue. However, it was shown that the simulator cannot be validated as the critic al data needed for the validati on could not be provided due to ITAR restrictions. This inability to validate the developed simulator prevented the completion of the development of the trajectory and link margin analysis tool required by the STARS as a support tool. Future work The future work would first include validating the simulator. Once validated, the simulator developed as a result of this research can be us ed as a basis/framework for the development of a future highfidelity sixdegreeoffreedom simu lator. The future work/modifications to the current simulator would include the following. The aerodynamic model used to compute quant ities such as the coefficient of drag, coefficient of lift, center of pr essure and other aerodynamic coefficients will be modified. The inertia model will be modified to handle nonsymme tric configurations of the launch vehicles. The developed sixdegreeoffreedom simulato r will be interfaced with the MATLAB based telemetry link margin analysis tool which is being developed by UCF to form the trajectory and link margin analysis tool. A GUI will be developed for the trajectory and link margin analysis tool and will be finally interfaced with Satellite Tool Kit (STK). PAGE 93 93 APPENDIX A MATLAB FUNCTIONS AND SCRIPT % % A SixDegreeofFreedom Launch Vehicle % % % Simulator for Range Safety Analysis % % % Author : Sharath Chandra Prodduturi, Graduate Student, SSG % University of Florida. % Date : 06232007 close all clear all clc % Constants % Angular velocity of Earth omega_earth = (2*pi)/(24*60*60); % Standard Gravitational Parameter for Earth mu_earth = 398600.4418*10^9; input( 'enter the year of launch (20XX) = ); input ( 'enter the month of launch (1 12) = ); input( 'enter the day of launch (1 31) = ); input( 'enter the hour of launch (1 24) = ); input( 'enter the minute of launch (1 60) = ); input( 'enter the second of launch (1 60) = ); % Function for checking the input dates and times % If the input data isnt valid, the program terminates. [year,month,day,hour,min,sec,check] = timecheck(year1,month1,day1,hour1, ... min1,sec1); if isempty(check)  (check == 0) error( 'Error in the Input. Execute the file again...!!' ) else % Function for Calculating the time elapsed (in sec) from JAN 1 2000 to the %input date [total_days,t0] = time(year,month,day,hour,min,sec); end disp( 'Enter the Launch site' ) disp( '1 for Kennedy Space Center, Merrit Island, Florida' ) disp([ '2 for Vandenberg Air Force Base, Space Command 30th Space Wing' ... ', California' ]) disp( '3 for Virginia Space Flight Center, Wallops Island, Virginia' ) disp( '4 for Edwards Air Force Base, California' ) disp( '5 for Poker Flat Research Range, Alaska' ) disp( '6 for Alaska Spaceport, Kodiak Island, Alaska' ) disp( '7 for Mojave Civilian Aerospace Test Center, California' ) disp( '8 for White Sands Space Space Harbor, Las Cruces, New Mexico' ) input( 'enter the launch site = ); % Initial Euler Angles input( 'Enter the initial Euler angles [th_1; th_2; th_3] = ); eu_in = [0; 90; 0]; PAGE 94 94 eu_in = pi/180 eu_in; % Function for extracting the Latitude and Longitude of the selected Launch Site [long,lat] = place(place_1); load vehicle.m vehicle_parameters = vehicle; %ODE45, Outputs of the function are x,y,z components of position,velocity of the initial Center of mass of the %rocket in the Inertial frame and Euler angles % % Initial conditions for the ODE solver initial = initial_conditions(lat, long, eu_in, omega_earth, t0); options = []; % odeset('RelTol',3e2,'AbsTol',3e5); [t,y] = ode113( 'position' [0:1:1000], initial, options, t0, omega_earth, ... lat, long, mu_earth, vehicle_parameters); format long % plots of position of the rocket in the Inertial frame (Nonrotating % geocentric equatorial reference frame) figure(1) subplot(2,2,1); plot(t,y(:,1)) xlabel( 't' ); ylabel( 'x' ); subplot(2,2,2); plot(t,y(:,2)) xlabel( 't' ); ylabel( 'y' ); subplot(2,2,3); plot(t,y(:,3)) xlabel( 't' ); ylabel( 'z' ); subplot(2,2,4); plot3(y(:,1),y(:,2),y(:,3)) xlabel( 'x' ); ylabel( 'y' ); zlabel( 'z' ); % plots of velocity of the rocket in the Inertial frame (Nonrotating % geocentric equatorial reference frame) figure(2) subplot(2,2,1); plot(t,y(:,4)) xlabel( 't' ); ylabel( 'velocity_x' ); subplot(2,2,2); plot(t,y(:,5)) xlabel( 't' ); ylabel( 'velocity_y' ); subplot(2,2,3); PAGE 95 95 plot(t,y(:,6)) xlabel( 't' ); ylabel( 'velocity_z' ); function [year,month,day,hour,min,sec,check] = timecheck(year1,month1,day1,hour1,min1,sec1); % Function for checking the input dates and times % If the input data isnt valid then the program terminates. if ((year1 < 2000)  (year1 > 2099)  (month1 > 12)  (day1 > 31)  (day1 < 0)  (hour1 > 24)  (hour1 < 0) (min1 > 60)  (min1 < 0) (sec1 > 60)  (sec1 < 0)) year = 0; month = 0; day = 0; hour = 0; min = 0; sec =0; elseif ((rem((year12000),4) == 0) & (month1 == 2) & (day1 > 29)) year = 0; month = 0; day = 0; hour = 0; min = 0; sec =0; elseif ((rem((year12000),4) ~= 0) & (month1 == 2) & (day1 > 28)) year = 0; month = 0; day = 0; hour = 0; min = 0; sec =0; elseif ((month1 == 4  6  9  11) & (day1 > 30)) year = 0; month = 0; day = 0; hour = 0; min = 0; sec =0; else year = year1; month = month1; day = day1; hour = hour1; min = min1; sec = sec1; end if ((year == 0)  (month == 0)  (day == 0)) check = 0; else check = 1; end function [total_days,t0] = time(year,month,day,hour,min,sec) %Function for Calculating the time elapsed (in sec) from JAN 1 2000 to the %input date days = 0; y1 = (year2000)/4; y = (year2000)+(ceil(y1))/365; y2 = rem((year2000),4); for k=1:month1 switch k case {1,3,5,7,8,10,12} days = days + 31; case {4,6,9,11} days = days + 30; case 2 switch y2 case 0 days = days + 29; otherwise days = days + 28; PAGE 96 96 end end end total_days = y*365 + days + (day1); % Time in seconds from JAN 1 2000 00:00:00.00 to the input time given. t0 = total_days*86400 + hour*3600 + min*60 + sec; % Time in seconds from JAN 1 2000 12:00:00.00 to the input time given. t0 = t0 43200; function [long,lat] = place(place_1) % % Function for extracting the Latitude and Longitude of the selected Launch Site % If the input data isnt valid, the program terminates. % % calling sequence: % [long,lat] = place(place_1) % Define Variables: % place_1 Launch Site % long Longitude of the Launch Site % lat Latitude of the Launch Site switch place_1 case 1 long = 81.0; lat = 28.5; case 2 long = 120.35; lat = 34.4; case 3 long = 75.5; lat = 37.8; case 4; long = 118; lat = 35; case 5 long = 147.8; lat = 64.9; case 6 long = 146; lat = 67.5; case 7 long = 118.2; lat = 35; case 8 long = 106.8; lat = 32.3; case 9 disp( 'You have not entered a proper choice' ) disp( 'You can select the Launch site of your choice' ) long = input( 'Enter the geographic longitude of the Launch site = ); lat = input( 'Enter the geocentric latitude of the Launch site = ); otherwise error( 'Error in the Input. Execute the file again...!!' ) PAGE 97 97 end long = long*pi/180; % in rad lat = lat*pi/180; % in rad % [long,lat] = place(place); % Ref : http://www.spacetoday.org/Rockets/Spaceports/LaunchSites.html function initial = initial_conditions(lat, long, eu_in, omega_earth, t0) % % This function calculates the initial conditions of the launch vehicle. % These initial conditions are required by the ODE solver. The initial % conditions produced by this function are: % (i) Initial position of the launch vehicle in the Inertial frame % (ii) Initial velocity of the launch vehicle in the Inertial frame % (iii) Initial quaternion parameters (i.e., q0, q1, q2 and q3) % Initial position of the rocket in the inertial frame e = 0.0818; a = 6378137.0; N = a/sqrt(1(e*sin(long))^2); % x_g, y_g, z_g are the components of initial postion of the launch vehicle % in the rotating geocentric frame x_g = N*cos(lat)*cos(long); y_g = N*cos(lat)*sin(long); z_g = N*(1e^2)*sin(lat); % Initial angle between X_I and X_G H_G = omega_earth*t0; % Initial transformation matrix from E_I to E_G C_GI = [cos(H_G) sin(H_G) 0; sin(H_G) cos(H_G) 0; 0 0 1]; % P_I is the initial position of the launch vehicle in the inertial frame P_I = C_GI'*[x_g; y_g; z_g]; % V_I is the initial velocity of the launch vehicle in the Inertial frame w_e = [0; 0 ; omega_earth]; V_I = cross(w_e, P_I); % Initial Quaternion Parameters q0 = cos(eu_in(3)/2)*cos(eu_in(2)/2)*cos(eu_in(1)/2) + ... sin(eu_in(3)/2)*sin(eu_in(2)/2)*sin(eu_in(1)/2); q1 = cos(eu_in(3)/2)*cos(eu_in(2)/2)*sin(eu_in(1)/2) ... sin(eu_in(3)/2)*sin(eu_in(2)/2)*cos(eu_in(1)/2); q2 = cos(eu_in(3)/2)*sin(eu_in(2)/2)*cos(eu_in(1)/2) + ... sin(eu_in(3)/2)*cos(eu_in(2)/2)*sin(eu_in(1)/2); q3 = sin(eu_in(3)/2)*cos(eu_in(2)/2)*cos(eu_in(1)/2) ... cos(eu_in(3)/2)*sin(eu_in(2)/2)*sin(eu_in(1)/2); %Initial conditions for the ODE solver PAGE 98 98 initial = [P_I(1), P_I(2), P_I(3), V_I(1), V_I(2), V_I(3), q0, q1, q2, ... q3,0,0,0]; function ydot = position(t, y, initial, t0, omega_earth, lat, long, ... mu_earth, vehicle_parameters) % % ODE45, Outputs of the function are x,y,z components of position, velocity % vectors of the initial Center of mass of the rocket in the Inertial % frame and Euler angles ydot = zeros(size(y)); % % Coordinate frame Transformations % Initial angle between X_I(Y_I) and X_G(Y_G) H_G = omega_earth*t0; % Instantaneous angle between X_I(Y_I) and X_G(Y_G) H_G_instant = omega_earth*(t + t0); % Instantaneous Transformation matrix from E_I to E_G C_GI_instant = [cos(H_G_instant) sin(H_G_instant) 0; sin(H_G_instant) cos(H_G_instant) 0; 0 0 1]; % Initial Transformation matrix from E_I to E_G C_GI = [cos(H_G) sin(H_G) 0; sin(H_G) cos(H_G) 0; 0 0 1]; % Transformation matrix from E_G to E_V C_VG = [sin(lat)*cos(long) sin(lat)*sin(long) cos(lat); sin(long) cos(long) 0; cos(lat)*cos(long) cos(lat)*sin(long) sin(lat)]; % Transformation matrix from E_V to E_R C_RV = [2*y(7)^2+2*y(8)^21, 2*y(8)*y(9)+2*y(7)*y(10), ... 2*y(8)*y(10)2*y(7)*y(9); 2*y(8)*y(9)2*y(7)*y(10), 2*y(7)^2+2*y(9)^21, ... 2*y(9)*y(10)+2*y(7)*y(8); 2*y(8)*y(10)+2*y(7)*y(9), 2*y(9)*y(10)2*y(7)*y(8), ... 2*y(7)^2+2*y(10)^21]; % Transformation matrix from E_I to E_R C_RI = C_RV*C_VG*C_GI; % % This function generates a flag which is used in calculation of the mass % properties, forces and moments. flag = flag_generator(t); %  PAGE 99 99 % Mass properties % Instantaneous Center of mass and Instantaneous mass calculations [r_c, v_c, a_c, m_r] = cm_function(t, flag); % Instantaneous Mass Moment of Inertia Tensor/Matrix Calculations [I_xx I_xy I_xz I_yy I_yz I_zz] = Inertia_function(t, flag); % % Forces and Moments % Thrust [Thrust_Force Thrust_Moment] = Thrust_function(t, flag, C_RI, ... vehicle_parameters); % Aerodynamic forces [Drag_Force Drag_Moment Lift_Force Lift_Moment] = ... AeroForces_function(t, y, C_GI, C_RI, flag, vehicle_parameters); % Gravity force [g_Force, g_Moment] = Gravity_function(y(1:3), C_GI_instant, m_r); % Moment_External = Thrust_Moment + Drag_Moment + Lift_Moment + g_Moment; % Rotational Equation of motion ydot(11) = (Moment_External(1) (I_zz I_yy)*y(12)*y(13))/I_xx; ydot(12) = (Moment_External(2) (I_xx I_zz)*y(11)*y(13))/I_yy; ydot(13) = (Moment_External(3) (I_yy I_xx)*y(11)*y(12))/I_zz; omega_r = [y(11); y(12); y(13)]; omega_dot_r = [ydot(11); ydot(12); ydot(13)]; % % Relation between the quaternion rates and the angular velocity M = 0.5*[0 y(11) y(12) y(13); y(11) 0 y(13) y(12); y(12) y(13) 0 y(11); y(13) y(12) y(11) 0] [y(7); y(8); y(9); y(10)]; ydot(7) = M(1); ydot(8) = M(2); ydot(9) = M(3); ydot(10) = M(4); % % Instantaneous Center of mass of the Rocket in the Vehicle frame R_c = r_c; % Instantaneous Velocity of the Center of mass of the Rocket in the Vehicle % frame dR_c = v_c; PAGE 100 100 % Instantaneous Acceleration of the Center of mass of the Rocket in the % Vehicle frame d2R_c = a_c; % % y(1),y(2),y(3) are the x,y,z components of the position of the Initial % center of mass of the rocket in the Inertial frame % y(4),y(5),y(6) are the x,y,z components of the velocity of the Initial % center of mass of the rocket in the Inertial frame % ydot(4),ydot(5),ydot(6) are the x,y,z components of the acceleration of % the Initial center of mass of the rocket in the Inertial frame ydot(1) = y(4); % y(1) = R_c1(1) ydot(2) = y(5); % y(2) = R_c1(2) ydot(3) = y(6); % y(3) = R_c1(3) % Translational Equation of motion Force_External = Thrust_Force + Drag_Force + Lift_Force + g_Force; YDOT = (1/m_r)*(Force_External C_RI'*(d2R_c + 2*cross(omega_r,dR_c) ... + cross(omega_dot_r,R_c) + cross(omega_r,(cross(omega_r,R_c))))); ydot(4) = YDOT(1); ydot(5) = YDOT(2); ydot(6) = YDOT(3); % function flag = flag_generator(t) % % This function generates a flag based on the thrust profile on the Launch % Vehicle This function checks if the thrust is zero during the immediate % past instant or immediate future instant in order to avoid interpolation % load Thrust_data/Thrust.txt load Thrust_data/Nozzle.txt % Checking if the thrust is zero during the immediate past instant or % immediate future instant in order to avoid interpolation [tmi, te, tpi] = numerical_rounding(t, 1); t_x_check1 = interp1(Thrust(:,1), Thrust(:,2), tpi); t_y_check1 = interp1(Thrust(:,1), Thrust(:,3), tpi); t_z_check1 = interp1(Thrust(:,1), Thrust(:,4), tpi); t_x_check2 = interp1(Thrust(:,1), Thrust(:,2), te); t_y_check2 = interp1(Thrust(:,1), Thrust(:,3), te); t_z_check2 = interp1(Thrust(:,1), Thrust(:,4), te); if (norm([t_x_check1; t_y_check1; t_z_check1]) == 0) flag = 1; elseif (norm([t_x_check2; t_y_check2; t_z_check2]) == 0) flag = 2; PAGE 101 101 else flag = 0; end function [xmi, xe, xpi] = numerical_rounding(x, n) % % This function rounds the input x to the specified n decimal points. n >=1 % The outputs of this function are the rounded x and two numerical values % which are greater and lesser than the rounded x by (10)^(n) % if (n >= 1) x1 = x fix(x); x1 = x1*(10)^n; x1 = fix(x1); x1 = x1*(10)^(n); xe = fix(x) + x1; xmi = xe (10)^(n); xpi = xe + (10)^(n); else error( 'n should be atleast 1' ) end function [r_c, v_c, a_c, m_r] = cm_function(t, flag) % % This function calculates/outputs the following: % (a) Position Vector % (b) Velocity Vector % (c) Acceleration Vector % of the Instantaneous Center of Mass of the Launch Vehicle % w.r.t the Initial Center of Mass of the the Launch Vehicle and % (d) Instantaneous Mass of the Launch Vehicle % % The input to this function is the time 't' of interest and the flag % generated by the flag_generator. % % This function requires the following data in text files in a folder % named 'Mass_properties' (without quotes) : % % (1) Time and the components of the Position, Velocity and Acceleration % vectors of the Instantaneous Center of Mass of the Launch Vehicle % w.r.t the Initial Center of Mass of the Launch Vehicle in a text % file named 'CM_data' (without quotes) in the following format: % % t1 Px(t1) Py(t1) Pz(t1) Vx(t1) Vy(t1) Vz(t1) Ax(t1) Ay(t1) Az(t1) % t2 Px(t2) Py(t2) Pz(t2) Vx(t2) Vy(t2) Vz(t2) Ax(t2) Ay(t2) Az(t2) % . . % . . % % where: % Px, Py and Pz are the X, Y, Zcomponents of the Position Vector % of the Instantaneous Center of Mass of the Launch Vehicle w.r.t. the % Initial Center of Mass respectively % % vx, Vy and Vz are the X, Y, Zcomponents of the Velocity Vector % of the Instantaneous Center of Mass of the Launch Vehicle w.r.t. the PAGE 102 102 % Initial Center of Mass respectively % Ax, Ay and Az are the X, Y, Zcomponents of the Acceleration Vector % of the Instantaneous Center of Mass of the Launch Vehicle w.r.t. the % Initial Center of Mass respectively % % (2) Time and Instantaneous Center of Mass of the Launch Vehicle in a % text file named 'Mass_data' (without quotes) in the following format: % % t1 M(t1) % t2 M(t2) % % % load Mass_properties/CM_data.txt load Mass_properties/Mass_data.txt [tmi, te, tpi] = numerical_rounding(t, 1); switch flag case 0 time = t; case 1 time = tpi; case 2 time = te; otherwise error([ 'Invalid value of "flag" in cm_function, Please check' ... the value of "flag" generated by flag_generator' ]) end Px = interp1(CM_data(:,1), CM_data(:,2), time); Py = interp1(CM_data(:,1), CM_data(:,3), time); Pz = interp1(CM_data(:,1), CM_data(:,4), time); Vx = interp1(CM_data(:,1), CM_data(:,5), time); Vy = interp1(CM_data(:,1), CM_data(:,6), time); Vz = interp1(CM_data(:,1), CM_data(:,7), time); Ax = interp1(CM_data(:,1), CM_data(:,8), time); Ay = interp1(CM_data(:,1), CM_data(:,9), time); Az = interp1(CM_data(:,1), CM_data(:,10), time); m_r = interp1(Mass_data(:,1), Mass_data(:,2), time); r_c = [Px; Py; Pz]; v_c = [Vx; Vy; Vz]; a_c = [Ax; Ay; Az]; function [I_xx I_xy I_xz I_yy I_yz I_zz] = Inertia_function(t, flag) % % This function calculates/outputs the components of the Inertia tensor % of the launch vehicle. The input to this function is the time 't' of % interest and the flag generated by the flag_generator. % % This function requires the following data in a text file named % 'Inertia_data' (without quotes) placed in a folder named % 'Mass_properties' (without quotes) : % % (i) Time and the components of the Inertia tensor of the launch vehicle PAGE 103 103 % about the instantaneous center of the launch vehicle in the % following format: % % t1 I_xx(t1) I_xy(t1) I_xz(t1) I_yy(t1) I_yz(t1) I_zz(t1) % t2 I_xx(t2) I_xy(t2) I_xz(t2) I_yy(t2) I_yz(t2) I_zz(t2) % . % . % load Mass_properties/Inertia_data.txt [tmi, te, tpi] = numerical_rounding(t, 1); switch flag case 0 time = t; case 1 time = tpi; case 2 time = te; otherwise error([ 'Invalid value of "flag" in Inertia_function, Please' ... check the value of "flag" generated by flag_generator' ]) end I_xx = interp1(Inertia_data(:,1), Inertia_data(:,2), time); I_xy = interp1(Inertia_data(:,1), Inertia_data(:,3), time); I_xz = interp1(Inertia_data(:,1), Inertia_data(:,4), time); I_yy = interp1(Inertia_data(:,1), Inertia_data(:,5), time); I_yz = interp1(Inertia_data(:,1), Inertia_data(:,6), time); I_zz = interp1(Inertia_data(:,1), Inertia_data(:,7), time); function [Thrust_Force Thrust_Moment] = Thrust_function(t, flag, C_RI, ... vehicle_parameters) % % This function calculates the Total Thrust & Thrust Moment (about the % Instantaneous Center of Mass of the Launch Vehicle) of the Launch Vehicle % (including Boosters). The output is cordinatized in the vehicle frame % % The input parameters (to be passed) to this function are the time 't' % of interest, number of Boosters, the flag generated by the % flag_generator, the transformation matrix C_RI and the vehicle_parameters % array % % This function requires the following data in text files in a folder % named 'Thrust_data' (without quotes) : % % (1) Time and the Thrust vector (coordinatized in the Vehicle frame) % of the Launch Vehicle in a text file named 'Thrust.txt' (without % quotes) in the following format: % % t1 Thrust_x(t1) Thrust_y(t1) Thrust_z(t1) % t2 Thrust_x(t2) Thrust_y(t2) Thrust_z(t2) % . % . % % (2) Time and the (Position) vector from the instantaneous center of mass % of the launch vehicle (coordinatized in the Vehicle frame) to the PAGE 104 104 % nozzle of the launch vehicle through which the consumed fuel / gases % are expelled in a text file named 'Nozzle.txt' (without quotes) in % the following format: % % t1 Nozzle_x(t1) Nozzle_y(t1) Nozzle_z(t1) % t2 Nozzle_x(t2) Nozzle_y(t2) Nozzle_z(t2) % . % . % % (3) Time and Thrust vector (coordinatized in the Booster / % Vehicle frame) of the i'th booster in a text file named % 'Thrust_bi.txt' (without quotes) where 'i' in '_bi' represents the % i'th booster (i = 1,2,...,n) in the following format. (For example, % 1st booster data should be in a text file named 'Thrust_b1.txt' % (without quotes), 2nd booster data in 'Thrust_b2.txt' (without % quotes) and so on..) % % t1 Thrust_x_bi(t1) Thrust_y_bi(t1) Thrust_z_bi(t1) % t2 Thrust_x_bi(t2) Thrust_y_bi(t2) Thrust_z_bi(t2) % . % . % % (4) Time and the (Position) vector from the instantaneous center of mass % of the launch vehicle (coordinatized in the vehicle frame) nozzle of % the i'th booster through which the consumed fuel / gases are expelled % in a file named 'Nozzle_bi.txt' (without quotes) where 'i' in '_bi' % represents the i'th booster (i = 1,2,...,n) in the following format. % (For example, 1st booster data should be in a text file named % 'Nozzle_b1.txt' (without quotes), 2nd booster data in 'Nozzle_b2.txt' % (without quotes) and so on..) % % t1 Nozzle_x_bi(t1) Nozzle_y_bi(t1) Nozzle_z_bi(t1) % t2 Nozzle_x_bi(t2) Nozzle_y_bi(t2) Nozzle_z_bi(t2) % . % . % % Initialize to zero vector(s) Thrust_Force = [0;0;0]; Thrust_Moment = [0;0;0]; % Miscellaneous Parameters Calculations [L_vehicle L_strap D_vehicle D_strap A_e_vehicle A_e_strap Ln_vehicle ... Ln_strap Ref_A_vehicle Ref_A_strap boosters] = ... misc_calculations_function(t, flag, vehicle_parameters); % Boosters Thrust Data for i = 1:boosters dummy = int2str(i); load (strcat( 'Thrust_data/' 'Thrust_b' ,dummy, '.txt' )) load (strcat( 'Thrust_data/' 'Nozzle_b' ,dummy, '.txt' )) eval(sprintf( 't_x = interp1(Thrust_b%d(:,1), Thrust_b%d(:,2), t);' ... i, i)); eval(sprintf( 't_y = interp1(Thrust_b%d(:,1), Thrust_b%d(:,3), t);' ... i, i)); eval(sprintf( 't_z = interp1(Thrust_b%d(:,1), Thrust_b%d(:,4), t);' ... i, i)); PAGE 105 105 eval(sprintf( 'n_x = interp1(Nozzle_b%d(:,1), Nozzle_b%d(:,2), t);' ... i, i)); eval(sprintf( 'n_y = interp1(Nozzle_b%d(:,1), Nozzle_b%d(:,3), t);' ... i, i)); eval(sprintf( 'n_z = interp1(Nozzle_b%d(:,1), Nozzle_b%d(:,4), t);' ... i, i)); eval(sprintf( 'thrust_b%d_x = t_x;' i)); eval(sprintf( 'thrust_b%d_y = t_y;' i)); eval(sprintf( 'thrust_b%d_z = t_z;' i)); eval(sprintf( 'nozzle_b%d_x = n_x;' i)); eval(sprintf( 'nozzle_b%d_y = n_y;' i)); eval(sprintf( 'nozzle_b%d_z = n_z;' i)); end for i = 1:boosters eval(sprintf([ 'Thrust_Force = Thrust_Force + [thrust_b%d_x; ... 'thrust_b%d_y; thrust_b%d_z];' ],i,i,i)); eval(sprintf([ 'Thrust_Moment = Thrust_Moment + ... 'cross(([nozzle_b%d_x; nozzle_b%d_y; nozzle_b%d_z])' ... ,[thrust_b%d_x; thrust_b%d_y; thrust_b%d_z]);' ] ... ,i,i,i,i,i,i)); end % Launch Vehicle Data load Thrust_data/Thrust.txt load Thrust_data/Nozzle.txt switch flag case {1, 2} t_x = 0; t_y = 0; t_z = 0; n_x = 0; n_y = 0; n_z = 0; case 0 t_x = interp1(Thrust(:,1), Thrust(:,2), t); t_y = interp1(Thrust(:,1), Thrust(:,3), t); t_z = interp1(Thrust(:,1), Thrust(:,4), t); n_x = interp1(Nozzle(:,1), Nozzle(:,2), t); n_y = interp1(Nozzle(:,1), Nozzle(:,3), t); n_z = interp1(Nozzle(:,1), Nozzle(:,4), t); otherwise error([ 'Invalid value of "flag" in Thrust_function, Please' ... check the value of "flag" generated by flag_generator' ]) end Thrust_Force = Thrust_Force + [t_x;t_y;t_z]; Thrust_Moment = Thrust_Moment + cross(([n_x; n_y; n_z]), ... [t_x;t_y;t_z]); % Coordinatizing the thrust force in the inertial frame Thrust_Force = transpose(C_RI)*Thrust_Force; function [L_vehicle L_strap D_vehicle D_strap A_e_vehicle A_e_strap ... PAGE 106 106 Ln_vehicle Ln_strap Ref_A_vehicle Ref_A_strap N_s] ... = misc_calculations_function(t, flag, vehicle_parameters) % % This function calculates/outputs various miscellaneous parameters: % 1. Length of the Launch Vehicle (and StrapOn Boosters) % 2. Diameter of the Launch Vehicle (and StrapOn Boosters) % 3. Nozzle exit area of the Launch vehicle (and StrapOn Boosters) % 4. Nose length of the Launch Vehicle (and StrapOn Boosters) % 5. Reference area of the Launch Vehicle (and StrapOn Boosters) % 6. Position vector of the Initial Center of Mass of the Launch Vehicle % w.r.t the nose of the Launch Vehicle expressed in the Vehicle Frame % (Body Coordinate Frame) % % The input parameters (to be passed) to this function are the time 't' % of interest, the flag generated by the flag_generator and the % vehicle_parameters array % % This function requires the following data in a text file named % 'Length_data' (without quotes) placed in a folder named % 'Mass_properties' (without quotes) : % % (i) Time and the length of the Launch Vehicle (excluding the StrapOn % Boosters) in the following format: % % t1 L(t1) % t2 L(t2) % % % % Length of the Launch Vehicle load Mass_properties/Length_data.txt [tmi, te, tpi] = numerical_rounding(t, 1); switch flag case 0 L_vehicle = interp1(Length_data(:,1), Length_data(:,2), t); case 1 L_vehicle = interp1(Length_data(:,1), Length_data(:,2), tpi); case 2 L_vehicle = interp1(Length_data(:,1), Length_data(:,2), te); otherwise error([ 'Invalid value of "flag" in misc_caclculations_' ... 'function, Please check the value of "flag" generated' ... by flag_generator' ]) end % Length of the Strapon Booster(s) L_strap = vehicle_parameters(1); % Diameter of the Launch Vehicle D_vehicle = vehicle_parameters(2); % Diameter of the Strapon Booster(s) D_strap = vehicle_parameters(3); PAGE 107 107 % Nozzle exit area of the Launch Vehicle A_e_vehicle = vehicle_parameters(4); % Nozzle exit area of the Strapon Booster(s) A_e_strap = vehicle_parameters(5); % Nose length of the Launch vehicle Ln_vehicle = vehicle_parameters(6); % Nose length of the Strapon Booster(s) Ln_strap = vehicle_parameters(7); % Reference area of the Launch Vehicle Ref_A_vehicle = (pi/4)*D_vehicle^2; % Reference area of the StrapOn Booster(s) Ref_A_strap = (pi/4)*D_strap^2; % Number of boosters N_s = vehicle_parameters(8); function [Drag_Force Drag_Moment Lift_Force Lift_Moment] = ... AeroForces_function(t, y, C_GI, C_RI, flag, vehicle_parameters) % % This function calculates/outputs the Atmospheric Forces and Moments: % 1. Drag Force (coordinatized in the inertial frame) % 2. Drag Moment (coordinatized in the vehicle frame) % 3. Lift Force (coordinatized in the inertial frame) % 4. Lift Moment (coordinatized in the vehicle frame) % % The input parameters (to be passed) to this function are the time 't' % of interest, the flag generated by the flag_generator and the vehicle % parameter array % % This function requires the following data in text files in a folder % named 'Mass_properties' (without quotes) : % % (1) Position Vector of the Instantaneous Center of Mass of the ith % Strapon Booster w.r.t the Instantaneous Center of Mass of the % entire Launch Vehicle expressed in the Vehicle Frame in a text file % named 'cm_bi.txt' (without quotes) where 'i' in '_bi' represents % the i'th booster (i = 1,2,...,n) in the following format. % (For example, 1st booster data should be in a text file named % 'cm_b1.txt' (without quotes), 2nd booster data in 'cm_b2.txt' % (without quotes) and so on..) position_I = y(1:3); velocity_I = y(4:6); % Position vector of the Launch vehicle in the ECEF frame position_ECEF = C_GI*position_I; % Velocity vector of the Launch Vehicle in the Vehicle Reference frame velocity_R = C_RI*velocity_I; % Miscellaneous Parameters Calculations PAGE 108 108 [L_vehicle L_strap D_vehicle D_strap A_e_vehicle A_e_strap Ln_vehicle ... Ln_strap Ref_A_vehicle Ref_A_strap boosters] = ... misc_calculations_function(t, flag, vehicle_parameters); % Density of air and Coefficient of Drag Calculation if t == 0 if boosters ~= 0 C_D_strap = 0; end C_D_main = 0; rho_air = 1.225; else [geod_lat, geod_long, h] = wgs84(position_ECEF); if h < 600000 rho_air = AtmDens2(h/1000); if boosters ~= 0 C_D_strap = C_Drag(rho_air, h, velocity_I, Ln_strap, ... D_strap/2, A_e_strap, Ref_A_strap, L_strap); end C_D_main = C_Drag(rho_air, h, velocity_I, Ln_vehicle, ... D_vehicle/2, A_e_vehicle, Ref_A_vehicle, L_vehicle); else rho_air = 0; if boosters ~= 0 C_D_strap = 0; end C_D_main = 0; end end [tmi, te, tpi] = numerical_rounding(t, 1); switch flag case 0 time = t; case 1 time = tpi; case 2 time = te; otherwise error([ 'Invalid value of "flag" in AeroForces_function, ... 'Please check the value of "flag" generated by ... 'flag_generator' ]) end % % Drag % Forces % Drag force acting on the Launch Vehicle (excluding boosters) Drag_Force_main = 0.5*C_D_main*rho_air*Ref_A_vehicle* ... norm(velocity_R)*velocity_R; % Drag force acting on the boosters Drag_Force_strap = 0.5*C_D_strap*rho_air*Ref_A_strap* ... norm(velocity_R)*velocity_R; PAGE 109 109 % Total Drag force acting on the Launch Vehicle (i.e., including boosters) Drag_Force = Drag_Force_main + boosters*Drag_Force_strap; % Moments % Moment due to the Drag force acting on the Launch Vehicle about the % instantaneous center of mass of the launch vehicle including the % strapon boosters. (The drag force is assumed to be acting at the center % of pressure of the launch vehicle Drag_Moment = cross([D_vehicle;0;0], Drag_Force); % Coordinatize the Drag Force in the Inertial Frame Drag_Force = transpose(C_RI)*Drag_Force; % % Lift % Forces % Lift force acting on the Launch Vehicle (excluding boosters) Lift_Force_main = [0;0;0]; % Lift force acting on the boosters Lift_Force_strap = [0;0;0]; % Total Lift force acting on the Launch Vehicle (i.e., including boosters) Lift_Force = Lift_Force_main + boosters*Lift_Force_strap; % Moments % Moment due to the Lift force acting on the Launch Vehicle about the % instantaneous center of mass of the launch vehicle including the % strapon boosters. (The Lift force is assumed to be acting at the center % of pressure of the launch vehicle) Lift_Moment = cross([D_vehicle;0;0], Lift_Force); % Coordinatize the Lift Force in the Inertial Frame Lift_Force = transpose(C_RI)*Lift_Force; % function [geod_lat, geod_long, alt] = wgs84(R_ECEF) % This function calculates the Geodetic Latitude, Longitude and the % Altitude from the rectangular ECEF coordinates. % % Form: [geod_lat, geod_long, alt] = wgs84(R_ECEF) % % Input to this function is the Position Vector of the point of interest in % the ECEF frame. % % Outputs of this function are: % 1. geod_lat = Geodetic Latitude of the Point of Interest. % 2. geod_long = Longitude of the Point of Interest. % 3. alt = Altitude (is the distance along the ellipsoidal normal, away from % the interior of the ellipsoid, between the surface of the ellipsoid and % the point of interest) PAGE 110 110 % Reference: The Global Positioning System & Inertial Navigation by Jay % A.Farell & Matthew Barth. x = R_ECEF(1); y = R_ECEF(2); z = R_ECEF(3); % WGS84 Ellipsoid parameters a = 6378137.0; % Semimajor axis length of the Ellipsoid (in meters) b = 6356752.3142; % Semiminor axis length of the Ellipsoid (in meters) f = (ab)/a; % Flatness of the Ellipsoid e = sqrt(f*(2f)); % Eccentricity of the Ellipsoid p = sqrt(x^2 + y^2); h = 0; % Initialization N = a; % Initialization lambda = 0; % Initialization epsilon = 1; % Initialization % Iteration while epsilon > 1e8 lambda_old = lambda; sine_lambda = z/(N*(1e^2) + h); lambda = atan((z + e^2*N*sine_lambda)/p); N = a/sqrt(1 (e*sin(lambda))^2); h = (p/cos(lambda)) N; epsilon = lambda lambda_old; end geod_lat = lambda; geod_long = atan2(y,x); alt = h; function [g_Force, g_Moment] = Gravity_function(R_I, C_GI, Mass) % % This function calculates/outputs the Gravity Force Vector and the Gravity % Force Moment (about the Instantaneous Center of Mass of the Launch % Vehicle) acting on the Launch Vehicle due to the gravitational attraction % of the Earth. The Force and Moment vectors calculated in this function % are coordinatized in the Inertial frame. % % The input parameters (to be passed) to this function are the position of % the Launch Vehicle coordinatized in the Inertial frame, the % (instantaneous) transformation matrix from the Earth Centered Inertial % frame to the ECEF frame (C_GI) and the instantaneous Mass of the % Launch Vehicle. % Position vector of the Launch Vehicle (coordinatized in the ECEF frame) R_ECEF = C_GI*R_I; % Acceleration due to gravity acting on the Launch Vehicle (coordinatized % in the ECEF Frame) g_ECEF = gravity(R_ECEF); % Acceleration due to gravity acting on the Launch Vehicle (coordinatized % in the Inertial Frame) g_I = transpose(C_GI)*g_ECEF; PAGE 111 111 % Gravity Force acting on the Launch Vehicle (coordinatized in the Inertial % Frame) g_Force = Mass*g_I; % Gravity Force Moment (about the Instantaneous Center of Mass of the % Launch Vehicle) acting on the Launch Vehicle (coordinatized in the % Inertial Frame) g_Moment = [0;0;0]; function [g_e] = gravity(R_g) % This function computes the Gravity vector in the ECEF frame. This % function is based on WGS84EGM96 Gravity model. % % The input to this function is the position vector of the point of % interest in the ECEF frame. (3 X 1 vector) load WGS84EGM96normalized % The coefficients Cnm, Snm, Jn in the file WGS84EGM96 are Normalized % Gravitational coefficients. C = nC; S = nS; J = nJ; n_max = length(J); x = R_g(1); y = R_g(2); z = R_g(3); r = norm(R_g); R_g_unitvector = R_g/r; mu_earth = 398600.4418*10^9; % Standard Gravitational Parameter for Earth; a = 6378137.0; s = x/r; t = y/r; u = z/r; % % rho_n calculation rho_0 = mu_earth/r; rho = a/r; for i = 1:(1+n_max) if i == 1 rho_n(i) = rho*rho_0; else rho_n(i) = rho*rho_n(i1); end end %  PAGE 112 112 % Anm_bar calculations A0n_bar(2) = 1/2*5^(1/2)*(3*u^21); A0n_bar(3) = 1/2*7^(1/2)*u*(5*u^23); Anm_bar = zeros(n_max,n_max+1); Anm_bar(2,1) = 15^(1/2)*u; Anm_bar(2,2) = 1/2*15^(1/2); Anm_bar(3,1) = 1/4*42^(1/2)*(5*u^21); Anm_bar(3,2) = 5.1235*u; Anm_bar(3,3) = 2.0917; for n = 4:n_max A0n_bar(n) = (1/n)*(u*sqrt((2*n+1)*(2*n1))*A0n_bar(n1) (n1)*sqrt((2*n+1)/(2*n3))*A0n_bar(n2)); end for n = 4:n_max Anm_bar(n,n) = sqrt((2*n+1)/(2*n))*Anm_bar(n1,n1); end for n = 4:n_max for m = n1:1:1 Anm_bar(n,m) = 2*(m+1)*u*sqrt(1/((n+m+1)*(nm)))*Anm_bar(n,m+1) + (u^21)*sqrt((n+m+2)*(nm1)/((n+m+1)*(nm)))*Anm_bar(n,m+2); end end % % a1_bar calculation a1_bar = 0; for n = 2:n_max a1_bar_1 = 0; for m = 0:n if m == 0 a1_bar_2 = 0; else a1_bar_2 = Anm_bar(n,m)*m*Enm_bar(n,m,s,t,C,S,J); end a1_bar_1 = a1_bar_1 + a1_bar_2; end a1_bar = a1_bar + a1_bar_1*rho_n(n+1)/a; end % % a2_bar calculation a2_bar = 0; for n = 2:n_max a2_bar_1 = 0; for m = 0:n if m == 0 a2_bar_2 = 0; PAGE 113 113 else a2_bar_2 = Anm_bar(n,m)*m*Fnm_bar(n,m,s,t,C,S,J); end a2_bar_1 = a2_bar_1 + a2_bar_2; end a2_bar = a2_bar + a2_bar_1*rho_n(n+1)/a; end % % a3_bar calculation a3_bar = 0; for n = 2:n_max a3_bar_1 = 0; for m = 0:n if m == 0 a3_bar_2 = A0n_bar(n)*Dnm_bar(n,m,s,t,C,S,J); else a3_bar_2 = Anm_bar(n,m+1)*Dnm_bar(n,m,s,t,C,S,J); end a3_bar_1 = a3_bar_1 + a3_bar_2; end a3_bar = a3_bar + a3_bar_1*rho_n(n+1)/a; end % % a4_bar calculation a4_bar = 0; for n = 2:n_max a4_bar_1 = 0; for m = 0:n if m == 0 a4_bar_2 = A0n_bar(n)*Dnm_bar(n,m,s,t,C,S,J); else a4_bar_2 = Anm_bar(n,m)*Dnm_bar(n,m,s,t,C,S,J); end a4_bar_1 = a4_bar_1 + a4_bar_2; end a4_bar = a4_bar + a4_bar_1*(n+1)*rho_n(n+1)/a; end a4_bar = rho_0/r a4_bar s*a1_bar t*a2_bar u*a3_bar; % % Acceleration due to gravity in ECEF g_e = [a1_bar; a2_bar; a3_bar] + a4_bar*R_g_unitvector; end % % Dnm_bar subfunction function sol = Dnm_bar(n,m,s,t,C,S,J) rmst = r_m(m,s,t); PAGE 114 114 imst = i_m(m,s,t); if m == 0 sol = J(n)*rmst; else sol = C(n,m)*rmst + S(n,m)*imst; end end % % Enm_bar subfunction function sol = Enm_bar(n,m,s,t,C,S,J) rm1st = r_m(m1,s,t); im1st = i_m(m1,s,t); if m == 0 sol = J(n)*rm1st; else sol = C(n,m)*rm1st + S(n,m)*im1st; end end % % Fnm_bar subfunction function sol = Fnm_bar(n,m,s,t,C,S,J) rm1st = r_m(m1,s,t); im1st = i_m(m1,s,t); if m == 0 sol = J(n)*im1st; else sol = S(n,m)*rm1st C(n,m)*im1st; end end % % r_m subfunction function sol = r_m(m,s,t) sol = real((s+i*t)^m); end % % i_m subfunction PAGE 115 115 function sol = i_m(m,s,t) sol = imag((s+i*t)^m); end %  PAGE 116 116 APPENDIX B SIMULATION CONFIGURATION Figure B1. The DELTA II Launch vehicle geometry Figure B2. Strapon booster geometry PAGE 117 117 Instantaneous Center of Mass of the Rocket From the geometry of the rocket and strapon boosters, instantaneous center of mass of the rocket is calculated Strapon booster parameters Lns = Length of the nosecone Xbs = Length from the nose tip to fin root leading edge Lc = Length from the base of the nosecone to the base of the strapon booster Ys= Length from the base of the nosec one to the top of the solidmotor Hs= Length of the solidmotor R3s= Outer radius of the cylindr ical strapon booster frame R2s= Inner radius of the cylindr ical strapon booster frame R1s = Inner radius of th e solidpropellant shell CRs = Fin root chord XRs= Length from fin root leading edge to fin ti p leading edge parallel to the strapon booster body CTs= Fin tip chord Lfs = Length of fin mid chord line Ss= Span of one fin tfs = Thickness of the fin Nfs = Number of fins 1 s = Density of the material of the strapon p = Density of the solid propellant Thrusts= Thrust of the solid motor Isp_s= Specific impulse of solid motor in sec Rocket Parameters Ln = Length of the nosecone Ls = Length from the nose tip of the rocket to the nose tip of the strapon booster L1 = Length of the first stage of the rocket L2 = Length of the second stage of the rocket L3 = Length of the third stage of the rocket R1 = Outer radius of the cylindrical rocket frame R2 = Inner radius of the cylindrical shell of the 1st stage rocket frame R3 = Inner radius of the cylindrical shell of the 2nd stage rocket frame R4 = Inner radius of the cylindrical shell of the 3rd stage rocket frame R5 = Inner radius of the solid propellant shell of the 3rd stage M1pi = Initial mass of the propellant of the first stage M2pi = Initial mass of the prope llant of the second stage M4 = Mass of the payload (parabolic nosecone) Thrust1 = Thrust of the 1st stage PAGE 118 118 Isp_1 = Specific impulse of the 1st stage in sec Thrust2 = Thrust of the 2nd stage Isp_2 = Specific impulse of the 2nd stage in sec Thrust3 = Thrust of the 3rd stage Isp_3 = Specific impulse of the 3rd stage in sec 1 = Density of the material of the rocket frame Ns = Number of strapon boosters Re = Rocket nozzle exit radius Instantaneous mass of the rocket Mr = M1 + M2 + M3 + M4 + NsMs Initial mass of the rocket = Mri = M1i + M2i + M3i + M4 + NsMsi Instantaneous center of mass of the rocket w.r.t intial center of mass 23112344icsss c r X MxMxMxMxNMX r M where ic X = Initial Center of mass of the rocket 11223344 s siics iii ic riX M xMxMxMxNMX M Center of mass of first stage = x1 = (L1/2 + L2 + L3 + Ln) Initial mass of the first stage M1i = M1_1 M1_2 + M1pi Instantaneous mass of the first stage = M1 = M1i t1 dm1p/dt M1_1 = 1 R1 2 L1 M1_2 = 1 R2 2 L1 Rate of change of mass of the first stage liquid engine = dm1p/dt = Thrust1/(Isp_1 g) g = Acceleration due to gravity on the surface of the Earth. t1 = time elapsed in sec from the moment the first stage was ignited Center of mass of second stage = x2 = (L2/2 + L3 + Ln) Initial mass of the second stage = M2i = M2_1 M2_2 + M2pi Instantaneous mass of the second stage = M2 = M2i t2 dm2p/dt M2_1 = 1 R1 2 L2 M2_2 = 1 R3 2 L2 Rate of change of mass of th e second stage liquid engine = dm2p/dt = Thrust2/(Isp_2 g) g = Acceleration due to gravity on the surface of the Earth. t2 = time elapsed in sec from the moment the second stage was ignited Center of mass of third stage = x3 = (L3/2 + Ln) Initial mass of the third stage = M3i = M3_1 M3_2 + M3_3 M3_4 Instantaneous mass of the third stage = M3 = M3i t3 dm3p/dt M3_1 = 1 R1 2 L3 M3_2 = 1 R4 2 L3 M3_3 = p R4 2 L3 M3_4 = p R5 2 L3 Rate of change of mass of th e third stage solid motor = dm3p/dt = Thrust3/(Isp_3 g) g = Acceleration due to gravity on the surface of the Earth. t3 = time elapsed in sec from the moment the third stage was ignited PAGE 119 119 Center of mass of nose cone/payload = x4 = 2 Ln/3 Initial mass of the strapon booster Msi = M1s M2s +M3is + Ns Mfs + M5s Instantaneous mass of the strapon booster Ms = M1s M2s +M3ps + Ns Mfs + M5s Instantaneous center of ma ss of the strapon booster 55 1122334 s s sssspsss fsfs s s M xMxMxNMxMx X M Initial center of mass of the strapon booster 55 1122334 s s ssssisss fsfs ics si M xMxMxNMxMx X M x1s = (Ls + Lns + Lc/2) M1s = 1 s R3s 2 Lc x2s = (Ls + Lns + Lc/2) M2s = 1 s R2s 2 Lc x3s = (Ls + Lns + Ys + Hs/2) M3is = M3s M4s M3s = p R2s 2 Hs M4s = p R1s 2 Hs M3ps = M3is ts dms/dt Rate of change of mass of the solid motor = dms/dt = Thrusts/(Isp_s g) ts = time elapsed in sec from the moment the strapon booster was ignited g = Acceleration due to gravity on the surface of the Earth. x4s = Ls (L1s 2XRs 2/3 L2s (L1s L2s/3))/(2 L1s XRs L2s) Mass of fin Mfs= M1fs M2fs M3fs M1fs = 1 s L1s tfs Ss/8 M2fs = 0.5 1 s Ss tfs Ss/8 M3fs = 1 s L1s tfs Ss/8 where L1s = XRs + CTs, L2s = XRs + CTs CRs x5s = Ls 2 Lns/3 M5s = 500; Mass Moment of Inertia Calculation The Mass moment of inertia of the rocket (and booster) is calculated about the Cartesian axes through the instantaneous center of mass of the rocket. The major elements of the rocket are 1. First stage liquid motor 2. Second stage liquid motor PAGE 120 120 3. Third stage solid motor 4. Payload 5. Solid strapon boosters Figure B3. Elements of DELTA II Launch vehicle and Strapon Booster Strapon booster: The booster is assume d to be made of 1. Cylindrical shell (Frame) 2. Propellant shell (the solid pr opellant is assumed to be dist ributed as a cylindrical shell) 3. Parabolic nose cone 4. Fins Cylindrical shell: A cylindrical shell of outer radius R3s, inner radius R2s and height Lc is basically a solid cy linder of radius R3s and height Lc from which another cy linder of radius R2s and height Lc is removed (the bases and the centroi dal axes of both cy linders coincide). Outer radius of the cylindrical shell = R3s Inner radius of the cylindrical shell = R2s Height of the cylindrical shell = Lc Density of the material of the frame = 1 s PAGE 121 121 M1s = 1 s R3s 2 Lc M2s = 1 s R2s 2 Lc Figure B4. Cylindrical shell Mass moments of inertias of the two cyli nders about their center of masses are (Ixxs)1 = 0.5 M1s R3s 2 (Iyys)1 = (1/12) M1s (3 R3s 2 + Lc 2) (Izzs)1 = (1/12) M1s (3 R3s 2 + Lc 2) (Ixxs)2 = 0.5 M2s R2s 2 (Iyys)2 = (1/12) M2s (3 R2s 2 + Lc 2) (Izzs)2 = (1/12) M2s (3 R2s 2 + Lc 2) Mass moments of inertias of the tw o cylinders about the center of mass of the solid booster are (Ixcxcs)1 = (Ixxs)1 (Iycycs)1 = (Iyys)1 + M1s ((Lc/2 + Lns) Ls Xics)2 (Izczcs)1 = (Izzs)1 + M1s ((Lc/2 + Lns) Ls Xics)2 (Ixcxcs)2 = (Ixxs)2 (Iycycs)2 = (Iyys)2 + M2s ((Lc/2 + Lns) Ls Xics)2 (Izczcs)2 = (Izzs)2 + M2s ((Lc/2 + Lns) Ls Xics)2 Propellant shell: A cylindrical shell of outer radius R2s, inner radius R1s and height Hs is basically a solid cy linder of radius R2s and height Hs from which another cy linder of radius R1s and height Hs is removed (the bases and the centroidal axes of both cylinders coincide). It can be seen that as the propellant burns, the radius of th e inner cylinder R1s increases. When all the propellant is burnt, R1s = R2s. Outer radius of the propellant shell = R2s PAGE 122 122 Initial inner radius of the propellant shell = R1si Height of the propellant shell = Hs Density of the solid propellant = p Figure B5. Propellant shell M3s = p R2s 2 Hs M4si = p R1si 2 Hs Initial mass of the strapon booster propellant = Mspi = M3s M4si Instantaneous mass of the st rapon booster propellant = Msp = M3s M4s where M4s = M4si + ts dms/dt Rate of change of mass of the solid motor = dms/dt = Thrusts / (Isp_s g) ts = time elapsed in sec from the moment the strapon booster was ignited g = Acceleration due to gravity on the surface of the Earth. Instantaneous inner radius of the propellant shell = R1s = 4s sM Hp Mass moments of inertias of the two cyli nders about their center of masses are (Ixxs)3 = 0.5 M3s R2s 2 (Iyys)3 = (1/12) M3s (3 R2s 2 + Hs 2) (Izzs)3 = (1/12) M3s (3 R2s 2 + Hs 2) (Ixxs)4 = 0.5 M4s R1s 2 (Iyys)4 = (1/12) M4s (3 R21s 2 + Hs 2) (Izzs)4 = (1/12) M4s (3 R1s 2 + Hs 2) Mass moments of inertias of the tw o cylinders about the center of mass of the solid booster are (Ixcxcs)3 = (Ixxs)3 (Iycycs)3 = (Iyys)3 + M3s ((Hs/2 + Ys + Lns) Ls Xics)2 (Izczcs)3 = (Izzs)3 + M3s ((Hs/2 + Ys + Lns) Ls Xics)2 (Ixcxcs)4 = (Ixxs)4 (Iycycs)4 = (Iyys)4 + M4s ((Hs/2 + Ys + Lns) Ls Xics)2 PAGE 123 123 (Izczcs)4 = (Izzs)4 + M4s ((Hs/2 + Ys + Lns) Ls Xics)2 Parabolic nose cone: The nose cone is assumed to be parabolic as shown in Fig. B6. Height of the nose cone = Lns Radius of the nose cone = R3s Mass of the nose cone = M5s Mass moments of inertia of the parabolic nose cone about it s center of mass are (Ixxs)5 = (1/3) M5s R3s 2 (Iyys)5 = (1/6) M5s (R3s 2 + Lns 2/3) (Izzs)5 = (1/6) M5s (R3s 2 + Lns 2/3) Mass moments of inertia of the parabolic nose co ne about the center of mass of the solid booster are: (Ixcxcs)5 = (Ixxs)5 (Iycycs)5 = (Iyys)5 + M5s (2 Lns/3 Ls Xics)2 (Izczcs)5 = (Izzs)5 + M5s (2 Lns/3 Ls Xics)2 Figure B6. Parabolic nose cone Fins: Each strapon booster is assumed to have four fins. The (geometry of the) fin is obtained by removing two triangular slabs from a rect angular slab as shown in figure B7. M1fs = 1 s L1s tfs Ss/8 M2fs = 0.5 1 s Ss tfs Ss/8 M3fs = 1 s L1s tfs Ss/8 where L1s = XRs + CTs L2s = XRs + CTs CRs The center of mass of the fin is given by XCm_s = (L1s 2XRs 2/3 L2s (L1s L2s/3))/(2 L1s XRs L2s) YCm_s = (L1s Ss2 XRs Ss/3 L2s Ss/3)/(2 L1s XRs L2s) (Ixx)1fs = (1/12) M1fs (Ss 2 + tfs 2) (Iyy)1fs = (1/12) M1fs (L1s 2 + tfs 2) (Izz)1fs = (1/12) M1fs (Ss 2 + L1s 2) (Ixx)2fs = M2fs (Ss 2/18 + tfs 2/12) PAGE 124 124 (Iyy)2fs = M2fs ( XRs 2/18 + tfs 2/12) (Izz)2fs = (1/18) M2fs (Ss 2 + L2s 2) (Ixx)3fs = M3fs (Ss 2/18 + tfs 2/12) (Iyy)3fs = M3fs ( L2s 2/18 + tfs 2/12) (Izz)3fs = (1/18) M3fs (Ss 2 + L2s 2) Figure B7. Fins The mass moments of inertia of the fi n about its centroida l axes are given by (Ixx)fs = (Ixcxc)1fs (Ixcxc)2fs (Ixcxc)3fs (Iyy)fs = (Iycyc)1fs (Iycyc)2fs (Iycyc)3fs (Izz)fs = (Izczc)1fs (Izczc)2fs (Izczc)3fs where (Ixcxc)1fs = (Ixx)1fs + M1fs (Ss/2 Ycm_s)2 (Iycyc)1fs = (Iyy)1fs + M1fs (L1s/2 Xcm_s)2 (Izczc)1fs = (Izz)1fs + M1fs ((L1s/2 Xcm_s)2 + (Ss/2 Ycm_s)2) (Ixcxc)2fs = (Ixx)2fs + M2fs (2 Ss/3 Ycm_s)2 (Iycyc)2fs = (Iyy)2fs + M2fs (XRs/3 Xcm_s)2 (Izczc)2fs = (Izz)2fs + M2fs ((2 Ss/3 Ycm_s)2 + (XRs/3 Xcm_s)2) (Ixcxc)3fs = (Ixx)3fs + M3fs (Ss/3 Ycm_s)2 (Iycyc)3fs = (Iyy)3fs + M3fs (L1s L2s/3 Xcm_s)2 (Izczc)3fs = (Izz)3fs + M3fs ((Ss/3 Ycm_s)2 + (L1s L2s/3 Xcm_s)2) Mass of fin Mfs= M1fs M2fs M3fs (Ixcxcs)6 = (Ixcxcs)8 = (Ixx)fs + Mfs (YCm_s + R3s)2 (Iycycs)6 = (Iycycs)8 = (Iyy)fs + Mfs ((XCm_s + Xbs) Ls Xics)2 (Izczcs)6 = (Izczcs)8 = (Izz)fs + Mfs ((YCm_s + R3s)2 + ((XCm_s + Xbs) Ls Xics)2) (Ixcxcs)7 = (Ixcxcs)9 = (Ixx)fs + Mfs (YCm_s + R3s)2 (Iycycs)7 = (Iycycs)9 = (Izz)fs + Mfs ((YCm_s + R3s)2 + ((XCm_s + Xbs) Ls Xics)2) (Izczcs)7 = (Izczcs)9 = (Iyy)fs + Mfs ((XCm_s + Xbs) Ls Xics)2 The Mass moments of inertia of the booster about its centroidal axes are (Ixx)s = (Ixcxcs)1 (Ixcxcs)2 + (Ixcxcs)3 (Ixcxcs)4 + (Ixcxcs)5 + (Ixcxcs)6 + (Ixcxcs)7 + (Ixcxcs)8 + (Ixcxcs)9 (Iyy)s = (Iycycs)1 (Iycycs)2 + (Iycycs)3 (Iycycs)4 + (Iycycs)5 + (Iycycs)6 + (Iycycs)7 + (Iycycs)8 PAGE 125 125 + (Iycycs)9 (Izz)s = (Izczcs)1 (Izczcs)2 + (Izczcs)3 (Izczcs)4 + (Izczcs)5 + (Izczcs)6 + (Izczcs)7 + (Izczcs)8 + (Izczcs)9 First stage: The first stage is assumed to be made of a cylindrical shell (frame) and liquid propellant. The liquid propellant is assumed to be distributed as a solid cylinder whose density decreases with time (i.e., as the propellant burns). Figure B8. First stage Cylindrical shell: A cylindrical shell of outer radius R1, inner radius R2 and height L1 is basically a solid cy linder of radius R1 and height L1 from which another cylinder of radius R2 and height L1 is removed (the bases and the centroi dal axes of both cylinders coincide). Outer radius of the cylindrical shell = R1 Inner radius of the cylindrical shell = R2 Height of the cylindrical shell = L1 Density of the material of the frame = 1 M1_1 = 1 R1 2 L1 M1_2 = 1 R2 2 L1 Mass moments of inertias of the two cyli nders about their center of masses are (Ixx)1_1 = 0.5 M1_1 R1 2 (Iyy)1_1 = (1/12) M1_1 (3 R1 2 + L1 2) (Izz)1_1 = (1/12) M1_1 (3 R1 2 + L1 2) (Ixx)1_2 = 0.5 M1_2 R2 2 (Iyy)1_2 = (1/12) M1_2 (3 R2 2 + L1 2) (Izz)1_2 = (1/12) M1_2 (3 R2 2 + L1 2) PAGE 126 126 Propellant cylinder: The liquid propellant is assumed to be distributed as a solid cylinder whose density decreases with time (i.e., as the propellant burns). Radius of the cylinder = R2 Instantaneous mass of the unburnt first stage liquid propellant M1_3 = M1pi t1 dm1p/dt Where M1pi = Initial mass of the first stage propellant t1 = time elapsed in sec from the moment the first stage is ignited dm1p/dt = Mass flow rate of the first stage propellant Rate of change of mass of the first stage liquid engine = dm1p/dt = Thrust1 / (Isp_1 g) g = Acceleration due to gravity on the surface of the Earth. (Ixx)1_3 = 0.5 M1_3 R2 2 (Iyy)1_3 = (1/12) M1_3 (3 R2 2 + L1 2) (Izz)1_3 = (1/12) M1_3 (3 R2 2 + L1 2) The mass moments of inertia of the first stage about its centroidal axes are given by (Ixx)1 = (Ixx)1_1 (Ixx)1_2 + (Ixx)1_3 (Iyy)1 = (Iyy)1_1 (Iyy)1_2 + (Iyy)1_3 (Izz)1 = (Izz)1_1 (Izz)1_2 + (Izz)1_3 Mass of the first stage =M1 = M1_1 M1_2 + M1_3 The mass moments of inertia of the first stage about the instantaneous centroidal axes of the rocket are given by: (Ixcxc)1 = (Ixx)1 (Iycyc)1 = (Iyy)1 + M1 (L1/2 + L2 + L3 + Ln Xic + rc)2 (Izczc)1 = (Izz)1 + M1 (L1/2 + L2 + L3 + Ln Xic + rc)2 Second stage: The second stage is also assumed to be made of a cylindrical shell (frame) and liquid propellant. The liquid propellant is assumed to be distributed as a solid cylinder whose density decreases with time (i.e., as the propellant burns). Cylindrical shell: A cylindrical shell of outer radius R1, inner radius R3 and height L2 is basically a solid cy linder of radius R1 and height L2 from which another cylinder of radius R3 and height L2 is removed (the bases and the centroi dal axes of both cylinders coincide). Outer radius of the cylindrical shell = R1 Inner radius of the cylindrical shell = R3 Height of the cylindrical shell = L2 Density of the material of the frame = 1 M2_1 = 1 R1 2 L2 M2_2 = 1 R3 2 L2 Mass moments of inertias of the two cyli nders about their center of masses are (Ixx)2_1 = 0.5 M2_1 R1 2 (Iyy)2_1 = (1/12) M2_1 (3 R1 2 + L2 2) PAGE 127 127 Figure B9. Second stage (Izz)2_1 = (1/12) M2_1 (3 R1 2 + L2 2) (Ixx)2_2 = 0.5 M2_2 R3 2 (Iyy)2_2 = (1/12) M2_2 (3 R3 2 + L2 2) (Izz) 2_2 = (1/12) M2_2 (3 R3 2 + L2 2) Propellant cylinder: The liquid propellant is assumed to be distributed as a solid cylinder whose density decreases with time (i.e., as the propellant burns). Radius of the cylinder = R3 Instantaneous mass of the unburnt second stage liquid propellant M2_3 = M2pi t2 dm2p/dt Where M2pi = Initial mass of the second stage liquid propellant t2 = time elapsed in sec from the moment the second stage is ignited dm2p/dt = Mass flow rate of the second stage propellant Rate of change of mass of the first stage liquid engine = dm2p/dt = Thrust2/(Isp_2 g) g = Acceleration due to gravity on the surface of the Earth. (Ixx)2_3 = 0.5 M2_3 R3 2 (Iyy)2_3 = (1/12) M2_3 (3 R3 2 + L2 2) (Izz)2_3 = (1/12) M2_3 (3 R3 2 + L2 2) The mass moments of inertia of the second st age about its centroida l axes are given by: (Ixx)2 = (Ixx)2_1 (Ixx)2_2 + (Ixx)2_3 (Iyy)2 = (Iyy)2_1 (Iyy)2_2 + (Iyy)2_3 (Izz)2 = (Izz)2_1 (Izz)2_2 + (Izz)2_3 Mass of the second stage =M2 = M2_1 M2_2 + M2_3 The mass moments of inertia of the second stag e about the instantaneous centroidal axes of the rocket are given by: (Ixcxc)2 = (Ixx)2 (Iycyc)2 = (Iyy)2 + M2 (L2/2 + L3 + Ln Xic + rc)2 (Izczc)2 = (Izz)2 + M2 (L2/2 + L3 + Ln Xic + rc)2 PAGE 128 128 Third stage: The third stage is assumed to be made of a cylindrical shell and cylindrical propellant shell. Cylindrical shell: A cylindrical shell of outer radius R1, inner radius R4 and height L3 is basically a solid cy linder of radius R1 and height L3 from which another cylinder of radius R4 and height L3 is removed (the bases and the centroi dal axes of both cylinders coincide). Figure B10. Third stage Outer radius of the cylindrical shell = R1 Inner radius of the cylindrical shell = R4 Height of the cylindrical shell = L3 Density of the material of the frame = 1 M3_1 = 1 R1 2 L3 M3_2 = 1 R4 2 L3 Mass moments of inertias of the two cyli nders about their center of masses are (Ixx)3_1 = 0.5 M3_1 R1 2 (Iyy)3_1 = (1/12) M3_1 (3 R1 2 + L3 2) (Izz)3_1 = (1/12) M3_1 (3 R1 2 + L3 2) (Ixx)3_2 = 0.5 M3_2 R4 2 (Iyy)3_2 = (1/12) M3_2 (3 R4 2 + L3 2) (Izz)3_2 = (1/12) M3_2 (3 R4 2 + L3 2) Propellant shell: A cylindrical shell of outer radius R4, inner radius R5 and height L3 is basically a solid cy linder of radius R4 and height L3 from which another cylinder of radius R5 and height L3 is removed (the bases and the centroi dal axes of both cylinders coincide). Outer radius of the propellant shell = R4 PAGE 129 129 Initial inner radius of the propellant shell = R5si Height of the Propellant shell = L3 Density of the solid propellant = p M3_3 = p R4 2 L3 M3_4i = p R5i 2 L3 Initial mass of the thirdstage solid propellant = M3pi = M3_3 M3_4i Instantaneous mass of the thirdstage propellant = M3p = M3_3 M3_4 where M3_4 = M3_4i + t3 dm3/dt Rate of change of mass of th e third stage solid motor = dm3/dt = Thrust3 / (Isp_3 g) t3 = time elapsed in sec from the moment the third stage was ignited g = Acceleration due to gravity on the surface of the Earth. Instantaneous inner radius of the propellant shell = R5 = 3_4 3M Lp Mass moments of inertias of the two cyli nders about their center of masses are (Ixx)3_3 = 0.5 M3_3 R4 2 (Iyy)3_3 = (1/12) M3_3 (3 R4 2 + L3 2) (Izz)3_3 = (1/12) M3_3 (3 R4 2 + L3 2) (Ixx)3_4 = 0.5 M3_4 R5 2 (Iyy)3_4 = (1/12) M3_4 (3 R5 2 + L3 2) (Izz)3_4 = (1/12) M3_4 (3 R5 2 + L3 2) The mass moments of inertia of the third st age about its centroidal axes are given by: (Ixx)3 = (Ixx)3_1 (Ixx)3_2 + (Ixx)3_3 (Ixx)3_4 (Iyy)3 = (Iyy)3_1 (Iyy)3_2 + (Iyy)3_3 (Iyy)3_4 (Izz)3 = (Izz)3_1 (Izz)3_2 + (Izz)3_3 (Izz)3_4 Mass of the third stage = M3 = M3_1 M3_2 + M3_3 M3_4 The mass moments of inertia of the third stage about the instantaneous centroidal axes of the rocket are given by: (Ixcxc)3 = (Ixx)3 (Iycyc)3 = (Iyy)3 + M3 (L3/2 + Ln Xic + rc)2 (Izczc)3 = (Izz)3 + M3 (L3/2 + Ln Xic + rc)2 Payload: The payload is assumed to be a parabolic (nose) cone Mass of the nosecone (pay load) = M4 Radius of the nosecone = R1 Length of the nosecone = Ln Mass moments of inertia of the nosecone about its centroida l axes are given by (Ixx)p = M4 R1 2/3 (Iyy)p = (M4/6) (R1 2 + Ln 2/3) (Izz)p = (M4/6) (R1 2 + Ln 2/3) Mass moments of inertia of the nos econe about the instantaneous centroidal axes of the rocket are given by: (Ixcxc)p = (Ixx)p (Iycyc)p = (Iyy)p + M4 (Ln/2 Xic + rc)2 PAGE 130 130 (Izczc)p = (Izz)p + M4 (Ln/2 Xic + rc)2 Figure B11. Payload The mass moments of inertias of Ns boosters about the instantane ous centroidal axes of the rocket are given by: (Ixcxc)s = Ns ((Ixx)s + Ms (R3s+Lcr+R1)2) (Iycyc)s = Ns ((Iyy)s + Ms (Xics + Ls Xic + rc)2) + Ms (R3s+Lcr+R1)2 Yres (Izczc)s = Ns ((Izz)s + Ms (Xics + Ls Xic + rc)2) + Ms (R3s+Lcr+R1)2 Zres Figure B12. Strapon boos ters around the Rocket where Yres = 1 2 0sin(*)sN ii Zres = 1 2 0cos(*)sN ii 360 s N The mass moments of inertia of the rocket a bout its instantaneous cen troidal axes are given by: Ixx = (Ixcxc)1 + (Ixcxc)2 + (Ixcxc)3 + (Ixcxc)p + (Ixcxc)s Iyy = (Iycyc)1 + (Iycyc)2 + (Iycyc)3 + (Iycyc)p + (Iycyc)s Izz = (Izczc)1 + (Izczc)2 + (Izczc)3 + (Izczc)p + (Izczc)s PAGE 131 131 The products of inertias of the rocket are zero because of the symmetry (i.e., the Cartesian axes through the instantaneous center of mass is the "principal axes of the rocket). The inertia tensor = 00 00 00xx yy zzI II I PAGE 132 132 LIST OF REFERENCES 1. J. Apt, T. Cochran, A. Eng, K. Florig, J. Lyngdal and A.H. Barber, Launching Safely in the 21st Century, Marion, Iowa: National Asso ciation of Rocketry, 2005. [Online]. Available: www.nar.org/pdf/launchsafe.pdf [Accessed May 15, 2007]. 2. J. Barrowman, Technical Information Report 33 Calculating the Center of Pressure of a Model Rocket, Phoenix, Arizona: Centur i Engineering Company, 1970. [Online]. Available: http://www.rockets4schools.org/ education/Calculating_Cp.pdf [Accessed May 15, 2007]. 3. J. Barrowman, The practical calculation of the aerodynamic characteristics of slender finned vehicles, M.S. thesis, The Catholic Un iversity of America, Washington, DC, USA, 1967. 4. J. Barrowman and J. A. Barrowman, The Theoretical Prediction of the Center of Pressure, Colorado Springs, Colorado: Apogee Components Inc., 2005. [Online]. Available: www.apogeerockets.com/education/ downloads/Barrowman_report.pdf [Accessed May 15, 2007]. 5. T. Benson, Rocket Stability, Cleveland, Ohio: Glenn Research Center, 2006. [Online]. Available: http://exploration.grc.nasa.gov/ education/rocket/rktstab.html [Accessed May 15, 2007]. 6. A. Bethencourt, J. Wang, C. Rizos, and A. H. W. Kearsley, Using personal computers in spherical harmonic synthesis of high de gree Earth Geopotential Models, Cairns, Australia: International Association of Geodesy, 2005. [Online]. Available: www.gmat.unsw.edu.au/snap/publica tions/bethencourt_etal2005a.pdf [Accessed May 15, 2007]. 7. Boeing Company, DELTA II Payload Planners Guide. Huntington Beach, California: The Boeing Company, 1996. 8. J. C. Chen, C. C. Wang, D. Taggart, and E. Ditata, Modeling an offnominal launch vehicle trajectory for range safety link analysis, in Proceedings of the 2003 IEEE Aerospace Conference, Vol. 7, March 2003, pp. 34193426. 9. J. W. Cornelisse, H. F. R. Schoyer, and K. F. Wakker, Rocket propulsion and spaceflight dynamics. London: Pitman, 1979. 10. E. Denson, L. M. Valencia, R. Sakahara, J. C. Simpson, D. E. Whiteman, S. Bundick, D. Wampler, and R. B. Birr, Spacebased Tele metry and Range Safety Flight Demonstration #1 Final Report, Kennedy Space Center, FL Rep. KSCYA6400 Revision Basic, 2004. 11. European Cooperation for Space Standardiza tion, System engineering coordinate systems, Noordwijk, Netherlands: ESA Publi cation Division, 2007. [Online]. Available: www.ecss.nl/.../dispatch.cgi/home/ showFile/100317/d20070507095706/No/ecsse1009ADraft1.1(3April2007).pdf [Accessed June 03, 2007]. PAGE 133 133 12. J. A. Farrell and M. Barth, The Global Positioning System & Inertial Navigation, 1st ed. McGrawHill Profession al Publication, 1998. 13. E. L. Fleeman, Maximizing Missile Flight Performance, Atlanta, Georgia: E. L. Fleeman, 2002. [Online]. Available: www.dbf.gatech.edu/performance.ppt [Accessed May 15, 2007]. 14. C. R. Kehler, Eastern and Western Range 1271, Range Safety Requirements. Patrick Air Force Base, Florida: Ra nge Safety Office, 1999. 15. Y. Y. Krikorian, D. A. Taggart, C. C. Wang, R. Kuma r, C. Chen, S.K. Do, D. L. Emmons, J. Hant, A. Mathur, M. Cutler, and N. El yashar, Dynamic link analysis tool for a telemetry downlink system, in IEEE Aerospace Conference Proceedings, Vol. 6, 613 March 2004, pp. 39944003 16. J. B. Lundberg and B. E. Schutz, Recursion formulas of Legendre functions for use with nonsingular geopotential models, Journal of Guidance, Control, and Dynamics, vol.11, no.1, pp. 3138, 1988. 17. M. Madden, Gravity Modeling for Vari able Fidelity Environments, in AIAA Modeling and Simulation Technologies Conference and Exhibit, 2006, pp. 10811096. 18. M. J. Muolo and R. A. Hand, Space Handbook A War Fighter's Guide to Space, Vol. 1. Maxwell Air Force Base, Alabama: Air Univer sity Press, 1993. [Online] Available: http://www.au.af.mil/au/awc/awcgate/au18/au180001.htm [Accessed: May 15, 2007]. 19. NASA, F15B SBRDC/ECANS, Edwards, Cali fornia: Dryden Flight Research Center, 2007. [Online]. Available: http://www1.dfrc.nasa.gov/ Gallery/Photo/F15B_SBRDCECANS/index.html [Accessed June 10, 2007]. 20. NASA, Emerging Technology Space Based Te lemetry and Range Safety 2006, Merritt Island, Florida: Kennedy Space Center, 2006. [Online]. Available: http://kscsma.ksc.nasa.gov/Range_Safety/ Annual_Report/2006/Printpages/ETSpaceBasedTelemetryimages.pdf [Accessed June 10, 2007]. 21. NASA, Flight Demonstration #2 (Technical Interchange Meeting), Space Based Range Demonstration & Certification, Kennedy Sp ace Center, Florida. Rep. STARS FD#2, 2006. 22. NASA, TDRSS Overview, Greenbelt, Maryla nd: Goddard Space Flight Center, 2005. [Online]. Available: http://msp.gsfc.nasa .gov/tdrss/oview.html [Accessed June 10, 2007]. 23. National Research Council, Streamlining Space Launch Range Safety. Washington, D.C.: National Academy Press, 2000. 24. F. Ohrtman and K. Roeder, WiFi Handbook: Building 802.11b wireless networks, 1st ed. McGrawHill Professional publication, 2003, pp.218 219. PAGE 134 134 25. S. Pines, Uniform Representation of the Gravitational Potential and its Derivatives, AIAA Journal, vol. 11, pp. 15081511, Nov. 1973. 26. L. C. Rabelo, J. Sepulveda, J. Compton, and R. Turner, Simulation of range safety for the NASA space shuttle, Aircraft Engineering & Aerospace Technology, vol. 78, no. 2, pp. 98, 2006. 27. G. Rosborough, Spherical Harmonic Representati on of the Gravity Field Potential, Palo Alto, California: cdeagle, 2006. [Online]. Available: http://www.cdeagle.com/pdf/gravity.pdf [Accessed May 15, 2007]. 28. F. Takahashi, Y. Takahashi, T. Kondo, and Y. Koyama, Very Long Baseline Interferometer. Amsterdam: IOS Press, 2000. 29. J. R. Wertz and W. J. Larson., Space Mission Analysis and Design, 3rd ed. Microcosm Press, October 1999, pp. 208. 30. D. E. Whiteman, L. M. Valencia, and R. B. Birr, SpaceBased Telemetry and Range Safety Project KuBand and KaBand Phas ed Array Antenna, NASA Dryden Flight Research Center, Edwards, California. Rep. H2603, NASA TM2005212872, 2005. 31. D. E. Whiteman, L. M. Valencia, and J. C. Simpson, SpaceBased Range Safety and Future Space Range Applications, NASA Dr yden Flight Research Center, Edwards, California. Rep. H2616, NASA TM2005213662, 2005. PAGE 135 135 BIOGRAPHICAL SKETCH Sharath Chandra Prodduturi was born in Hydera bad, India, in 1983. He graduated with a bachelors degree in mechanical engineering from Osmania University, India, in 2004. In 2005, he moved to Gainesville, Florid a, to pursue his Master of Sc ience degree in mechanical engineering at the University of Florida. 