Session 5: Review of Classical Astrodynamics

Session 5: Review of Classical Astrodynamics In previous lectures we described in detail the process to find the optimal specific impulse for a particular situation. Among the mission requirements that serve as constraints, is the change in velocity, Δv. The purpose of the propulsion system is to provide this Δv, which depends on the particular trajectory followed by the spacecraft and is determined by the orbital mechanics at play. Because of this, the use of concepts in astrodynamics becomes very relevant to space propulsion. In this Lecture, we present a short review of the fundamental laws in astrodynamics (Further reading: R. Battin, “An Introduction to the Mathematics and Methods of Astrodynamics”, AIAA Education Series, 1999). We start by analyzing the two-body problem with central forces. Such analysis provides the fundamentals of most orbital mechanics problems. Consider the situation depicted in the figure below: z m1  r1    r = r2 − r1  r2 m2 y x We have two point masses interacting through central (gravitational) forces. Under no external forces, the equations of motion for these particles are, d2 r 1 r2 − r1 = Gm1 m2 2 dt |r2 − r1 |3 d2 r 2 r2 − r1 m2 2 = −Gm1 m2 dt |r2 − r1 |3 m1 (1) Defining the relative position vector r = r2 − r1 , we subtract Eqs. (1) and obtain a single equation of motion, d2 r r +µ 3 =0 2 r dt or dv r +µ 3 =0 dt r (2) where µ = G(m1 + m2 ) and v = dr/dt. In most cases of interest in astrodynamics, one of the masses dominates over the other. For the particular situation in the figure, we could assume the motion of m1 is negligible for m1 » m2 , therefore µ = Gm1 . This is a perfectly acceptable approximation for orbital motion about planetary bodies. In the case of earth, µE = 398,601 km3 /s2 . This innocently looking equation (2) kept occupied several of the most prodigious minds in mathematics, like Gauss and Euler, ever since Newton proclaimed his gravitational law and in parallel with Leibniz, developed infinitesimal calculus. In general, there are two ways 1 to tackle this equation: i ) through numerical integration, by taking advantage of today’s computational power (a luxury, that neither Gauss of Euler ever hoped to tap into) and ii ) analytically, by uncovering a number of integrals of motion related to the structure of this non-linear differential equation. The analytical method, as expected, is much more enlightening and will be explored in some detail in this lecture. We begin by taking the cross product of vector ~r with Eq. (2),  ~r × d~v ~r +µ 3 dt r  = ~r × d~v d(~r × ~v ) d~r d(~r × ~v ) − × ~v = =0 = dt dt dt dt (3) where ~r × ~r = ~v × ~v = 0. From here, we see that the resulting vector ~r × ~v must be constant. We define this constant as, ~h = ~r × ~v (4) which represents a massless angular momentum, which is conserved when no external forces are applied to the masses. The angular momentum is our first integral of motion. We now take Eq. (2) and calculate its cross product with ~h,  d  ~r × ~h ~v × ~h + µ 3 = 0 dt r (5) We use the vector identity ~r × ~h = ~r × (~r × ~v ) = ~r(~r · ~v ) − r2~v and notice that ~r · ~v = rvr , where vr = dr/dt is the radial component of the velocity vector. The second term in Eq. (5) can be written as (excluding µ),   ~r × ~h vr ~v d ~r = 2 ~r − = − r3 r r dt r (6) Using this result, we can write Eq. (5) as,    d  d ~r ~ ~v × h − µ =0 dt dt r (7) This equation can be directly integrated, ~r ~v × ~h − µ = µ~e r (8) where the constant (vector) of integration was selected as µ~e for convenience, since later on will be identified as the eccentricity of a geometrical orbit (times µ). The eccentricity is our second integral of motion. We calculate the magnitude of the eccentricity vector, ~e · ~e = ~v × ~h ~r − µ r 2 !2 (9) Expanding the square on the right hand side, ~v × ~h ~r − µ r  2 ~v × ~h !2 =   ~v × ~h · ~r +1−2 µ2 (10) µr Here we note that (~v × ~h)2 = v 2 h2 since ~v ⊥ ~h. For the last term in Eq. (10), we use the vector identity ~v × ~h = ~v × (~r × ~v ) = ~rv 2 − rvr~v , and therefore, (~v × ~h) · ~r = r2 (v 2 − vr2 ) = r2 vθ2 = h2 The magnitude of the eccentricity vector becomes, e2 = v 2 h2 2h2 + 1 − µ2 µr → ET ≡
1 2 µ µ2 2 v − = e −1 2 2 r 2h (11) which is known as the vis-viva integral and represents the total energy ET (per unit mass), our third and last integral of motion. In this way, without explicitly solving the differential equation Eq. (2), it is possible to write expressions for the two-body problem in terms of integrals of motion. To see this, use Eq. (8) to calculate the dot product ~e · ~r, ~e · ~r = (~v × ~h) · ~r h2 −r = −r µ µ (12) Define f as the angle between the eccentricity and radial vectors (known as the true anomaly) and solve for r, r = h2 /µ 1 + e cos f or r = p 1 + e cos f (13) Eq. (13), known as the equation of orbit, is the solution in space (not in time) of the differential Eq. (2). It describes the relative planar motion of the two point masses considered in the problem. A closer look at Eqs. (11) and (13) reveals three motion regimes that depend on the magnitude of the eccentricity vector. 0 ≤ e < 1 ET < 0 bounded e=1 ET = 0 unbounded (zero velocity at ∞) e>1 ET > 0 unbounded elliptic parabolic hyperbolic Unbounded trajectories are interesting, especially in space propulsion where in some instances the objective is to achieve an escape trajectory for exploration missions. However, for the moment, we focus on bounded (elliptical) orbits. Eventually, we will look at the result of perturbing these orbits through small forces provided by electric propulsion. The planar motion in the bounded case (0 ≤ e < 1) defines an elliptical trajectory (circular, if e = 0), as shown in the figure below. From Eq. (13), we see that when f = π/2 the radius 3 ! r f ! e p 2a becomes h2 /µ. This quantity is typically represented by p and is known as the parameter of the orbit. The periaxis and apoaxis can be found directly from Eq. (13) with f = 0 and f = π, respectively, rp = p 1+e and ra = p 1−e (14) Since (twice) the semi-major axis of the orbit is 2a = ra + rp , we find, p = a(1 − e2 ) (15) Consequently, the orbital energy in Eq. (11) can be written as, ET = − µ 2a (16) Given the eccentricity and semi-major axis, we use Eq. (13) to find the orbital radius as a function of the true anomaly. Then, we use Eq. (11) to find the orbital velocity at that particular location. What remains is orbital timing. In other words, given an initial time and position, find where the orbiting object will be some time afterwards. Use the conservation of angular momentum, h = r2 (df /dt), and substitute in the equation of orbit, df = (1 + e cos f )2 � µ dt p3 (17) Appropriately, this equation is known as Kepler’s equation. The period of the orbit could be easily found by integrating Kepler’s law (orbit sweeps equal√areas in equal u times) since the area of the ellipse is πab and the semi-minor axis is b = a 1 − e2 = h a/µ, the familiar result is, � T = 2π a3 µ (18) We mentioned earlier that a number of brilliant mathematicians worked intensively in solving the two-body problem. In fact, most of that work was devoted to solve Kepler’s equation (17). Once the orbital timing is found, the solution of the two-body problem is complete. There are a number of different ways in which this equation could be written, each one with its own procedure to find an approximate solution. In general, Eq. (17) could be solved 4 numerically, but of course, little is gained that way since, after all, the original equation of motion Eq. (2) could also be solved numerically. At the end, approximate solutions using iterative methods, are extremely accurate. It is of course, more interesting for our particular application to solve the equations of motion Eq. (1) when there are forces acting on the masses (from thrusters, for instance). In particular, if we assume a force F~ is applied to m2 , the resulting equation of motion will be, d~v ~r + µ 3 = ~ac (~r, t) dt r (19) where ~ac = F~ /m2 is the acceleration due to external forces (other than gravity) acting on the moving mass. This new equation brings no additional difficulty if the problem is solved numerically. However, our analytical approach is crippled since it is no longer possible to find integrals of motion (neither the angular momentum, eccentricity or energy are longer constant). It is said that the resulting trajectories are non-Keplerian. Despite this, if the acceleration is a small quantity (which is typical for the type of maneuvers using electric propulsion), the constants of motion will drift slowly. This allows us to follow these changes and provide accurate predictions of the actual trajectory. Not surprisingly, this lies in the realm of the method of variation of parameters. For now, let us finalize this review by writing the vectorial equation of motion, Eq. (19), as two differential equations in polar coordinates. This will be useful later on to find analytical solutions to some simple problems where small accelerations are imparted and motion is confined to a constant orbital plane: d2 r −r dt2  dθ dt 2 + µ = ar r2 and 5 d2 θ 2 dr dθ aθ + = dt2 r dt dt r (20)