Question? Leave a message!




Orbital Mechanics!

Orbital Mechanics!
Orbital Mechanics Space System Design, MAE 342, Princeton University Robert Stengel Conic section orbits Equations of motion Momentum and energy Kepler’s Equation Position and velocity in orbit Copyright 2016 by Robert Stengel. All rights reserved. For educational use only. http://www.princeton.edu/stengel/MAE342.html 1 Orbits 101 Satellites Escape and Capture (Comets, Meteorites) 2TwoBody Orbits are Conic Sections 3 Classical Orbital Elements Dimension and Time a : Semimajor axis e : Eccentricity t : Time of perigee passage p Orientation :Longitude of the Ascending/Descending Node i : Inclination of the Orbital Plane ": Argument of Perigee 4Orientation of an Elliptical Orbit i ir rs st t P Po oi in n"" o of f A Ar ri ie es s 5 Orbits 102 (2Body Problem) • e.g., – Sun and Earth or – Earth and Moon or – Earth and Satellite • Circular orbit: radius and velocity are constant • Low Earth orbit: 17,000 mph = 24,000 ft/s = 7.3 km/s • Supercircular velocities – Earth to Moon: 24,550 mph = 36,000 ft/s = 11.1 km/s – Escape: 25,000 mph = 36,600 ft/s = 11.3 km/s • Near escape velocity, small changes have huge influence on apogee 6nd Newton’s 2 Law Particle of fixed mass (also called a point mass) acted upon by a force changes velocity with acceleration proportional to and in direction of force Inertial reference frame Ratio of force to acceleration is the mass of the particle: F = m a d dv t ( ) mv t = m = ma t = F ( ) ( ) " v t ( ) f dt dt x x d f x f m v (t) = y y dt f F = = force vector y f v t ( ) z z f " z" " 7 Equations of Motion for a Particle nd Integrating the acceleration (Newton’s 2 Law) allows us to solve for the velocity of the particle f m a x x dv t 1 ( ) f m a =v(t)= F = = y y dt m f m a z z "" T T T dv t 1 ( ) v(T)= dt +v(0)= a(t)dt +v(0)= Fdt +v(0) 0 0 0 dt m 3 components of velocity v (T) a (t) v (0) f (t) m v (0) x x x x x T T v (T)= a (t)dt + v (0)= f (t) mdt + v (0) y y y y y '' 0 0 v (T) a (t) v (0) f (t) m v (0) z z z z z 8 """""Equations of Motion for a Particle Integrating the velocity allows us to solve for the position of the particle x t v (t) ( ) x dr t ( ) = r t = v t = y t = v (t) ( ) ( ) ( ) y dt z(t) v t ( ) z "" T T dr(t) r(T)= dt +r(0)= vdt +r(0) 0 0 dt 3 components of position v t x(T) ( ) x(0) x T v t y(T)= ( ) dt + y(0) y ' 0 z(T) v t z(0) ( ) z 9 """ Spherical Model of the Rotating Earth Spherical model of earth s surface, earthfixed (rotating) coordinates x cosL cos' o E E R = y = cosL sin' R E o E E z sinL o E "" E L : Latitude (from Equator), deg E : Longitude (from Prime Meridian), deg E R : Radius (from Earth's center), deg Earth's rotation rate, , is 15.04 deg/hr 10NonRotating (Inertial) Reference Frame for the Earth Celestial longitude, , measured C from First Point of Aries on the Celestial Sphere at Vernal Equinox = +" tt = +"t ( ) C E epoch E 11 Transformation Effects of Rotation Transformation from inertial frame, I, to Earth’s rotating frame, E ' x '' o cos"t sin"t 0 cos"t sin"t 0 ) )) R = R = y sin"t cos"t 0sin"t cos"t 0 E I o) )) ) 0 0 1 0 0 1 )) z (( o ( I Location of satellite, rotating and inertial frames "" cosL cos cosL cos E E E C '' r = cosL sin R+Altitude ; r = cosL sin R+Altitude ( ) ( ) E E' E C' E I '' sinL sinL E E Orbital calculations generally are made in an inertial frame of reference 12Gravity Force Between Two Point Masses, e.g., Earth and Moon Magnitude of gravitational attraction Gmm 1 2 F= 2 r "11 2 2 G : Gravitational constant = 6.6710 Nm /kg st 24 m : Mass of 1 body= 5.9810 kg for Earth 1 nd 22 m : Mass of 2 body= 7.3510 kg for Moon 2 r : Distance between centers of mass of m and m , m 1 2 13 Acceleration Due To Gravity Gmm 1 2 F =m a = 2 2 1on2 2 “Inversesquare r Law” Gm µ 1 1 a = 1on2 2 2 r r µ =Gm 1 1 st Gravitational parameter of 1 mass At Earth’s surface, acceleration due to gravity is 14 3 2 µ 3.9810 m s 2 E a g = = = 9.798m s g o 2 2 Earth R 6,378,137m ( ) surface 14Gravitational Force Vector of the Spherical Earth Force always directed toward the Earth’s center (+ x " µ r µ E I E y F =m =m (vector), as r = r g I I 2 3 ' r r r I I I z ), I (x, y, z) establishes the direction of the local vertical x y cosL cos' E I z r " I I cosL sin' = = E I 2 2 2 r x + y +z I I I I sinL 15 E " Equations of Motion for a Particle in an InverseSquareLaw Field nd Integrating the acceleration (Newton’s 2 Law) allows us to solve for the velocity of the particle (+ x " dv(t) 1 µ r µ E I E =v t = F = = y ( ) g 2 3 ' dt m r r r I I I z ), T T T dv t 1 ( ) v(T)= dt +v(0)= a(t)dt +v(0)= Fdt +v(0) 0 0 0 dt m 3 components of velocity 3 v (T) v (0) x r x x I T 3 v (T)='µdt+ v (0) y r y E y I ( 0 3 v (T) z r v (0) z I z 16 """Equations of Motion for a Particle in an InverseSquareLaw Field As before; Integrating the velocity allows us to solve for the position of the particle x t v (t) ( ) x dr t ( ) = r t = v t = y t = v (t) ( ) ( ) ( ) y dt z(t) v t ( ) z "" T T dr(t) r T = dt +r 0 = vdt +r 0 ( ) ( ) ( ) 0 0 dt 3 components of position v t x(T) ( ) x(0) x T v t y(T)= ( ) dt + y(0) y ' 0 z(T) v t z(0) ( ) z "" 17 " Dynamic Model with Inverse SquareLaw Gravity No aerodynamic or thrust force Neglect motions in the z direction m m satellite Earth Dynamic Equations Example: Initial Conditions at Equator 3 v (t)=µ x (t) r (t) x E I I v 0 = 7.5,8,8.5 km/s ( ) 3 x v t =µ y t r t ( ) ( ) ( ) y E I I v 0 = 0 ( ) y x t = v t ( ) ( ) I x x 0 = 0 ( ) y t = v t ( ) ( ) I y y 0 = 6,378 km = R ( ) 2 2 where r t = x t + y t ( ) ( ) ( ) I I I 18Equatorial Orbits Calculated with InverseSquareLaw Model 19 Work “Work” is a scalar measure of change in energy With constant force, In one dimension W =F rr =F"r ( ) 12 2 1 In three dimensions T T W =F rr =F"r ( ) 12 2 1 With varying force, work is the integral " dx r r 2 2 ' T W = F dr = f dx+ f dy+ f dz , dr = dy ( ) 12 x y z ' r r 1 1 ' dz 20Conservative Force Iron Filings Assume that the 3D force Around Magnet field is a function of position F =F r ( ) Force Emanating The force field is conservative if from Source r r 2 1 T T F r dr+ F r dr =0 ( ) ( ) r r 1 2 … for any path between r and r and back 1 2 21 Gravitational Force is Gradient of a Potential, V(r) Gravity potential, V(r), is a function only of position µ" µ" E E F =m r = m V(r) g I 3 ( r"r r"r ' I 22Gravitational Force Field Gravitational force field µ E F =m r g I 3 r I Gravitational force field is conservative because r r 2 1 V r dr V r dr = ( ) ( ) I I "" rr r r 1 2 r r 2 2 µ µ E E m r dr + m r dr =0 3 I I 3 I I "" r r I I r r 1 1 … for any path between r and r and back 1 2 23 Potential Energy in Gravitational Force Field Potential energy, V or PE, is defined with respect to a reference point, r 0 PE r =V r =V"U ( ) ( ) ( ) 0 0 0 0 µ µ µ µ PEV(r )"V(r )=" m +V + m +V ="m +m 2 1 0 0 (( r r r r '' 2 1 2 1 24Kinetic Energy nd Apply Newton’s 2 Law to the definition of Work dr dv = v; dr= vdt F= m dt dt r t T t 2 2 2 dv " T 1 1 2 2 2 W = F dr = m vdt 12' = mv = m"v tv t ( ) ( ) 2 1 dt r t 1 1 2 2 t 1 t t 2 2 1 d 1 d T 2 1 = m v v dt = m v dt 2 2 ( ) ( ) = m"vvTTKE 2 dt 2 dt 2 1 2 1 t t 1 1 2 st nd Work = integral from 1 to 2 time T (= KE) is the kinetic energy of the point mass, m 25 Total Energy in PointMass Gravitational Field Potential energy of mass, µ m, depends only on the V =PE=m r gravitational force field Kinetic energy of mass, m, 1 2 depends only on the velocity TKE= mv 2 magnitude measured in an inertial frame of reference E =PE+KE Total energy, E, is the µ 1 2 =m + mv sum of the two: r 2 = Constant 26Interchange Between Potential and Kinetic Energy in a Conservative System E E =0 2 1 "" µ 1 µ 1 2 2 m + mvm + mv =0 2 1 '' r 2 r 2 2 1 P 2 " µ µ 1 1 " 2 2 m +m = mv mv 2 1' ' r r 2 2 2 1 PEPE =KEKE 2 1 2 1 P 1 27 Specific Energy… Energy per unit of the satellite’s mass P E =PE +KE 1 S S S 1 mµ 1 " 2 = + mv ' m r 2 µ 1 P 2 2 = + v r 2 28Angular Momentum of a Particle (Point Mass) h = rmv =m rv =m rr ( ) ( ) ( ) 29 Angular Momentum of a Particle • Moment of linear momentum of a particle – Mass times components of the velocity that are perpendicular to the moment arm h = rmv =m rv ( ) ( ) • Cross Product: Evaluation of a determinant with unit vectors (i, j, k) along axes, (x, y, z) and (v , v , x y v ) projections on to axes z i j k x y z rv = = yv"zv i+ zv" xv j + xv" yv k ( ) ( ) ( ) z y x z y x v v v x y z 30Cross Product in Column Notation Cross product identifies perpendicular components of r and v i j k x y z rv = = yv"zv i+ zv" xv j + xv" yv k ( ) ( ) ( ) z y x z y x v v v x y z " yv"zv ( ) v z y" 0z y Column x ( ' ' notation ( v rv= zv" xv = ( ) z 0x ' x z y ' ( ' ' (y x 0 v xv" yv z ( ) ' y x ( ' 31 Angular Momentum Vector is Perpendicular to Both Moment Arm and Velocity yv"zv ( ) z y ( ( h = mrv= m zv"xv ( ) x z ( ( xv"yv ( ) y x ( ' v 0"z y x ( ( v = m = mrv z 0"x ( y ( ( ( "y x 0 v z '( ' 32Specific Angular Momentum Vector of a Satellite … is the angular momentum per unit of the satellite’s mass, referenced to the center of attraction m h = rv=rv=rr S m Perpendicular to the orbital plane 33 Equations of Motion for a Particle in an InverseSquareLaw Field Acceleration is " µ r t µ ( ) I a t =v t =r t = = r t ( ) ( ) ( ) ( ) 2 3 ' r (t) r(t) r (t) … or µ r+ r= 0 3 r 34Cross Products of Radius and Radius Rate Then µ " rr+ r = 0 3 ' r rr= 0 rr= 0 … because they are parallel Chain Rule for Differentiation d rr = rr + rr = rr ( ) ( ) ( ) ( ) dt 35 Specific Angular Momentum µ µ " rr+ r = rr + rr ( ) ( ) 3 3 ' r r 0 d dh S = rr = =0 ( ) dt dt Consequently h = Constant S h h = r r" (Perpendicular to the plane of motion) ( ) S Orbital plane is fixed in inertial space 36Eccentricity Vector is a Constant of Integration µ µ r+ r'h=r'h+ r'h=0 3 3 r r " µ µ rh=" rh=" r rr ( ) 3 3 r r With triple vector product identity (see Supplement) µ µ d r rh=" r(rr)=" (rr"rr)= µ ( 3 2 ' r r dt r Integrating r (rh)dt =rh= µ +e ( " ' r e = Eccentricity vector Constant of integration ( ) 37 Significance of Eccentricity Vector T r r ),), rh"µ +e h=0 because rh"µ +e =0 (( +.+. '' r r T µr h T T (r"h) hµe h=0 r 0 0 T "µe h=0 e is perpendicular to angular momentum, which means it lies in the orbital plane Its angle provides a reference direction for the perigee 38General Polar Equation of a Conic Section ) r, T r rh"µ +e =0 ( +. ' r st 1 term is angular momentum squared T T T 2 r (rh)=h (rr)=h h =h Then T " r r 2 T 2 T h = µ r+r e = µ r+recos = 0 hµ +r e =0 ( ) ( ) ' r 2 h µ r= 1+ecos 39 Elliptical Planetary Orbits Assume satellite mass is negligible compared to Earth’s mass Then Center of mass of the 2 bodies is at Earth’s center of mass Center of mass is at one of ellipse’s focal points Other focal point is “vacant” 2 p h µ r= = , m or km 1+e cos 1+e cos : True Anomaly= Angle from perigee direction, deg or rad 40Properties of Elliptical Orbits Eccentricity can be determined from apogee and perigee radii rr rr a p a p e= = r +r 2a a p r = a /(1e) r = a /(1+e) p a Semimajor axis is the average of the two r +r a p a= 2 41 Properties of Elliptical Orbits Semilatus rectum, p, can be expressed as a function of h or a and e 2 2 p =h µ =a 1e ( ) Semiminor axis, b, can be expressed as a function of r a and r p b= rr a p Area of the ellipse, A, is 2 2 2 A= ab=a 1"e , m 42Energy is Inversely Proportional to the SemiMajor Axis At the periapsis, r p 2 2 r = 0 and v= r p =h µ =a 1e p p p ( ) 2 2 1 µ 1 h µ E" E = r" =" ( ) S p p 2 r = a 1e ( ) p 2 r 2 r r p p p 1 µ 2 " E = µp2µr = 1e2 1e ( ) ( ) ( ) p 2 2 2r 2a(1e) p µ(1e) = 2a 1e ( ) µ E = 43 2a Classification of Conic Section Orbits Orbit Shape Eccentricity, e Energy, E SemiMajor SemiLatus Axis, a Rectum, p Circle 0 0 0 a 2 Ellipse 0 e 1 0 0 a(1– e ) Parabola 1 0Undefined 2r p (") 2 Hyperbola 1 0 0 a(1– e ) 44“Vis Viva (Living Force) Integral” Velocity is a function of radius and specific energy 1 µ µ 2 v = + E v= 2 + E " 2 r r Specific total energy, E, is µ inversely proportional to the E = semimajor axis 2a Velocity is a function of 2 1 " v= µ radius and semimajor axis ' r a 45 Maximum and Minimum Velocities on an Ellipse Velocity at periapsis " " 2 1 2 1 v = µ = µ p ' ' r a a(1e) a p µ 1+e ( ) = a (1e) Velocity at apoapsis " 2 1 µ (1e) v = µ = a ' r a a 1+e ( ) a 46Velocities at Periapses and Infinity of Parabola and Hyperbola Parabola (a "" ") 2µ v = v =0 p r p Hyperbola (a 0) " 2 1 µ v = µ v =" p ' r a a p 47 Relating Time and Position in Elliptical Orbit Rearrange and integrate angular momentum over time and angle d d( Ellipse Area) 2 h = µp = r = 2 dt dt t"" f f f 2 3 r p d" dt = d" = 2 µ µp 1+ecos" ( ) t"" o o o … but difficult to integrate analytically 48Anomalies of an Elliptical Orbit Angles measured from last periapsis (or "): True Anomaly E (or ): Eccentric Anomaly M: Mean Anomaly 49 Kepler’s Equation for the Mean Anomaly M =EesinE From the diagram, ae+ pr e ( ) ar cosE = = a ae r =a 1cosE ( ) r = E aesinE 50Relationship of Time to Eccentric Anomaly 2 " µ 1 h 2 E = = r +µr 2 ' 2a 2 r 2 2 where h = µa 1e ( ) ( 2 1+ " 2 2 2 2 then r r = µ ra 1e ( ) ' r a ), 2 3 a dE 2" leading to (1ecosE) = 1 ' µ dt or µ dt =(1ecosE)dE 3 a 51 Integrating the Prior Result … t E 1 1 µ dt = 1"ecosE dE ( ) 3 a t E 0 0 E µ 1 (t"t )=(E"esinE) = M"M 3 1 0 1 0 E 0 a Time is proportional to Mean Anomaly MM MM 1 0 1 0 tt = or t =t + ( ) 1 0 1 0 3 3 µ a µ a Orbital Period, P µ 2" P0 = EesinE = 2" ( ) ( ) 3 0 a 3 P = 2" a µ 52Orbital Period Orbital period is related to the total energy 3 2 a"µ P= 2 = where E 0 for an ellipse 3 µ 2E Mean Motion, n, is the inverse of the Period 3 a 2 P= 2 where n is the Mean Motion µ n 53 Position and Velocity in Orbit at Time, t Mean Anomaly, given time from perigee passage, t p µ M(t)= tt ( ) p 3 a Eccentric Anomaly, E(t), given Mean Anomaly, M(t) E(t)esinE(t)= M(t) Newton’s method of successive approximation, using M(t) as starting guess for E(t) E (t)= M(t)+esinM(t) o Iterate until M Tolerance i M = M(t)"E (t)"esinE (t) i i i E =M 1"ecosE t ( ) i+1 i i E (t)= E (t)+E 54 i+1 i i+1Position and Velocity in Orbit at Time, t Calculate True Anomaly, given Eccentric Anomaly 1+e E(t) "1 (t)= 2tan tan ( 1"e 2 ' Compute magnitude of radius 2 a 1e ( ) r t = ( ) 1+ecos"(t) 55 Position and Velocity in Orbit at Time, t Radius vector, in the orbital plane x t r t cos' t ( ) ( ) ( ) r t = = ( ) y t r t sin' t ( ) ( ) ( ) "" Velocity vector, in the orbital plane v (t)'sin( t ( ) x µ v t = = ( ) p v t ( ) e+cos((t) y " " see Weisel, Spaceflight Dynamics, 1997, pp. 6466 56e ex xt t T Ti im me e:: ""lla an ne et ta ar ry y D De ef fe en ns se e 57 u up pp plle em me en nt ta all M Ma ar ri ia all 58First Point of Aries (Ecliptic Intercept at Right)" 59 Dimension of energy" Scalar (1 x 1)" Dimension of linear momentum" Vector (3 x 1)" Dimension of angular momentum" Vector (3 x 1)" 60SubOrbital (Sounding) Rockets 1945 Present NASA Wallops Island Control Center Canadian LTV Black Brant XII Scout 61 MATLAB Code for Flat Earth Trajectories Script for Script for Numerical Solution Analytic Solution tspan = 40; Time span, s xo = 10;100;0;0; Init. Cond. g = 9.8; t = 0:0.1:40; t1,x1 = ode45('FlatEarth',tspan,xo); vx0 = 10; Function for vz0 = 100; Numerical Solution x0 = 0; z0 = 0; function xdot = FlatEarth(t,x) x(1) = vx vx1 = vx0; x(2) = vz vz1 = vz0 – gt; x(3) = x x1 = x0 + vx0t; x(4) = z z1 = z0 + vz0t 0.5gt. t; g = 9.8; xdot(1) = 0; xdot(2) = g; xdot(3) = x(1); xdot(4) = x(2); xdot = xdot'; 62 endTrajectories Calculated with FlatEarth Model Constant gravity, g, is the only force in the model, i.e., no aerodynamic or thrust force Can neglect motions in the y direction Dynamic Equations Analytic (ClosedForm) Solution v t = v ( ) x x 0 v t =g z positive up ( ) ( ) z v (T)= v x x 0 x(t)= v (t) x T z(t)= v (t) v T = v gdt = vgT ( ) z z z z " 0 0 0 Initial Conditions x(T)= x +v T 0 x 0 v (0)= v x x 0 T 2 v (0)= v z z z T = z +v T gtdt = z +v TgT 2 ( ) 0 0 z 0 z " 0 0 0 x 0 = x ( ) 0 63 z(0)= z 0 Trajectories Calculated with FlatEarth Model v 0 =10m/s ( ) x v (0)=100,150,200m/s z x 0 = 0 ( ) z 0 = 0 ( ) 64MATLAB Code for SphericalEarth Trajectories Script for Numerical Solution R = 6378; Earth Surface Radius, km tspan = 6000; seconds options = odeset('MaxStep', 10) xo = 7.5;0;0;R; t1,x1 = ode15s('RoundEarth',tspan,xo,options); for i = 1:length(t1) v1(i) = sqrt(x1(i,1)x1(i,1) + x1(i,2)x1(i,2)); r1(i) = sqrt(x1(i,3)x1(i,3) + x1(i,4)x1(i,4)); end function xdot = RoundEarth(t,x) x(1) = vx x(2) = vy x(3) = x x(4) = y Function for mu = 3.98105; km2/s2 r = sqrt(x(3)2 + x(4)2); Numerical Solution xdot(1) = mu x(3) / r3; xdot(2) = mu x(4) / r3; xdot(3) = x(1); xdot(4) = x(2); xdot = xdot'; 65 end Equations that Describe Ellipses 2 2 x y + = 1 2 2 a b b a : Semi major axis, m or km y o x a b : Semi minor axis, m or km x = a cos ( ) ( ) y =b sin ( ) ( ) : Angle from x axis (origin at center) rad 66Constructing Ellipses FP +FP =FP +FP =2a 1 1 2 1 1 2 2 2 Foci (from center), x 2 2 2 2 f ' a'b a'b = , y 0 0 f "" " 1,2 String, Tacks, and Pen Archimedes Trammel Hypotrochoid 67 Ellipses Semi latus rectum ("The Parameter"), 2 b p= , m a 2 b Eccentricity,e= 1 2 a 2 b =1e; b=a 1e 2 a 68How Do We Know that Gravitational Force is Conservative Because the force is the derivative (with respect to r) of a scalar function of r called the potential, V(r): µ µ V(r)=m +V =m +V o o 1 2 T r r r ( ) " Vx " x ' V(r) µ ' y =Vy = m =(F ' g 3 ' r r ' ' z Vz This derivative is also called the gradient 69 of V with respect to r Conservation of Energy Energy is conserved in an elastic collision, i.e. no losses due to friction, air drag, etc. “Newton’s Cradle” illustrates interchange of potential and kinetic energy in a gravitational field 70Examples of Circular Orbit Periods for Earth and Moon Period, min Altitude above Surface, km Earth Moon 0 84.5 108.5 100 86.5 118 1000 105.1 214.6 10000 347.7 1905 71 Typical Satellite Orbits GPS Constellation Sun 26,600 km Synchronous Orbit 72GeoSynchronous Ground Track Geo Synchronous Ground Track 42, 164 km Marco Polo1 2 73 Background Math" Triple Vector Product Identity" a bc" a ic b a ib c ( ) ( ) ( ) T T = a c b a b c ( ) ( ) Dot Product of Radius and Rate " dr T r r = r ir = r dt 74
Website URL
Comment
sharer
Presentations
Free
Document Information
Category:
Presentations
User Name:
Dr.TomHunt
User Type:
Teacher
Country:
United States
Uploaded Date:
23-07-2017