Small Satellites

Home » Space Flight/Orbital Mechanics Excercise

Category Archives: Space Flight/Orbital Mechanics Excercise

Advertisements

Bi-Elliptic Hohmann Transfer

In this example we calculate the total change in speed required for a bi-elliptic Hohmann transfer from a geocentric circular orbit of 7200 km radius to circular orbit of 125000 km radius. The apogee of the first transfer ellipse is 190000 km.

 clc;
 clear all;

 R_i  = 7200;     % [km]
 R1_a = 190000;   % [km]
 R_f  = 125000;   % [km]

mu   = 398600;   % [km^3/s^2] Earth’s gravitational parameter
% For initial circular orbit
V_i = (mu/R_i)^0.5;
% Speed at apogee and perigee for the first transfer ellipse
V1_a = (2*mu*R_i/(R1_a*(R1_a+R_i)))^0.5;
V1_p = (2*mu*R1_a/(R_i*(R1_a+R_i)))^0.5;
% Semimajor axes of the first transfer ellipse
a1 = (R_i + R1_a)/2;
% Speed at apogee and perigee for the second transfer ellipse
V2_a = (2*mu*R_f/(R1_a*(R1_a+R_f)))^0.5;
V2_p = (2*mu*R1_a/(R_f*(R1_a + R_f)))^0.5;
% Semimajor axes of the second transfer ellipse
a2 = (R_f + R1_a)/2;
% For target circular orbit
V_f = (mu/R_f)^0.5;

% For bi-elliptic maneuver the total speed change required
dV =  abs(V_i - V1_p)+ abs(V1_a - V2_a) + abs(V_f - V2_p); %[km/s]

% Time required for transfer
t_bi = pi/(mu)^0.5*(a1^1.5+a2^1.5);     %[s]

fprintf('Total speed change = %4.4f [km/s]\n',dV);
fprintf('Time required for transfer = %4.2f [hours]\n',t_bi/3600);
Total speed change = 3.9626 [km/s]
Time required for transfer = 129.19 [hours]

 

Advertisements

Simple orbital mechanics problems

An unmanned satelite orbits the Earth

clear all;
clc;
rp = 7000;                      %[km] Perigee radius
ra = 70000;                     %[km] Apogee radius
mu = 398600;                    %[km^3/s^2]

ecc = (ra-rp)/(ra+rp)           %Eccentricity of the orbit
a   = (ra+rp)/2                 %[km] Semimajor axis
T   = 2*pi/mu^0.5*a^1.5/3600    %[h] Period of the Orbit
E_spec = -mu/(2*a)              %[km^2/s^2]Specific energy
ecc =

    0.8182

a =

       38500

T =

   20.8833

E_spec =

   -5.1766

Find true anomaly when Alt = 1000 km;

R_earth = 6378;
theta = acos(a*(1 -ecc^2)/((R_earth+1000)*ecc)-1/ecc)*180/pi %[deg]
theta =

   27.6069

Velocity components at this anomaly

hmu = (a*(1 -ecc^2)/mu)^0.5;
vr = 1/hmu *ecc*sin(theta)          %radial component
vt = 1/hmu*(1 + ecc*cos(theta))     %tangential component
vr =

    2.8342

vt =

    2.0001

Velocities at Apogee and Perigee

vp = 1/hmu*(1 + ecc)
va = 1/hmu*(1 - ecc)
vp =

   10.1751

va =

    1.0175

The spacecraft is in a 250 km by 300km low eart orbit. How long does it take to coast from perigee to appogee. Semimajor axis of this eliptical orbit

a = R_earth + (250+300)/2; % [km]
% Time which takes to fly from Perigee to Appogee equals to half of the
% Orbit period
Th = pi/mu^0.5*a^1.5/60 %[min]
Th =

   45.0045

The altitude of a satelite in an elliptical orbit around the earth is 1600 km at apogee and 600 km at perigee.

 ra = R_earth + 1600;
 rp = R_earth +600;
% Eccentricity of the Orbit
ecc = (ra - rp)/(ra + rp)
% Semimajor axis;
a = (ra + rp)/2
% The orbital Speed at perigee and apogee
hmu = (a*(1 -ecc^2)/mu)^0.5;
vp = 1/hmu*(1 + ecc)
va = 1/hmu*(1 - ecc)
% The period of the Orbit
T = 2*pi*a^1.5/mu^0.5/60 %[min]
ecc =

    0.0669

a =

        7478

vp =

    7.8065

va =

    6.8280

T =

  107.2601

A satellite is palced into an orbit at perigee at an altitude of 1270 km with a speed of 9 km/s. Calculate the flight path angle gamma and altitude of the satelite at a true anomaly of 100deg.

altp  = 1270;    %[km]
vp    = 9;       %[km/s]
theta = 100;     %[deg]
rp = altp + R_earth;
h = rp * vp
ecc = h^2/(rp*mu) - 1
gamma = atan( ecc*sind(theta)/(1 + ecc*cosd(theta)));
gamma = 180/pi*gamma    %[deg]
h =

       68832

ecc =

    0.5542

gamma =

   31.1256

Altitude of the satelite at a true anomaly of 100 [deg].

rd = h^2/mu*(1/(1 + ecc* cosd(theta)));
altd = rd - R_earth
altd =

  6.7738e+003

A satelite is launched into the orbirt at an altitude of 640 km with a speed 9.2 km/s and flight path angle of 10 deg. Calculate the true anomaly of the lunch poinnt and the period of the orbit

alt = 640;           %[km]
rs  = R_earth + alt; %[km]
vel = 9.2;           %[km/s]
gamma = 10;          %[deg]
vp = vel*cosd(gamma);
vr = vel*sind(gamma);
h  = rs * vp;
ac = (h^2/(mu*rs)-1);
bc = (vr*h/mu);
theta = atan(bc/ac)*180/pi  %Anomaly of the lunch point
%Period of the orbit
%Eccentricity
ecc = ac/cosd(theta);
T = 2*pi/mu^2*(h/(1 - ecc^2)^0.5)^3  %[sec]
theta =

   29.7830

T =

  1.6075e+004

A satellite has a perigee and apogee altitudes of 250 km and 42000km. Calculate the orbit period, eccentricity and the maximum speed.

rp = R_earth + 250; %km
ra = R_earth + 42000;%km

%Semimajor axis of this eliptical orbit
a =( rp + ra )/2;
T = 2*pi/mu^0.5 * a ^(3/2)/60       %[min]
ecc = (ra - rp)/(ra+rp)
h = (mu*(1 + ecc)*rp)^0.5;
vp = h/rp                           %[km/s]
T =

  756.5368

ecc =

    0.7590

vp =

   10.2852

A rocket lunched from the surface of the earth has a speed 8.85 km/s when powered flight ends at an altitude of 550 km. The flight path angel at this time is 6 degree.Determinine eccentricity

vel  = 8.85 ;           %[km/s]
ra = R_earth + 550;     %[km]
gamma = 6;              %[deg]
vd = vel * cosd(gamma);
h  = ra * vd;
ect = h^2/(mu*ra)-1;
vr = vel* sind(gamma);
est = h*vr/mu;
theta = atan(est/ect);
ecc = ect/cos(theta)
ecc =

    0.3742

Published with MATLAB® 7.10

Euler rotation example, Rotation matrix, Quaternion, Euler Axis and Principal Angle

A classical Euler rotation involves first a rotation about e3 axis, then one about the e1 axis and finally a rotation about the e3 axis.

% Compute the elements of rotation matrix R
clc;
psi   = pi/4;
theta = pi/4;
fi    = pi/6;

R3_psi   = [cos(psi) sin(psi) 0;
            -sin(psi) cos(psi) 0;
                0  0        1]
R1_theta = [1 0 0;
            0 cos(theta)  sin(theta);
            0 -sin(theta) cos(theta)]
R3_fi    = [cos(fi)   sin(fi) 0;
             -sin(fi) cos(fi) 0;
            0 0 1]
R3_psi =

    0.7071    0.7071         0
   -0.7071    0.7071         0
         0         0    1.0000

R1_theta =

    1.0000         0         0
         0    0.7071    0.7071
         0   -0.7071    0.7071

R3_fi =

    0.8660    0.5000         0
   -0.5000    0.8660         0
         0         0    1.0000

Multiplying each rotation matrix

 R = R3_fi*R1_theta*R3_psi
R =

    0.3624    0.8624    0.3536
   -0.7866    0.0795    0.6124
    0.5000   -0.5000    0.7071

From Rotation Matrix Compute Quaternion

q4 = 0.5*(1 + R(1,1)+ R(2,2) + R(3,3))^0.5;
q1 = (R(2,3) - R(3,2))/(4*q4);
q2 = (R(3,1) - R(1,3))/(4*q4);
q3 = (R(1,2) - R(2,1))/(4*q4);
q  = [q1 q2 q3 q4]
norm_q = norm(q)            % Checking that the norm of q = 1
qv     = [q1 q2 q3];
q =

    0.3794    0.0500    0.5624    0.7330

norm_q =

    1.0000

Compute Euler axis and its components along the axis E1,E2,E3 unit vector. Compute the Euler principle angle.

Euler_angle = 2*acos(q4)*180/pi            % [deg]
Euler_axis  = qv/sind(Euler_angle/2)       % [E1,E2,E3]
norm_E      = norm(Euler_axis)
Euler_angle =

   85.7293

Euler_axis =

    0.5577    0.0734    0.8268

norm_E =

    1.0000

 

Radial Release

Example 10.8, Orbital Mechanics for Engineering Students, 2nd Edition.

A satellite is to be completely despin using a two-mass yo-yo device with tangential release. Assume the spin axis of moment of inertia of the satellite is C = 200 kg · m^2 and the initial spin rate w0 = 5 [rad/s]. The total yo-yo mass is 4 kg, and the radius of the spacecraft is 1 meter.

clc;
C  = 200;       % [kg*m^2]
w0 = 5;         % [rad/s]
m  = 4;         % [kg]
R  = 1;         % [m]
K = 1 + C/(m*R^2); %     Nondimensional factor K

The required cord length l

l = R*(K)^0.5       %[m]
% Time t to despin
t = (K)^0.5/w0      %[s]
% Tension in the yo-yo cables
t = 0:0.001:5;
N = 2*K*C*w0^3*t./(R*(K + w0^2*t.^2).^2);
plot(t,N);
xlabel('time [s]');
ylabel('Tension N [N]');
grid on;
[maxN,i] = max(N);
l =

    7.1414

t =

    1.4283

The maximum tension in the yo-yo cables [N];

Tension in the yo-yo cables
maxN
% which ocurs at
t = i*0.001  % [s]
maxN =

  454.7542

t =

    0.8260

The speed of the masses at release;

v = R*w0*(K)^0.5  % [m/s]
v =

   35.7071

The angle rotated by the satellite during the despin;

theta = (pi/2 - 1)*(K)^0.5   % [rad]
theta =

    4.0763

Cord length required for complete despin

l = R*(K^0.5 - 1)            % [m]
l =

    6.1414

Published with MATLAB® 7.10

 

Inertial components of the angular acceleration

Example 10.4, Orbital Mechanics for Engineering Students, 2nd Edition.

The inertial components of the angular momentum of a torque-free rigid body are

Hg = [320; -375; 450 ];       %[kg*m^2/s] IJK
% the Euler angles [deg] are
fi      =  20;
theta   =  50;
psi     =  75;
% The inertia tensor in the body-fixed principal frame is
Ig = [1000, 0,    0;
         0, 2000, 0;
         0, 0,    3000]; %[kg*m^2]

Obtain the inertial components of the (absolute) angular acceleration Matrix of the transformation from body-fixed frame to inertial frame

QxX = [-sind(fi)*cosd(theta)*sind(psi) + cosd(fi)*cosd(psi), ...
 -sind(fi)*cosd(theta)*cosd(psi) - cosd(fi)*sind(psi),sind(fi)*sind(theta);
 cosd(fi)*cosd(theta)*sind(psi) + sind(fi)*cosd(psi),...
  cosd(fi)*cosd(theta)*cosd(psi) - sind(fi)*sind(psi),-cosd(fi)*sind(theta);
    sind(theta)*sind(psi),   sind(theta)*cosd(psi),    cosd(theta)
       ]
QxX =

    0.0309   -0.9646    0.2620
    0.6720   -0.1740   -0.7198
    0.7399    0.1983    0.6428

Matrix of the transformation from inertial frame to body-fixed frame

QXx = QxX'
QXx =

    0.0309    0.6720    0.7399
   -0.9646   -0.1740    0.1983
    0.2620   -0.7198    0.6428

Obtain the components of HG in the body frame

Hgx = QXx*Hg            % [kg*m^2/s]
Hgx =

   90.8616
 -154.1810
  643.0376

The components of angular velocity in the body frame

Ig_inv = inv(Ig);
wx = Ig_inv*Hgx         % [rad/s]
wx =

    0.0909
   -0.0771
    0.2143

From Euler ’s equations of motion we calculate angular acceleration in the body frame.

alfa_x = - Ig_inv*cross(wx,Ig*wx)               % [rad/s^2]
alfa_x =

    0.0165
    0.0195
    0.0023

Angular acceleration in the inertial frame

alfa_X = QxX*alfa_x                             % [rad/s^2] IJK
alfa_X =

   -0.0177
    0.0060
    0.0176

Published with MATLAB® 7.10

 

%d bloggers like this: