Download MATLAB Lab: Numerically Confirming Kepler's Laws and Simulating Planetary Motion and more Exercises Classical Mechanics in PDF only on Docsity!
16.07 Lab I
Issued: October 16, 2009 Due: Monday, October 26, 2009
Introduction
This laboratory is the first in a series that deal with celestial and spacecraft dynamics. In this lab, you will create a MATLAB code that will allow you to numerically confirm Kepler’s laws as well as simulate the motion of the earth and the other planets in the solar system as well as spacecraft in orbit about or on the way to one of them. Since this lab will focus on the earth and moon system, the presence of the sun will be ignored. In order to accomplish these tasks, your MATLAB code will integrate the equations of motion in time, starting from an initial condition. In later labs, you will add features like spacecraft thrusting, which may allow you to simulate orbital transfer, rendezvous, descent and ascent from planetary surfaces and other more complicated situations.
Figure 1: 3 Body coordinate system.
Figure 1 shows the coordinate system for a generic problem involving the three planetary/spacecraft bodies. The bodies are labeled as Body 1, Body 2 and Body 3, but these names are just an example. These bodies will turn into the Earth, the Moon and a spacecraft. The important thing is to choose a convention that you understand and can be consistent with. The first step is to formulate the equations we would like to solve with the MATLAB code. In the general case, we desire to calculate the positions and velocities of each body as a result of its interaction with the other members of the system. For our purposes, we only consider the gravitational forces resulting from the mutual attraction of body pairs. Thus, the motion of any of the three moving bodies must consider the gravitational influence of the other two moving bodies. For instance, the attractive force that Body 2 exerts on Body 1 is given by,
F 12 =
μ 2 m 1 3 r^12.^ (1) r 12
Here, μ 2 = Gm 2 , the gravitational parameter of Body 2, where G is the universal gravitational constant^1 and m 2 is the mass of Body 2. The mass of Body 1 is denoted by m 1 , and r 12 is the magnitude of the distance vector that points from Body 1 to Body 2, r 12 = r 2 − r 1.
Celestial bodies move in three dimensions, so in general you would write your equations for the x, y,and z positions of each body. However, for this lab, we will take the position of the earth, moon and spacecraft system and their motions as confined to the x, y plane. Therefore, the z equations are not required. For numerical purposes it is convenient to write the equation of motion r¨ = F /m for each body as a system of first order equations. That is, we write,
r˙ = v, (2a)
v˙ = F /m. (2b)
Therefore, to describe the motion of a body in two dimensions, four state variables are required – two for position and two for velocity. For three moving bodies confined to planar motion in the x, y plane,we need only four for each body for a total of twelve. We will use a Cartesian coordinate system which means that the equations should be separated into Cartesian x,and y components. The state vector X contains all variables that completely describe the state of the system. For purposes of uniformity, you should order the states as follows,
X = [x 1 y 1 x 2 y 2 x 3 y 3 u 1 v 1 u 2 v 2 u 3 v 3 ]T^ , (3)
where x,and y are position components, u is the x-component of velocity, and v is the y-component of velocity. You will be solving several different problems in this lab and building on your code in future labs. Thus, it will be beneficial to code in a flexible way. For example, you might want to define variables for the masses of the objects, instead of hard-coding the values. Some problems in this or future labs may only deal with two bodies such a planet in orbit around the sun, or in others, the gravitational pull of one of the bodies may be neglected^2. However, your code must have the capability to handle up to three bodies. An effective way to remove the effect of a body on the other bodies in your code is to set its gravitational parameter to zero.
For later work, we present physical parameters for 5 celestial bodies: the sun, the earth, Mars, Jupiter and the moon. Because of the extremely large mass of the sun, it is usual to consider it to be stationary in our coordinate system. The planets, including the earth, orbit around the sun; and the moon orbits around the earth. When considering planetary orbits around the sun or the moon’s orbit around the earth, the effect of other celestial bodies are ignored. These simplifications are referred to as the two-body and the three-body problems. The two body problem has an analytic solution in term of the Kepler planetary orbit problem. The three-body solution has no analytic solution and as we shall later see can exhibit quite wild behavior.
(^1) G = 6. 673 × 10 − (^11) m (^3) /kg-s 2 (^2) For example, the pull of the space shuttle on the Earth could certainly be neglected.
solving equations numerically, it is standard practice to change variables such that the unit length and unit time are characteristic of your system. We will use the bar to denote non-dimensionalized variables. Thus, ¯r and t¯^ will denote non-dimensionalized distance and time, respectively. Note that in this gravitational problem we do not need to scale the mass since we are working with gravita tional parameters, μ. In the general F = ma problem, we would need to scale length, time, and mass. We will use the radius of the earth as our length scale, and 1 day (24 hours) as our time scale.
In your report, you need only write the equations of motion for the three-body earth- moon plus spacecraft system, and their derivation, in the non-dimensionalized bar variables. You should obtain μearth = 11472 and μmoon = 141. 3 in your non-dimensionalized equations.
2. Observing Kepler’s Laws
Once you have the equations of motion, you should write your solver using MATLAB. There is a function in MATLAB called ode45 which nicely solves differential equations of the form given in Equation 4. For details on the use of this function, type help ode45 into MATLAB. Another command associated with ode45 is a command to set the error tolerances. The ode45 function uses a 4-stage Runge-Kutta integration technique, and has a feature to vary the time step to control the error. For example, if the curvature of an orbit gets very high, ode45 will reduce the time step in order to carefully resolve the trajectory. These error tolerances are set with the command:
options=odeset(‘RelTol’,1e-6,‘AbsTol’,1e-6,‘MaxStep’,.001);
The values shown here are a good starting point. If your calculation takes too long to run you can loosen the tolerances, and if you need more detail in the solution, you can tighten them.
Your solver should consist of a MATLAB script file that sets the solver options, sets the ini tial state, and then calls ode45. One of the input parameters to ode45 is a function that calculates the state derivative. In addition to the script file, you will need to write a function to define your system of equations.
Once you have a working solver, test it by simulating the orbit of a geostationary satellite in a planar x, y orbit around a fixed Earth. In simulating the motion of a satellite, we ignore the effect of the satellite on the earth, and take the earth to be stationary at the center of our coordinate system. A geostationary satellite orbits at a radius of approximately 42,042 km from the center of the earth so that it remains fixed above a point on the rotating earth; it rotates in inertial space in a period of 23 hrs, 56 minutes and 4.09 sec ( because a the earth rotates about the sun and in a time that is close but not exactly equal to 1 year.) You need to transform these initial values into your non-dimensional framework setting x(0), y(0), u(0) and v(0) equal to the non-dimensionalized form of the geostationary satellite parameters. Use the solution to the satellite orbit problem as your initial conditions. Simply initialize the satellite to start at (R 0 , 0) with velocity (0, v 0 ), where R 0 and v 0 have been non-dimensionalized. The orbit should be nearly circular (try for good accuracy so that the later parts of the lab make sense; strive for say .001) and if the simulation is run for one day, the planet should complete just about one orbit. This is a good quick check that your solver is working.
Kepler’s First Law
The orbits of the planets are ellipses, with the Sun at one focus of the ellipse. (We will use the earth and a satellite to validate Kepler’s Laws; Kepler didn’t have satellites.)
Once your solver is working, you will simulate the eccentric (elliptical) orbit of a satellite around a fixed Earth. Let’s call it object A. Start the object at (10, 0) with an initial velocity of (0, 20). From your simulation, determine the eccentricity and period of object A’s orbit, and describe how you determined these values from your numerical results. Compare these results to the theoretical values you would expect, given the initial conditions used. The lecture notes from lectures L15 and L16 will give you the details of the theory needed. Additionally, the Appendix at the end of this document gives some key formulas and labels critical parts of an elliptical orbit.
Kepler’s Second Law
The line joining a satellite and the fixed Earth sweeps out equal areas in equal times as the object follows its orbit. (Again, we will use a satellite in orbit about the earth.)
Now, using the same simulation of object A, calculate the area swept by the object as it completes one orbit. You can do this as a simple post-processing operation once you have determined the position of the object as a function of time. The function ode45 will return a vector t with time values, and an array X with the value of the state variables at each time. From this you should be able to compute a new vector Area that contains the value of the area swept by the object for each time. Plot the integrated area swept as a function of time, and verify that the area increases linearly.
Kepler’s Third Law
The ratio of the squares of the orbital periods of two objects is equal to the ratio of the cubes of their semimajor axes.
Here we consider a different planet, call it object B. The initial velocity for this object is the same as for object A, but its starting position is twice as distant from the earth. Verify the ratio in Kepler’s Third Law, for object A and B using the results of your numerical simulation. Turn in plots showing both orbits, and label the semimajor axes for both. Your plots should show the results of the independent simulations of objects A and B; you should not simulate them simultaneously.
The ratio of a orbiting object’s period squared to its semimajor axis cubed can be shown to equal 4 π^2 /μ. Verify this ratio from your numerical simulations. Remember to use non-dimensional values.
4. Earth-Moon Simulation-Rotating Coordinates
Now program the equations for the earth moon relative motion in a rotating coordinate system. You have already derived the expression for the inertial acceleration of a point in rotating coordi nates; the gravitation forces between the bodies are added without modification. Turn in your derivation of these governing equations. Now, if you equations are correct, you should be able to place the earth and the moon at their initial positions and see no motion: i.e they will be fixed in your coordinate system. You may have to adjust the values of Ω slightly to get them to stay absolutely fixed in the rotating system; shouldn’t be a big adjustment.
Report Preparation
There is no need to produce a fancy report. You should include all the figures requested and answer all the questions posed. If your write-up is neat, organized, and readable, then there is no need to type it.
Turn in the value of Ω that leads to stationary positions of the earth and moon in the rotating coordinate system. Submit both the script file and the function file as a
zip file to the course website.
Appendix
Figure 2 labels key parts of the elliptical orbit, where π is the periapsis, or point of closest approach, α is the apoapsis, or the furthest point in the orbit, O is the location of the gravitational center, the planet being orbited, a is the semimajor axis, b is the semiminor axis, p is called the parameter of the conic, and e is the eccentricity of the orbit. The following are a few relationships you might need:
b = a 1 − e^2 , (5)
rπ = 1 +
p e
= a(1 − e), (6)
rα =
p = a(1 + e). (7) 1 − e
To find the relationship between the velocity of the orbiting planet and its distance from the gravitational center, the energy equation can be used,
1 2 μ μ 2
v − r
2 a
Figure 2: An elliptical orbit with key values labeled.
