EN2231 – Laboratório de Guiagem, Navegação e Controle

Atividade 2: Elementos orbitais

André Mitsuo Ravazzi RA 11002511

Gyslla Danielle Bento da Silva RA 11049811

Vagner Pascualinotto Junior RA 11050907

Relatório apresentado à

Universidade Federal do ABC

como parte dos requisitos para

aprovação na disciplina EN2231

– Laboratório de Guiagem,

Navegação e Controle

Prof., Dr. Leandro Baroni

SANTO ANDRÉ – 2014

Introdução

Neste relatório apresentaremos uma outra forma de obtenção da órbita de

um V/E partindo-se dos vetores posição 𝑹 e velocidade 𝑽. A órbita será obtida pela definição de seus elementos orbitais. Estes elementos servem para caracterizar a órbita, definem o seu tamanho, a sua forma e a sua orientação, além de definir a localização do V/E.

Orbita elíptica: De acordo com a primeira lei de Kepler, toda órbita é uma cônica (circulo,

elipse, parábola, hipérbole) onde a Terra ocupa um dos focos. Nas duas últimas, o satélite só passa uma vez perto da Terra. Nos interessam as órbitas circular e elíptica. O círculo é apenas um caso particular de elipse com excentricidade e = 0, ou seja, os focos ocupam o mesmo lugar no centro do círculo. Segue abaixo, uma visualização de uma órbita elíptica:

Figura 1: Órbita elíptica.

Onde: a = semi-eixo maior; b = semi-eixo menor; c = meio distância entre focos F e F´; p = semi-latus rectum; r = raio do satélite ao centro da Terra; R = raio de Terra;

𝜃 = anomalia verdadeira; rp = raio do perigeu; ra = raio do apogeu; hp = altura do V/E no perigeu ha = altura do V/E no apogeu;

Elementos Orbitais Clássicos: Os EOCs que definem a órbita de um satélite são:

Semi-eixo maior (𝑎)

Excentricidade (𝑒) Inclinação (𝑖) Ascensão reta do nodo ascendente ou ARNA (𝛺) Argumento do perigeu (𝜔0)

Anomalia verdadeira (𝜃)

Figura 2: Definição da ARNA.

Figura 3: Definição da inclinação, argumento e anomalia verdadeira da órbita.

Órbita Geoestacionária: Para que o satélite permaneça numa posição fixa em relação a um ponto na Terra, é preciso que esteja em uma órbita geoestacionária GSO, com as seguintes características: - órbita deve ser geossíncrona, ou seja, ter o mesmo período de revolução da Terra, que é de um dia sideral = 23 h 56 min 04s. Portanto, deverá ter semi-eixo maior de 42164 km; - órbita deve ter excentricidade zero, ou seja, deve ser circular. O semi-eixo maior é o raio r da órbita. A altura do satélite será, portanto, de 35786 km; - órbita deve ter inclinação zero, ou seja, o satélite deve girar no plano do equador e no mesmo sentido da rotação da Terra

Objetivos

Revisão dos EOC (Elementos Orbitais Clássicos) que determinam e definem uma órbita. Determinação dos elementos, logo, das órbitas, a partir de medições de R e V. Apresentar gráficos da órbita determinada, utilizando o Matlab.

Metodologia A caracterização da órbita (determinação do tamanho, da forma, da

orientação, e a localização do Veículo Espacial), será dada pelos Elementos Orbitais Clássicos (EOC), que serão detalhados abaixo. Conforme proposta da atividade, a determinação dos elementos orbitais será dada a partir do vetor

Posição 𝑹 e do vetor Velocidade 𝑽.

Semi-eixo maior, 𝒂: É a meia distância ao longo do eixo principal de uma elipse. Especifica o

tamanho da órbita e se relaciona com a sua energia. É obtido a partir da relação da energia mecânica especifica dado pelas equações abaixo:

𝜀 =1

2𝑉2 −

𝜇

𝑅

𝑎 = −𝜇

2 𝜀

Onde,

𝑉 é 𝑜 𝑣𝑒𝑡𝑜𝑟 𝑣𝑒𝑙𝑜𝑐𝑖𝑑𝑎𝑑𝑒 𝑒𝑚 (𝑘𝑚/𝑠)

𝑅 é 𝑜 𝑣𝑒𝑡𝑜𝑟 𝑝𝑜𝑠𝑖çã𝑜 𝑒𝑚 (𝑘𝑚)

𝜀 é 𝑎 𝑒𝑛𝑒𝑟𝑔𝑖𝑎 𝑒𝑠𝑝𝑒𝑐í𝑓𝑖𝑐𝑎

𝜇 = 𝐺 𝑀𝑇 ≅ 3,986.105𝑘𝑚3/𝑠2

Excentricidade, 𝒆: Especifica a forma da órbita e diz que tipo de seção cônica ela

representará. O vetor excentricidade é calculado da seguinte forma:

𝑒 =1

𝜇[(𝑉2 −

𝜇

𝑅) �⃗⃗� − (�⃗⃗�. �⃗⃗�)�⃗⃗�]

O vetor 𝑒 é sem unidade, tendo direção apontada para o perigeu e que tem magnitude igual a excentricidade da órbita, 𝑒.

Inclinação, 𝒊: Especifica a orientação/deslocamento angular do plano orbital com

relação ao plano fundamental, o plano do equador. Lembrando-se que 𝑖, é o

ângulo entre ℎ ⃗⃗⃗ ⃗(𝑣𝑒𝑡𝑜𝑟 𝑚𝑜𝑚𝑒𝑛𝑡𝑜 𝑎𝑛𝑔𝑢𝑙𝑎𝑟 𝑒𝑠𝑝𝑒𝑐í𝑓𝑖𝑐𝑜, 𝑒𝑚 𝑘𝑚2/𝑠) e

�⃗⃗�(𝑣𝑒𝑡𝑜𝑟 𝑢𝑛𝑖𝑡á𝑟𝑖𝑜 𝑛𝑎 𝑑𝑖𝑟𝑒çã𝑜 𝑑𝑜 𝑝𝑜𝑙𝑜 𝑛𝑜𝑟𝑡𝑒 𝑐𝑒𝑙𝑒𝑠𝑡𝑒), tem-se :

ℎ⃗⃗ = �⃗⃗� × �⃗⃗�

𝑖 = 𝑎𝑟𝑐𝑐𝑜𝑠 (�⃗⃗� . ℎ⃗⃗

𝑘 ℎ)

Ascensão reta do nodo ascendente, 𝜴: Especifica a orientação do plano orbital com relação à direção principal do

SGI (𝐼 − 𝑝𝑜𝑛𝑡𝑜 𝛾). Definindo o vetor �⃗⃗� na direção do nodo ascendente (linha dos nodos), temos:

�⃗⃗� = �⃗⃗� × ℎ⃗⃗ Assim,

𝛺 = 𝑎𝑟𝑐𝑐𝑜𝑠 (𝐼 . �⃗⃗�

𝐼 𝑛)

Nesse tipo de elemento orbital clássico, deve-se realizar uma verificação

com relação à posição do ângulo de ascensão reta do Nodo Ascendente, 𝛺 :

𝑆𝑒 𝑛𝑗 ≥ 0 ; 𝑒𝑛𝑡ã𝑜 0 ≤ 𝛺 ≤ 180°

𝑆𝑒 𝑛𝑗 < 0 ; 𝑒𝑛𝑡ã𝑜 180° ≤ 𝛺 ≤ 360°

Argumento do perigeu, 𝝎: Especifica a orientação da órbita dentro do plano orbital, assim como

resposta ele dá a localização do perigeu, ou seja, o ponto mais próximo da órbita de um astro. O argumento do perigeu é dado pelo ângulo entre o nodo

ascendente (direção dada pelo 𝑣𝑒𝑡𝑜𝑟 �⃗⃗� ) e o perigeu (direção dada pelo 𝑣𝑒𝑡𝑜𝑟 𝑒), conforme segue abaixo:

𝜔 = 𝑎𝑟𝑐𝑐𝑜𝑠 (�⃗⃗⃗� . �⃗⃗⃗�𝑛 𝑒

)

Nesse tipo de elemento orbital clássico, deve-se realizar uma verificação

com relação à posição do ângulo do argumento do perigeu, 𝜔 :

𝑆𝑒 𝑒𝑘 ≥ 0 ; 𝑒𝑛𝑡ã𝑜 0 ≤ 𝜔 ≤ 180° 𝑆𝑒 𝑒𝑘 < 0 ; 𝑒𝑛𝑡ã𝑜 180° ≤ 𝜔 ≤ 360°

Anomalia verdadeira, 𝝂: Especifica a localização do veiculo espacial dentro do plano Orbital. Assim,

é formado pelo ângulo entre o perigeu (𝑣𝑒𝑡𝑜𝑟 𝑒) e o (𝑣𝑒𝑡𝑜𝑟 �⃗⃗�), conforme segue abaixo:

𝜈 = 𝑎𝑟𝑐𝑐𝑜𝑠 (�⃗⃗⃗� . �⃗⃗⃗�𝑒 𝑅

)

Nesse tipo de elemento orbital clássico, deve-se realizar uma verificação

com relação à posição do ângulo da anomalia verdadeira, 𝜈:

𝑆𝑒 (�⃗⃗�. �⃗⃗�) ≥ 0 ; 𝑒𝑛𝑡ã𝑜 0 ≤ 𝜈 ≤ 180°

𝑆𝑒 (�⃗⃗�. �⃗⃗�) < 0 ; 𝑒𝑛𝑡ã𝑜 180° ≤ 𝜈 ≤ 360°

Em seguida, utilizando-se a equação da órbita abaixo e os elementos orbitais calculados, será possível obter a relação (𝑅, 𝜐) , e assim, plotar graficamente a órbita em coordenadas polares (2D).

𝑹 = 𝑎(1 − 𝑒2)

1 + 𝑒 cos (𝝂)

Figura 4: Sistema de coordenadas polares (2D).

Para obter a representação gráfica em coordenadas cartesianas, será preciso submeter a relação (𝑅, 𝜐) a uma série de operações. Mas para isso, faremos uma consideração: - Sabe-se que os EOC são calculados com base no SGI (R e V vetores no SGI). Assim, a relação entre ambos é direta e a passagem da descrição da órbita de um sistema para outro pode ser obtida facilmente. Para simplificação deste problema, será proposto que o Sistema de Coordenadas Cartesianas Terrestres (x,y,z, com x em Greenwich) permanecerá coincidente com o SGI (X,Y,Z, com X no ponto de áries). A posição do V/E no plano xy da órbita, onde: x na direção do perigeu e y a 90 graus (anti-horário); e z na direção do momento angular.

𝑃𝑝𝑜𝑙𝑎𝑟 = [𝑅 ∙ cos 𝜐 𝑅 ∙ sin 𝜐 0]

Neste caso, três rotações levam o vetor 𝑃𝑝𝑜𝑙𝑎𝑟 , descrito no sistema do

plano da órbita, para coordenadas cartesianas inerciais (SGI),

𝑃𝑥𝑦𝑧 = [𝑟𝑥 𝑟𝑦 𝑟𝑧].

Obs.: todas são no sentido horário, ou seja, em ângulos negativos.

1) 𝝎 em torno do eixo z; 2) 𝒊 em torno do eixo x;

3) 𝜴 em torno de z.

Figura 5: Visualização das rotações.

As matrizes relativas a estas rotações, respectivamente, são dadas a seguir.

𝑟𝑜𝑡𝑧1 = [cos −𝜔 sin −𝜔 0

− sin −𝜔 cos −𝜔 00 0 1

]

𝑟𝑜𝑡𝑥2 = [1 0 00 cos −𝑖 sin −𝑖0 sin −𝑖 cos −𝑖

]

𝑟𝑜𝑡𝑧3 = [cos −𝛺 sin −𝛺 0

− sin −𝛺 cos −𝛺 00 0 1

]

Com uso destas três rotações, na ordem citada, obtém-se a posição SGI do V/E, a partir da posição no plano da órbita:

𝑃𝑥𝑦𝑧 = 𝑟𝑜𝑡𝑧3 ∙ 𝑟𝑜𝑡𝑧2 ∙ 𝑟𝑜𝑡𝑧1 ∙ 𝑃𝑝𝑜𝑙𝑎𝑟

Assim, para cada vetor posição em coordenadas no plano da órbita, tem-se o mesmo vetor descrito em coordenadas cartesianas geocêntricas inerciais. Como consideramos aqui que estas são coincidentes com o SCT (sistema cartesiano terrestre), temos condição de obter, a partir destas, as coordenadas de latitude e longitude do V/E para cada posição

𝑃𝑥𝑦𝑧 = [𝑟𝑥 𝑟𝑦 𝑟𝑧]

que ele ocupe.

Resultados

O MATLAB foi a ferramenta utilizada para a obtenção dos resultados. O algoritmo pode ser dividido em três partes. A primeira possui a função de calcular os elementos orbitais. A segunda possui a função de gerar as projeções da órbita em coordenadas polares (2D, no plano da órbita) e em coordenadas cartesianas (3D, no Sistema Geocêntrico Inercial). A terceira possui a função de projetar a trajetória da órbita calculada, na primeira parte, sobre a superfície terrestre. Esta projeção será dada pelo método de Mercator.

Parte A

Os valores dos vetores posição 𝑹 e velocidade 𝑽 foram pré-definidos na proposta da atividade, como sendo da ISS.

𝑹 = 4890.7 𝑖 − 5224.8 𝑗 − 850.1�⃗⃗�

𝑽 = −1.4 𝑖 + 0.1𝑗 − 7.3 �⃗⃗�

Assim, obtivemos os elementos keplerianos da órbita da ISS:

Figura 6: elementos keplerianos da órbita da ISS

Parte B

A partir da relação abaixo e dos valores de semi-eixo maior, excentricidade e anomalia verdadeira obtidos no MATLAB, pudemos plotar a

relação do raio 𝑹 versus o ângulo 𝝂. E assim obtivemos o traçado da órbita em coordenadas polares sobre o plano da órbita.

𝑅 = 𝑎(1 − 𝑒2)

1 − 𝑒 cos (𝝂)

Figura 7: Órbita da ISS em coordenadas polares.

Parte C A seguir, através de operações matriciais o traçado em coordenadas

polares foi transformado em coordenadas cartesianas inerciais. Obtive-se então:

Figura 8: Órbita da ISS em coordenadas cartesianas no Sistema Geocêntrico Inercial.

As mesmas técnicas e scripts foram aplicados para se obter os elementos keplerianos e órbitas dos satélites StarOne C2 e Molniya 1-91, que seguem abaixo.

StarOne C2

Figura 9: elementos keplerianos da órbita do satélite StarOne C2.

Figura 10: Órbita do satélite StarOne C2 em coordenadas polares.

Figura 11: Órbita do satélite StarOne C2 em coordenadas cartesianas no Sistema Geocêntrico Inercial.

Molniya 1-91

Figura 12: elementos keplerianos da órbita do satélite Molniya 1-91.

Figura 13: Órbita do satélite Molniya 1-91 em coordenadas polares.

Figura 15: Órbita do satélite StarOne C2 em coordenadas cartesianas no

Sistema Geocêntrico Inercial.

Conclusão

Neste relatório, apresentamos uma forma de obtenção da órbita de um

veículo espacial partindo-se dos vetores posição 𝑹 e velocidade 𝑽. A órbita foi obtida pela definição de seus elementos orbitais. Estes elementos tem a função de caracterizar uma órbita, definindo o seu tamanho, a sua forma, sua orientação, e a localização do veículo espacial.

Com uma programação relativamente simples, obtivemos os valores dos elementos orbitais requeridos para a obtenção da órbita em coordenadas cartesianas, no sistema SGI.

Anexo 1: código-fonte do script utilizado clear all

%Declaração de variáveis r = zeros(360,3); n = zeros(360,3); rx = zeros(360,3); ry = zeros(360,3); rz = zeros(360,3); la = zeros(360,3); lo = zeros(360,3);

%Dados de Entrada % ISS % vetorR = [4890.7 -5224.8 -850.1]; %vetor Posição % vetorV = [-1.4 -0.1 -7.3]; %vetor Velocidade

% StarOne C2 % vetorR = [3010.33 -42067.38 -0.59]; %vetor Posição % vetorV = [3.07 0.22 0.001]; %vetor Velocidade

% Molniya 1-91 vetorR = [10016.34 -17012.52 7899.28]; %vetor Posição vetorV = [2.50 -1.05 3.88]; %vetor Velocidade

%constante mi: produto entre a constante da gravitação universal (G) e

a massa da Terra (M) mi = 3.986*10^5; %Semi-eixo Maior

E = 0.5*(norm(vetorV)^2) - (mi/norm(vetorR)); %equação da energia

a = -mi/(2*E); %semi-eixo maior

%Excentricidade

vetor_e = (1/mi)*(((norm(vetorV))^2 - (mi/(norm(vetorR))))*vetorR -

(dot(vetorR,vetorV)*vetorV));

e = norm(vetor_e); %excentricidade

fprintf('e = %i\n\n\n\n', e);

%Inclinação da Órbita

vetorK = [0 0 1]; % vetor unitário na direção do pólo norte celeste

vetorH = cross(vetorR,vetorV); %vetor momento angular específico

i = acos(dot(vetorK,vetorH)/(norm(vetorK)*norm(vetorH))); %inclinação

fprintf('i = %i\n\n\n\n', i);

%Ascensão Reta do Nodo Ascendente

vetorI = [1 0 0]; %vetor unitário na direção principal

vetorN = cross(vetorK,vetorH); %vetor N: produto vetorial entre K e H

arna = acos(dot(vetorI,vetorN)/(norm(vetorI)*norm(vetorN))); %ascensão

reta do nodo ascendente

%verificação da posição do ângulo de ascensão reta do nodo ascendente if vetorN(2)<0 arna = 2*pi - arna; %arna = 360 - arna end

fprintf('arna = %i\n\n\n\n', arna);

%Argumento do Perigeu

w = acos(dot(vetorN,vetor_e)/(norm(vetorN)*norm(vetor_e))) ; %angulo

do perigeu %w = (acos(dot(vetorN,vetor_e)/(norm(vetorN)*norm(vetor_e))))*180/pi;

%verificação da posição do ângulo do perigeu if vetor_e(3)<0 w = 2*pi - w; %w = 360 - w end

fprintf('w = %i\n\n\n\n', w);

%Anomalia Verdadeira

ni = acos(dot(vetor_e,vetorR)/(norm(vetor_e)*norm(vetorR))); %angulo

da anomalia verdadeira ni = ni*180/pi;

%verificação da posição do ângulo da anomalia verdadeira if dot(vetorR,vetorV) < 0 ni = 360 - ni; else ni = ni*180/pi; end

fprintf('ni = %i\n\n\n\n', ni);

%Cálculo para plotagem dos gráficos

for m=1:360; R = (a*(1-e^2))/(1+e*cosd(m)); %calculo de R em função de ni r(m,1) = R; %vetor com valores de R n(m,1)=m*pi/180; %vetor com valores de ni

Ppolar = [R*cosd(m);R*sind(m);0]; %representação do vetor posição

no plano da orbita em coordenadas cartesianas

%operações de mudança da base para SGI rotz1 = [cos(-w) sin(-w) 0; -sin(-w) cos(-w) 0; 0 0 1]; %rotação w

sobre o eixo Z rotx2 = [1 0 0;0 cos(-i) sin(-i);0 -sin(-i) cos(-i)]; %rotação i

sobre o eixo X rotz3 = [cos(-arna) sin(-arna) 0;-sin(-arna) cos(-arna) 0; 0 0

1]; %rotação arna sobre o eixo Z

Pxyz = rotz1*rotx2*rotz3*Ppolar; %mudança de base para SGI

rx(m,1) = Pxyz(1,1); ry(m,1) = Pxyz(2,1); rz(m,1) = Pxyz(3,1);

end

% Plot da Órbita Elíptica em Coordenadas Polares (2D)

figure polar(n(:,1),r(:,1)); grid on; %plotando trajetória em coordenadas

polares title('Traçado da Órbita em Coordenadas Polares') xlabel('ni') ylabel('R')

% Plot da Órbita Elíptica em Coordenadas Cartesianas (3D)

figure plot3(rx,ry,rz),axis square;grid on;hold on; %plotando trajetória no

SGI title('Traçado da Órbita no SGI')

%Função da elipsoide que simula o globo terrestre

Req = 6378.137; % raio do equador Rp = 6356.7523; % raio polar

xc = 0; yc = 0; zc = 0; % índice c representa centro de massa, xc, yc

e zc são as coordenadas do centro de massa xr= Req; yr = Req; zr = Rp; % GRS ( Geodetic Reference System) o eixo

z é no sentido polar n = 24; % número de divisões feita na linha do equador e na linha

polar ou número de células ou densidade de malha % No caso 24, simbolizando 360º/24 = 15°. Uma divisão das

horas locais [x,y,z] = ellipsoid(xc,yc,zc,xr,yr,zr,n); % Comando que irá calcular

as coordenadas em 3D(x,y,z) o elipsóide em 3D(x,y,z) surfl(x,y,z) % Comando que desenha o gráfico em 3D(x,y,z)

xlabel('eixo x','FontSize', 16) ylabel('eixo y','FontSize', 16) zlabel('eixo z','FontSize', 16) colormap copper axis equal