Software Development
Low Thrust Orbit Propagation Software
There are two well-known methods one can use to propagate an orbit in the presence of a perturbing force, i.e., the small thrust produced by an electric engine. The first way uses the Lagrange Planetary Equations (LPE); the second way uses an integrator to calculate the future orbital positions based on the forces imbedded in the integrator’s conditions. We implemented the later method because we felt it would be more robust if we wished to add additional perturbing forces in the future, such as solar radiation pressure and third body forces from the planets. Furthermore, the Langrange Equations typically require conservative forces and the low thrust of an electric engine is not conservative. We will describe both processes here, as they are both very relevant to our study. But first we need to give a brief description of the Keplerian orbital elements.
Keplerian Orbital Elements
A spacecraft’s orbit requires six elements to describe its position and velocity at a point in time. The two main conventions in use are the standard Cartesian orbital set [x, y, z, vx, vy, vz] and the Keplerian orbital set [a, e, i, ω, Ω, ν], where
There are other Keplerian elements that can be used as well; common examples are E the Eccentric Anomaly and M, the Mean Anomaly, often replacing ν in the orbital set. A graphical representation of these elements is shown below.
In an orbit with no perturbing forces only the element ν changes with time. The other elements strictly describe the orbit’s shape and position relative to the central body. This is the advantage over the Cartesian orbital set where all six of the elements change with time. However, when perturbing forces are added to the system, more of the elements change with time. For example, adding a J2 effect from an oblate central body changes the three elements ω, Ω, and ν over time on a long time-scale and every element in the set over time on a short time-scale. This is where the Lagrange Planetary Equations come in.
Lagrange Planetary Equations
If we have a potential function described by R, we can determine how each of the orbital elements changes with time. R can be the gravitational potential of a central body, a third body force, or any other conservative force, F, modeled by a potential function. The temporal changes of the orbital elements are found by the following LPEs:
|
|
|
|
where h = na2(1-e2)1/2. If we have a force, either conservative or non-conservative, that can be expressed as:
such that r, t, and n are unit vectors that represent coordinates in the Radial, Along-Track, and Cross-Track coordinate system, then the temporal changes of the orbital elements can be written in the Gaussian form of the Lagrange Planetary Equations:
|
|
where,
|
|
In this case, one needs to know the accelerations of the orbiting mass in each direction (radial, along-track, and cross-track). In order to propagate the Keplerian orbital elements forward in time, one simply integrates these equations over the appropriate time interval. This is simple until the accelerations of the object become very complex, taking many perturbing forces into account. We have opted to use the simpler approach using a numerical method instead of an analytical method.
Numerical Propagation
We have developed and written our propagation software in Matlab. Matlab’s libraries include many integrators that would be appropriate to our study; we have chosen to use ODE45 as our propagation integrator. This routine takes in the initial conditions of a spacecraft and outputs the orbital characteristics of that spacecraft at given future times. The integrator also requires knowledge of how the various components change with time. For our purposes we had two conditions: the spacecraft always experienced the central gravitational force from the host body but the spacecraft did not always experience the perturbing force of the electric propulsion. In Cartesian coordinates, the acceleration due to the host body’s gravitational attration can be expressed as:
where, once again, μ is equal to the product of G – the gravitational constant – and M, the mass of the body supplying the gravitational force. This suffices for the relations for the change of velocity with respect to time of the object without the presence of any other force. To add the perturbing electric propulsion, one simply adds each component of the electric acceleration to the appropriate component of the central-body’s acceleration. For our considerations, we only considered the case such that the low-thrust vector T pointed in the exact opposite direction of the velocity vector:
,
where the bars represent the magnitude of the vector.
Matlab Code
This algorithm was designed to determine the orbit of a spacecraft as it transferred from a Low-Earth orbit to a Geosynchronous orbit using a low-thrust constant propulsion. The algorithm set the spacecraft’s initial position to be at an altitude of 300 km above the surface of the earth and its velocity to be about 7.73 km/s in the y-direction, typical of an object in LEO. We assumed the object’s initial LEO orbit and final GEO orbit would be coplanar and circular. The software then requests a thrust magnitude from the user – appropriate values are between 10 m/s and 0.01 m/s. The algorithm then propagates the spacecraft’s orbit forward in time until the moment that its apogee radius coincides with the radius of its target GEO orbit. The magnitude of the spacecraft’s thrust then goes to zero and the program continues to propagate the spacecraft’s cruise phase forward in time. At the moment when the spacecraft attains GEO radius, the program determines how much ΔV is required to circularize the orbit (using an impulsive maneuver). Finally, it outputs all of the relevant information it acquired from the procedure.
A few more comments about the propagation procedure… The program calls ODE45 numerous times during its procedure. At the beginning of a propagation phase the integrator only propagates the orbit forward in time for a fraction of the total duration required to reach GEO – on the order of 10-20 minutes of simulated time. As the spacecraft approaches its phase’s stopping condition, the algorithm instructs ODE45 to integrate at smaller and smaller increments, approaching a propagation period of 1 second. In this way, the error values of the numerical integration are reduced and the simulation does not require any excessive amount of ΔV to reach its goal.
Required subroutines: calc_ra.m, two_body_low.m
This algorithm runs exactly the same as LEO2GEO.m, except that it does not output anything to the screen except the plots, which it also saves to appropriate files. These plots are shown in LEO to GEO Transfers Section.
Required subroutines: calc_ra.m, two_body_low.m
<<< Back to LEO to GEO Using Low-Thrust Maneuvers
This algorithm is identical to LEO2GEO.m except that it works to produce interplanetary transfers instead of transfers between earth orbits. The user chooses which planetary orbit to transfer to (without interacting with any other planet) and how much thrust to use. In the case of Mercury and Venus, the user must input a negative thrust, corresponding to a thrust-vector aligned with the spacecraft’s velocity vector. The spacecraft’s initial conditions are equivalent to Earth’s position and velocity about the sun: 1 AU from the sun traveling at a velocity close to 30 km/s. As above, this program assumes every planetary orbit to be circular and coplanar. The program still alters the integration increments to reflect how close the spacecraft is to a stopping condition.
Required subroutines: calc_ra.m, calc_rp.m, two_body_low.m
This algorithm works the same as Low_Thrust.m except as a function to record all of the information about an interplanetary transfer orbit implementing many hardwired thrust-values.
Required subroutines: calc_ra.m, calc_rp.m, two_body_low.m
This algorithm works the same as Low_Thrust.m except it is hardwired to produce the 72 plots shown in the Interplanetary Transfer section.
Required subroutines: calc-ra.m, calc_rp.m, two_body_low.m
<<< Back to Low-Thrust Transfers
This algorithm combines the program Low_Thrust.m with gravity assists from other planets en-route to the destination. While it propagates the spacecraft’s orbit forward in time in each of the two phases, it continually checks to see if the spacecraft is approaching another planet’s orbit. If it finds a planetary orbit to be near, the program reduces the integrator’s propagation interval until it reaches the planet. Then, assuming the planet itself is in the intersecting part of the orbit (this is one area that could be improved in a project-extension), the program asks the user whether or not a gravity assist is desired. If it is not desired, the program increases the propagation interval once again and continues on its merry way; if an assist is desired, the program completes one. It does so by assuming the user wishes to obtain the greatest ΔV possible from the planet (this is another area that could be improved in an extension). Hence, the spacecraft’s new velocity will point in the same direction as the planet’s velocity. The ΔV is calculated and the new velocity is recorded; thereafter the program continues to propagate the orbit until a stopping condition, or another planet, is encountered.
Required subroutines: calc_ra.m, calc_rp.m, two_body_low.m, Check_Encounter.m, Set_increment_in.m, Set_increment_out.m
<<< Back to Gravity Assists
