Download MATLAB Assignment: Two-Body Problem and Torque-Free Attitude Dynamics and more Assignments Electrical and Electronics Engineering in PDF only on Docsity!
2020/21 EEE3039 - MATLAB Assignment
N. Baresi, H. Holt & N. Bernardini
- Introduction: The aim of this assignment is to facilitate and demonstrate understanding of the two-body problem and torque-free attitude dynamics. The coursework is divided into two parts. The first part of the assignment focuses on orbital mechanics and the equations of the two-body problem. You will be challenged with the task of integrating the equations of motion in the Earth-Centered Inertial reference frame (ECI) and comparing your solution with the analytical propagations enabled by Kepler’s equation. You will also be asked to derive the equations of motion in the Earth-Centered Earth-Fixed reference frame (ECEF) and to generate ground tracks using a provided MATLAB function that outputs the Earth sub-satellite point of a spacecraft. Finally, you will be tasked with designing a Sun-synchronous orbit and the orbital maneuver required to reach the desired orbit configuration. The second part of the coursework deals with attitude kinematics and Euler’s equations. You will be asked to represent the orientation of a satellite with respect to different reference frames and to propagate both Euler’s equations and the kinematics equations for a desired time interval. This numerical experiment will provide you with the opportunity to investigate spacecraft orientations and witness the advantages and disadvantages of the different attitude representations that will be introduced throughout the module.
- Assignment Task: Write a report describing your results, derivations, and observations. The report should be no longer than 20 pages including Figures and without counting the Appendix. When writing your report, please make an effort to explain your reasoning and the equations being used. Your MATLAB code should be included in the Appendix and as part of your submission on SurreyLearn (either in a “.zip” file or as “.m” files). Do not include screenshots, as these do not comply with the professional standards of industry and research reports.
- Due Date: Please submit a PDF of your final report online on SurreyLearn by Tuesday, November 24, 16:00 (BST). A “Coursework Assignment” folder is available for your submissions on SurreyLearn under “Assessment”.
Coursework Data
- μ = 398600.4418 [km^3 /s^2 ];
- R⊕ = 6378.137 [km];
- J 2 = 1082. 63 × 10 −^6 ;
- Ixx = 2500 [kg m^2 ];
- Iyy = 5000 [kg m^2 ];
- Izz = 6500 [kg m^2 ];
- Spacecraft inertia matrix in principal axes frame: I =
Ixx 0 0 0 Iyy 0 0 0 Izz
1 The Two-Body Problem
A satellite is being launched into a highly-inclined Earth orbit with initial orbit elements:
a = 7178. 137 km e = 0. 07 i = 96 deg Ω = 30 deg ω = 65 deg M 0 = 182 deg
1.1 Solve Kepler’s Equation using Newton’s Method and Find Initial True Anomaly
Your lectures have introduced Kepler’s equation as a way to link the mean anomaly of a satellite (M ) with its eccentric anomaly (E) and eccentricity (e):
M = E − e sin E. (1)
Later on, we learned how to relate the mean anomaly of a satellite to the time along its orbit using the time equation M = M 0 + n (t − t 0 ), (2)
where n =
GM/a^3 is the mean motion of the satellite. Using Newton’s method, Eq. (1) can be solved numerically to find the value of E that satisfies Eq. (1) for any suitable pair of e ∈ [0, 1) and M values. Once the value of E is found, the true anomaly of the spacecraft can be calculated via
tan
θ 2
1 + e 1 − e
tan
E
Newton’s method is an iterative method for solving nonlinear equations of the form f (x) = 0 via successive approximations. Starting from a reliable initial guess, x(0), one can expand f (x(0)^ + δx(0)) using a Taylor’s series approximation and find
f (x(0)^ + δx(0)) = f (x(0)) + f ′(x(0)) δx(0)^ + H.O.T. (4)
Neglecting high-order terms and assuming f (x(0)^ + δx) = 0 yields
δx(0)^ = −
f (x(0)) f ′(x(0))
Using Eq. (5), it is possible to update the value of the initial guess as x(1)^ = x(0)^ + δx(0)^ and repeat the process until f (x(k)) < T ol, i.e, f (x(k)) is smaller than a given threshold (we recommend using T ol = 10−^10 )
- Create a MATLAB function [E]=Kepler(e, M, tol) to solve Kepler’s equation using Newton’s method. Use E(0)^ = M to initialize the iterative algorithm and output the value of the Eccentric anomaly in radians; [6 marks]
- Report the value of the initial Eccentric Aanomaly of the satellite found via Newton’s method as well as the number of iterations needed to achieve convergence. Then, apply Eq. (3) to find the initial value of θ and comment on whether the satellite is closer to its apoapsis or periapsis point. [4 marks]
1.5 Ground Tracks
Starting from the trajectory of the spacecraft in ECEF over 10 orbital periods,
- Calculate the longitude and latitude of the spacecraft in the ECEF frame and plot their values as a function of time (Tip: consider using the yyaxis command in MATLAB to save space in your report); [ marks]
- Pass these values to the provided function PlotGroundTracks.m (type help PlotGroundTracks for info on inputs and outputs) and generate a 2D plot of the satellite ground-track. Enclose your figure in the report, comment your findings, and discuss whether the satellite roughly passes over Guildford, UK. [ marks]
1.6 Orbital Perturbations and Maneuvers
Given the semi-major axis and eccentricity of a satellite as in Eq. (1),
- Calculate the required inclination to place the satellite in a Sun-synchronous orbit and report your value. Apart from J2, what would be the other major perturbation experienced by a satellite at these altitudes? [4 marks]
- Determine the velocity of the spacecraft at the ascending and descending node and report your values; [ marks]
- Explain where you would execute an inclination change maneuver and provide reasons for your conclusion; At your chosen point, calculate the ∆V required to change the inclination of the spacecraft from i = 96 deg to the Sun-synchronous orbit value found in 1.6.1; [5 marks]
- Calculate the ∆V expenditure required for the inclination change if you were to implement the maneuver in the more expensive point. What is the difference with respect to the ∆V cost observed in 1.6.3? [ marks]
2 Attitude Representations and Torque-Free Motion
One hour before the inclination change maneuver, the spacecraft has ECI position and velocity coordinates
r =
(^) [km], v =
(^) [km s−^1 ] (7)
2.1 Direction Cosine Matrices
Given the state of the satellite as in Eq. (7),
- Find the direction cosine matrix
[
I O
]
from the orbital frame of the satellite O to the ECI Frame I and report its values in the report; [3 marks]
- Knowing that the current direction cosine matrix from the orbital frame to the principal axes frame of the satellite can be represented as a 3-2-1 sequence of single axis rotations via [ B O
]
= R 1 (γ) R 2 (β) R 3 (α),
where α = 7 deg, β = 4 deg, and γ = 2 deg, compute
[
B I
]
, the direction cosine matrix from the ECI frame to the principal axes frame of the spacecraft, and report its values in your report. [3 marks]
2.2 Attitude Representations
Starting from the values of
[
B I
]
found in 2.1.2,
- Calculate and report the values of Euler’s Principal Axis and Angle; [3 marks]
- Calculate and report the values of Euler’s Parameters (Quaternions); [3 marks]
- Calculate and report the values of a 3-2-1 Euler Angles sequence (yaw-pitch-roll). [3 marks]
2.3 Kinematics and Euler’s Equations
Let θ 1 , θ 2 , and θ 3 be the yaw, pitch, and roll angles found in 2.2.3, respectively. Let us also denote
ωB/I =
ω 1 ω 2 ω 3
= [B O]
n
as the angular velocity of the spacecraft as seen from the principal axes frame B one hour before the in- clination change maneuver. Given that n =
μ/a^3 is the mean motion of the satellite and defining X = [ θ 1 , θ 2 , θ 3 , ω 1 , ω 2 , ω 3
]T
as the attitude state of the spacecraft
- Write a MATLAB function [dXdt] = AttitudeDynamics(t, X, I) that calculates the time derivative of X assuming that no external torques are acting on the satellite (torque-free motion). Include the equations for X˙ in your report while providing explicit formulas for your kinematics equation (i.e., do not simply state that θ˙ = B(θ) ω, but provide explicit expressions for B(θ)); [7 marks]
- Use ode45 or ode113 to integrate the attitude dynamics of the spacecraft for 3600 s and output its state every 10 s. Plot and report the time histories of Euler angles and angular velocity components. Comments on whether this attitude representation approaches singularities; [4 marks]
- Calculate the magnitude of the angular momentum and rotational kinetic energy of the spacecraft at every step of the numerical integration. Plot these values as a function of time and comment on the conservation of these two quantities; [4 marks]
- Create a Polhode plot of the angular momentum sphere and energy ellipsoid using MATLAB’s sphere and ellipsoid functions. Evaluate your plot and describe whether the satellite is close to a maximum or minimum energy configuration. [4 marks]
