EN2231 – Laboratório de Guiagem, Navegação e Controle Atividade 2: Elementos orbitais
André Mitsuo Ravazzi 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 ( ) 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:
= 12 − =− 2 Onde,
éé çã / é= í ≅3,986.10 /
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 − −(.) .
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 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− = 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.
cos− si n − 0 1=[−sin0 − cos− ] 0 0 1 1 0 0 2=[00 cos− ] si n − sin − cos− cos− si n − 0 3=[−si0n − cos− ] 0 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, temse 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 proposta da atividade, como sendo da ISS.
foram pré-definidos na
=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− = 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