UNIVERSIDADE ESTADUAL PAULISTA "JÚLIO DE MESQUITA FILHO" CAMPUS DE SÃO JOÃO DA BOA VISTA ARNALDO PIERONI STECCA FILHO Aplicação de estratégias de controle low-thrust em veículos espaciais ao redor de Vênus São João da Boa Vista 2022 Arnaldo Pieroni Stecca Filho Aplicação de estratégias de controle low-thrust em veículos espaciais ao redor de Vênus Trabalho de Graduação apresentado ao Conselho de Curso de Graduação em Engenharia Aeronáutica do Campus de São João da Boa Vista, Universidade Estatual Paulista, como parte dos requisitos para obtenção do diploma de Graduação em Engenharia Aeronáutica . Orientadora: Profª Drª. Rita de Cássia Domingos Coorientador: Profº Dr. Antonio Fernando Berta- chini de Almeida Prado São João da Boa Vista 2022 S811a Stecca Filho, Arnaldo Pieroni Aplicação de estratégias de controle low-thrust em veículos espaciais ao redor de Vênus / Arnaldo Pieroni Stecca Filho. -- São João da Boa Vista, 2022 59 p. : il., tabs. Trabalho de conclusão de curso (Bacharelado - Engenharia Aeronáutica) - Universidade Estadual Paulista (Unesp), Faculdade de Engenharia, São João da Boa Vista Orientadora: Rita de Cássia Domingos Coorientador: Antonio Fernando Bertachini de Almeida Prado 1. Engenharia aeroespacial. 2. Órbitas. 3. Veículos espaciais. 4. Perturbação (Astronomia). I. Título. Sistema de geração automática de fichas catalográficas da Unesp. Biblioteca da Faculdade de Engenharia, São João da Boa Vista. Dados fornecidos pelo autor(a). Essa ficha não pode ser modificada. UNIVERSIDADE ESTADUAL PAULISTA “JÚLIO DE MESQUITA FILHO” FACULDADE DE ENGENHARIA - CÂMPUS DE SÃO JOÃO DA BOA VISTA GRADUAÇÃO EM ENGENHARIA AERONÁUTICA TRABALHO DE CONCLUSÃO DE CURSO APLICAÇÃO DE ESTRATÉGIAS DE CONTROLE LOW-THRUST EM VEÍCULOS ESPACIAIS AO REDOR DE VÊNUS Aluno: Arnaldo Pieroni Stecca Filho Orientador: Prof.ª Dr.ª Rita de Cássia Domingos Banca Examinadora: - Rita de Cássia Domingos (Orientadora) - Murilo Sartorato (Examinador) - Priscilla Andressa de Sousa Silva (Examinadora) A ata da defesa com as respectivas assinaturas dos membros encontra-se no prontuário do aluno (Expediente nº 041/2021) São João da Boa Vista, 04 de julho de 2022 Dedico este trabalho a meus pais. AGRADECIMENTOS Primeiramente, agradeço aos meus pais por todo amor, carinho, apoio e compreensão durante esta etapa da minha vida. Agradeço à professora Rita pelo apoio, paciência, convivência e pelos conselhos, fundamentais para o meu crescimento acadêmico e pessoal. Agradeço ainda ao professor Bertachini pela orientação essencial para a conclusão deste trabalho. De forma geral, agradeço aos professores e colaboradores da UNESP de São João da Boa Vista por todo o comprometimento e dedicação. Agradeço aos meus amigos que sempre estiveram presentes direta ou indiretamente em todos os momentos de minha formação. Ademais, ficam os meus agradecimentos ao apoio financeiro ao processo nº2016/024561-0, Fun- dação de Amparo à Pesquisa do Estado de São Paulo (FAPESP); processos nº123686/2019-9 e nº310317/2016-9, Conselho Nacional de Desenvolvimento Científico e Tecnológico (CNPq), funda- mentais para a realização deste trabalho. Por fim, o presente trabalho foi realizado com apoio da Coordenação de Aperfeiçoamento de Pessoal de Nível Superior – Brasil (CAPES) – Código de Financiamento 001. Este trabalho contou com o apoio das seguintes entidades: FAPESP - Fundação de Amparo à Pesquisa do Estado de São Paulo (Processo: 2016/024561-0) CAPES - Coordenação de Aperfeiçoamento de Pessoal de Nível Superior (Código 001) CNPq - Conselho Nacional de Desenvolvimento Científico e Tecnológico (Processos: 123686/2019-9 e 310317/2016-9) GDOP - Grupo de Dinâmica Orbital e Planetologia da Unesp de Guaratinguetá. “I may not have gone where I intended to go, but I think I have ended up where I needed to be.“ (Douglas Adams) RESUMO Dentre os diversos corpos do Sistema Solar, Vênus se destaca por ser o planeta com maior temperatura média superficial, além de possuir em sua atmosfera uma densa camada de nuvens que dificultam o seu estudo. Sendo assim, é fundamental a adoção de missões que permitam a observação contínua do planeta de modo a criar uma série histórica de levantamento de dados. Para isso, um possível modelo orbital a ser adotado é o de órbitas congeladas, uma vez que, essas órbitas apresentam os elementos orbitais excentricidade e argumento do periapsis praticamente constantes ao longo de um determinado período de tempo, ou seja, a altitude será sempre a mesma em um determinado ponto geográfico do planeta. Entretanto, veículos espaciais com missões definidas para este tipo de órbita sofrem perturbações que podem alterar a configuração orbital. Desse modo, faz-se necessária a adoção de correções orbitais através de manobras orbitais a fim de manter as características da missão proposta. No modelo considerado, o veículo espacial sofre perturbações devido ao potencial gravitacional do Sol e da não esfericidade de Vênus. Para fazer as correções orbitais na órbita do veículo, serão adotadas as estratégias de controle low-thrust contínuo visando encontrar as regiões em que ocorre o menor gasto energético para a missão. Por fim, para essa análise foram selecionadas condições orbitais relevantes ao problema como as condições críticas, regiões próximas à superfície e próxima à região em que as forças perturbativas são equivalentes. PALAVRAS-CHAVE: Controle low-thrust. Veículos espaciais. Manobras orbitais. Planeta Venus. ABSTRACT Among the bodies in the Solar System, Venus stands out for being the planet with the highest average surface temperature, in addition to having in its atmosphere a dense layer of clouds that make it difficult to study. Therefore, it is essential to adopt missions that allow continuous observation of the planet in order to create a historical series of data collection. For this, a possible orbital model to be adopted is frozen orbits model, since these orbits present the orbital elements eccentricity and periapsis argument practically constant over a certain period of time, that is, the altitude will always be the same in a certain geographic point of the planet. However, space vehicles with missions defined for this type of orbit suffer perturbations that can change the orbital configuration. Thus, it is necessary to adopt orbital corrections through orbital maneuvers in order to maintain the characteristics of the proposed mission. In the considered model, the space vehicle suffers perturbations due to the gravitational potential of the Sun and the non-sphericity of Venus. To make the orbital corrections in the orbit of the vehicle, continuous low-thrust control strategies will be adopted in order to find the regions where the lowest energy demanding for the mission occurs. Finally, for this analysis, orbital conditions relevant to the problem were selected, such as critical conditions, regions close to the surface and close to the region where the perturbative forces are equivalent. KEYWORDS: Low-thust control. Space vehicle. Orbital maneuvers. Venus. LISTA DE ILUSTRAÇÕES Figura 2.1 Elementos que caracterizam uma órbita no espaço: semieixo maior a, excentrici- dade e, inclinação i, argumento de periapsis ω, longitude do nodo ascendente Ω e anomalia verdadeira θ. Em geral, nos cálculos orbitais ao invés de θ, utiliza-se a anomalia média M . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18 Figura 2.2 Um sistema de referência inercial composto por três corpos. Vênus é o corpo central e possui massa MV , o veículo espacial de massa m orbita Vênus e o Sol com massa MSOL é o corpo perturbador. Os vetores posição R⃗V , R⃗ e R⃗Sol são dados em relação à origem inercial O do três corpos. Os r⃗ e r⃗Sol representam, respectivamente, os vetores posição do veículo espacial e do Sol em relação ao corpo central. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20 Figura 2.3 Condições iniciais de órbitas congeladas em torno de Vênus. . . . . . . . . . . 28 Figura 2.4 Direção das acelerações radial (S), tangencial (T ) e transversal (W ). . . . . . . 28 Figura 2.5 Modelo esquemático dos quatro métodos de controle de órbitas congeladas. . . 31 Figura 3.1 Fluxograma do método para o cálculo da aceleração de controle. . . . . . . . . 39 Figura 4.1 Variação da aceleração de controle, para os quatro métodos de análise, em função de a, e e i em (a) Marte ( resultados na coluna à esquerda) e (b) Vênus (resultados na coluna à direita). As condições iniciais da órbita são: a = 10.000 km, i = 30◦, e = 0, 2 , ω = 0◦, Ω = 0◦ e M = 0◦. . . . . . . . . . . . . . . . . . . . . . . . 41 Figura 4.2 Superfície da aceleração de controle em m/s2 para o método low-thrust e função de (a) a e i, (b) e e i, (c) e e a. . . . . . . . . . . . . . . . . . . . . . . . . . . 43 Figura 4.3 Evolução orbital de ω em 30 períodos orbitais de Vênus e i = 35◦ conside- rando (a) os coeficientes dos harmônicos esféricos até ordem e grau 15 e (b) os coeficientes dos harmônicos zonais J2, J3 e J4. . . . . . . . . . . . . . . . . . 46 Figura 4.4 Evolução orbital de ω em 30 períodos orbitais de Vênus e i = 45◦ para (a) os coeficientes dos harmônicos esféricos até ordem e grau 15 e (b) os coeficientes dos harmônicos zonais J2, J3 e J4. . . . . . . . . . . . . . . . . . . . . . . . . 46 Figura 4.5 Evolução orbital de e em 30 períodos orbitais de Vênus e i = 35◦ considerando (a) os coeficientes dos harmônicos esféricos até ordem e grau 15 e (b) coeficientes dos harmônicos zonais J2, J3 e J4. . . . . . . . . . . . . . . . . . . . . . . . . 48 Figura 4.6 Evolução orbital de e em 30 períodos orbitais de Vênus e i = 45◦ considerando (a) coeficientes dos harmônicos esféricos até ordem e grau 15 e (b) os coeficientes dos harmônicos zonais J2, J3 e J4. . . . . . . . . . . . . . . . . . . . . . . . . 48 Figura 4.7 Evolução orbital de ω em 30 períodos orbitais de Vênus e i = 62◦ para (a) harmônicos até ordem e grau 15 e (b) harmônicos zonais. . . . . . . . . . . . . 50 Figura 4.8 Evolução orbital de e em 30 períodos orbitais de Vênus e i = 62◦ considerando (a) coeficientes dos harmônicos esféricos até ordem e grau 15 e (b) os coeficientes dos harmônicos zonais J2, J3 e J4. . . . . . . . . . . . . . . . . . . . . . . . . 51 Figura 4.9 Evolução orbital de ω em 30 períodos orbitais de Vênus e i = 75◦ para (a) coeficientes dos harmônicos esféricos até ordem e grau 15 e (b) os coeficientes dos harmônicos zonais J2, J3 e J4. . . . . . . . . . . . . . . . . . . . . . . . . 53 Figura 4.10Evolução orbital de ω em 30 períodos orbitais de Vênus e i = 85◦ considerando (a) coeficientes dos harmônicos esféricos até grau e ordem 15 (b) os coeficientes dos harmônicos zonais J2, J3 e J4. . . . . . . . . . . . . . . . . . . . . . . . . 53 Figura 4.11Evolução orbital de e em 30 períodos orbitais de Vênus e i = 75◦ considerando (a) os coeficientes harmônicos esféricos até grau e ordem 15 e (b) os coeficientes dos harmônicos zonais J2, J3 e J4. . . . . . . . . . . . . . . . . . . . . . . . . 54 Figura 4.12Evolução orbital de e em 30 períodos orbitais de Vênus e i = 85◦ considerando (a) os coeficientes dos harmônicos esféricos até grau e ordem grau 15 e (b) os coeficientes dos harmônicos zonais J2, J3 e J4. . . . . . . . . . . . . . . . . . 55 LISTA DE TABELAS Tabela 2.1 – Taxa de variação média dos elementos orbitais por componente da aceleração perturbadora. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30 Tabela 3.1 – Elementos orbitais e parâmetros dos planetas Vênus e Marte. . . . . . . . . . . . 36 Tabela 3.2 – Variação de velocidade gerada no veículo devido a ação perturbativa do potencial gravitacional de Vênus (PIg) e da presença do Sol (PISol). . . . . . . . . . . . . 37 Tabela 3.3 – Comparação entre o resultado obtido por Wu, Jiang e Li (2014) e pelo programa implementado para o caso da Terra considerando J2, J3 e J4. . . . . . . . . . . . 38 Tabela 3.4 – Comparação entre o resultado obtido por Wu, Jiang e Li (2014) e pelo programa implementado para o caso de Marte considerando J2, J3 e J4. . . . . . . . . . . . 38 Tabela 4.1 – Condições iniciais das órbitas congeladas selecionadas para o estudo. . . . . . . . 45 Tabela 4.2 – Aplicação do método low-thrust para as órbitas da região de inclinação abaixo da crítica, considerando ∆ω = 3◦ e ∆ω = 5◦. . . . . . . . . . . . . . . . . . . . . . 47 Tabela 4.3 – Aplicação do método low-thrust para as órbitas da região de inclinação abaixo da crítica, considerando ∆e = 0, 001 e ∆e = 0, 01. . . . . . . . . . . . . . . . . . . 49 Tabela 4.4 – Aplicação do método low-thrust para a órbita da região de inclinação crítica, considerando ∆ω = 3◦ e ∆ω = 5◦. . . . . . . . . . . . . . . . . . . . . . . . . . 50 Tabela 4.5 – Aplicação do método low-thrust para as órbitas da região de inclinação crítica, considerando ∆e = 0, 001 e ∆e = 0, 01. . . . . . . . . . . . . . . . . . . . . . . 51 Tabela 4.6 – Aplicação do método low-thrust para as órbitas da região de inclinação polar, considerando ∆ω = 3◦ e ∆ω = 5◦. . . . . . . . . . . . . . . . . . . . . . . . . . 54 Tabela 4.7 – Aplicação do método low-thrust para as órbitas da região de inclinação polar, considerando ∆e = 0, 001 e ∆e = 0, 01. . . . . . . . . . . . . . . . . . . . . . . 55 SUMÁRIO 1 INTRODUÇÃO . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 2 FUNDAMENTAÇÃO TEÓRICA . . . . . . . . . . . . . . . . . . . . . . . . . 18 2.1 A órbita no espaço . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18 2.2 Problema elíptico restrito de três corpos . . . . . . . . . . . . . . . . . . . . . . . 19 2.2.1 Formulação matemática . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20 2.3 Modelo dinâmico . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 2.3.1 Aceleração P⃗G . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22 2.3.2 Coeficientes dos Harmônicos esféricos de Vênus . . . . . . . . . . . . . . . . . . 24 2.4 Integral da aceleração . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24 2.5 Órbitas congeladas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25 2.6 Condições iniciais . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 2.7 Estratégias de controle low-thrust . . . . . . . . . . . . . . . . . . . . . . . . . . 28 2.7.1 Estratégia de controle: método 1 . . . . . . . . . . . . . . . . . . . . . . . . . . 31 2.7.2 Estratégia de controle: método 2 . . . . . . . . . . . . . . . . . . . . . . . . . . 32 2.7.3 Estratégia de controle: método 3 . . . . . . . . . . . . . . . . . . . . . . . . . . 33 2.7.4 Estratégia de controle: método 4 . . . . . . . . . . . . . . . . . . . . . . . . . . 33 3 METODOLOGIA . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35 3.1 Integração numérica da equação de movimento . . . . . . . . . . . . . . . . . . . 35 3.2 Grau e ordem dos coeficientes dos harmônicos esféricos . . . . . . . . . . . . . . 35 3.3 Parâmetros de Vênus e Marte . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36 3.4 Definição da região de estudo . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36 3.5 Implementação numérica do método de controle . . . . . . . . . . . . . . . . . . 37 4 RESULTADOS E DISCUSSÕES . . . . . . . . . . . . . . . . . . . . . . . . . 40 4.1 Efeito dos elementos orbitais nos métodos de controle . . . . . . . . . . . . . . . 40 4.2 Avaliação das órbitas congeladas nas regiões de interesse ao redor de Vênus . . . 44 4.2.1 Região de inclinação abaixo da crítica . . . . . . . . . . . . . . . . . . . . . . . . 45 4.2.1.1 Limitando ω . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45 4.2.1.2 Limitando e . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47 4.2.2 Região de inclinação crítica . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49 4.2.2.1 Limitando ω . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50 4.2.2.2 Limitando e . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51 4.2.3 Região de inclinação polar . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52 4.2.3.1 Limitando ω . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52 4.2.3.2 Limitando e . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54 5 CONCLUSÃO . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56 REFERÊNCIAS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58 15 1 INTRODUÇÃO A exploração espacial representa um marco técnico e científico na história da humanidade, pois contribuiu para um melhor entendimento do Sistema Solar e de seus corpos celestes, bem como, no desenvolvimento de tecnologias de comunicação, monitoramento e vigilância da Terra. O início dessa exploração remonta ao lançamento do satélite soviético Sputinik 1, em 1957, que foi o primeiro satélite artificial a entrar em órbita (CARLEIAL, 1999). Nos anos seguintes, a exploração espacial se intensificou culminando em missões tripuladas, como o voo orbital de Yuri Gagarin, em 1961, e as missões Apollo, entre 1968 e 1972, que levaram o homem à Lua. Contudo, uma das maiores realizações dos programas espaciais foi a possibilidade da exploração do Sistema Solar por meio de veículos não tripulados. O primeiro corpo celeste visitado por uma missão espacial não tripulada foi a Lua, em 1959, através das missões soviéticas Luna 1, 2 e 3 (WILLIAMS, 2005), e da missão americana Pioneer 4 (WILLIAMS, 2021). Já o primeiro planeta a ser visitado foi Vênus, através da missão americana Mariner 2 que, em 1962, fez uma passagem bem-sucedida pelo planeta (SONETT, 1963). Nas décadas seguintes, outros corpos do Sistema Solar foram visitados por diversas missões e, em 2012, atingiu-se a fronteira interestelar com a nave Voyager 1, o mesmo ocorreu com a Voyager 2 em 2018 (NASA, 2021a). Essas missões permitiram o acesso às informações que contribuíram de maneira significativas para entendermos melhor o processo de formação do Sistema Solar, bem como, a interação entre os corpos que o compõem, entretanto, ainda há muitas perguntas a serem respondidas sobre os corpos celestes e suas origens, e isso motiva o planejamento e a implantação de novas missões. Muitas dessas questões são feitas em relação a Vênus que, apesar de ter sido o primeiro planeta a ser visitado, ainda é desconhecido em diversos aspectos. Vênus é o segundo planeta rochoso mais próximo ao Sol e possui grandes semelhanças com a Terra em relação à massa, volume, tamanho e composição química, o que fez com que Vênus fosse conhecido como planeta quase gêmeo da Terra (DE LUCA, 1982). Conforme apresentado por NASA (2021b), Vênus não apresenta satélites naturais, seu período orbital é de aproximadamente 224,7 dias, sua atmosfera é composta majoritariamente por dióxido de carbono e, em sua superfície, a temperatura média é de 464◦C e a pressão média é de 90,8 atm, o que dificulta a implantação de missões de exploração da superfície do planeta. Sendo assim, Vênus foi foco de missões interplanetárias apenas no início da exploração espacial. Entretanto, nos últimos anos, com os avanços tecnológicos dos veículos espaciais e a descoberta de possíveis evidências de vida em sua atmosfera, elevou-se o interesse por esse corpo celeste, de modo que, diversas missões estão sendo planejadas para o planeta, como é o caso da missão DAVINCI+, que pretende coletar dados atmosféricos (KIEFER et al., 2020), e da missão VERITAS, que pretende estudar a topografia e a gravidade de Vênus (SMREKAR et al., 2020). No planejamento de uma missão espacial, dentre o múltiplos fatores que devem ser avaliados, destaca-se o estudo orbital, que garante a implementação do modelo de órbita ideal para cada missão. Dentre os diferentes modelos orbitais, alguns se destacam para aplicações em missões de estudo da at- 16 mosfera e da superfície de planetas rochosos, pois estabelecem condições adequadas de monitoramento contínuo do planeta, essas órbitas são conhecidas por órbitas heliossíncronas e órbitas congeladas. As órbitas heliossíncronas se aproximam de órbitas polares, ou seja, apresentam inclinações próximas a 90◦ com o plano orbital fixado em direção ao Sol, resultando em órbitas localizadas em regiões de média a baixa altitude com excentricidade quase circular (MACDONALD et al., 2010). Por sua vez, as órbitas congeladas apresentam valores de excentricidade e argumento de periapsis que permanecem praticamente constantes ao longo do tempo, fazendo com que a órbita mantenha a altitude em deter- minado ponto geográfico do planeta (HOOGLAND, 2015). Wu, Jiang e Li (2015) mostraram que as órbitas congeladas e heliossíncronas naturais se assemelham às artificiais fazendo com que as missões a planetas rochosos sejam beneficiadas pela aplicação desses tipos de órbitas. Na Terra, a órbita heliossíncrona é aplicada em missões como a dos satélites meteorológicos da série NOAA (National Oceanic and Atmosphere Administration, dos Estados Unidos) (DAVIS, 2007) e o satélite Sino-Brasileiro de observação da superfície terrestre CBERS (China-Brazil Earth-Resources Satellite) (SANTANA; COELHO, 2009) e, são exemplos de missões que empregam órbita congelada o satélite europeu de sensoriamento remoto ERS-2 (European Remote Sensing Satellite) (HOOGLAND, 2015), e a missão conjunta entre França e Estados Unidos de mapeamento da topografia oceânica TOPEX/Poseidon (FU et al., 1994). Em particular, neste trabalho estamos interessados nas órbitas congeladas ao redor de Vênus. De acordo com Schulz (1996), essa órbita é caracterizada por apresentar uma altitude praticamente constante em determinado ponto da superfície do planeta, isso se dá a partir da correta seleção dos parâmetros orbitais iniciais da órbita, especificamente, os valores de excentricidade (e), argumento do periapsis (ω), semieixo maior (a) e inclinação (i). Entretanto, em missões de monitoramento, mapeamento e imageamento, os valores de a e i são previamente definidos, uma vez que, esses valores devem satisfazer as características dos equipamentos embarcados na missão. Dessa forma, pode-se encontrar os valores de e e ω que satisfaçam as condições de órbitas congeladas. Tipicamente, essas órbitas apresentam pequenos valores de e (quase circulares) e a variação de e e ω da órbita inicial é praticamente constante ao longo do tempo, mesmo sob a ação de forças perturbativas (HOOGLAND, 2015). Sendo assim, missões de monitoramento, mapeamento e imageamento são beneficiadas pela manutenção da altitude sobre a área de estudo, pois permitem a criação de uma série de dados históricos a partir de um mesmo ponto de observação. Contudo, comparando a utilização dessas órbitas na Terra e em Vênus nota-se que há diferenças relevantes, principalmente com relação às perturbações sofridas pelo veículo espacial ao redor desses planetas. Como apresentado na literatura, por exemplo Roy (1998), as forças perturbadoras mais significativas sobre veículos espaciais são causadas pelo efeito da não esfericidade do corpo orbitado, pela pressão de radiação, pelo efeito do arrasto atmosférico e pela presença de um terceiro corpo (satélites naturais e o Sol). Há também, o efeito perturbativo de menor impacto causado por forças eletromagnéticas, precessão relativística e albedo (KUGA; CARRARA; RAO, 2011). Comparando as perturbações sofridas por um veículo espacial ao redor da Terra e de Vênus, pode-se ressaltar que um veículo espacial ao redor de Vênus não sofre perturbação devido a um satélite natural, além disso, Vênus possui uma não esfericidade mais pronunciada que a Terra e está muito mais próximo do Sol. 17 As perturbações que acometem os veículos espaciais fazem com que suas órbitas sejam alteradas ao longo do tempo. Dessa forma, faz-se necessário a aplicação de manobras de correção orbital que utilizam um sistema propulsivo fixo no veículo, responsável por modificar sua velocidade e, consequentemente, sua órbita (DOMINGOS; PRADO; GOMES, 2014). As manobras de correção orbital podem ser implementadas a partir de diferentes estratégias de controle, como por exemplo, o controle multi-impulsivo, ou também, o controle contínuo low-thrust (WU; JIANG; LI, 2014) e, normalmente, são empregados sistemas de propulsão elétricos que garantem uma maior força propulsiva quando comparados aos sistemas de combustão química (MORANTE; RIVO; SOLER, 2021). Dentro desse contexto, será estudada a estabilidade de órbitas de veículos espaciais em órbitas congeladas ao redor do planeta Vênus e as aplicações de manobras de correção orbital utilizando o método de controle contínuo low-thrust em órbitas congeladas ao redor desse planeta. Para isso, será avaliado o efeito da perturbação do potencial gravitacional do Sol e da não esfericidade de Vênus sobre as órbitas. O objetivo é selecionar as condições orbitais do veículo para a aplicação da correção orbital que resulte no mínimo gasto de combustível possível. Para avaliar o comportamento da órbita dos veículos espaciais e inferir a magnitude da variação da aceleração média devida às forças de perturbação, será empregado o programa numérico de Sanchez, Prado e Yokoyama (2014) que foi adaptado para o caso de Vênus. O modelo de correção para órbitas congeladas que será abordado aqui é baseado no trabalho de Wu, Jiang e Li (2014), onde são descritas as estratégias para obter o controle de menor consumo energético de um veículo espacial ao redor de um planeta. Dessa forma, é possível aplicar a manobra de uma maneira mais eficiente de forma que seja aumentada a vida útil da missão. A organização desse trabalho foi feita da seguinte forma: na seção 2 é apresentada uma fundamen- tação teórica com algumas definições e conceitos básicos necessários para o entendimento do trabalho; na seção 3 é apresentada a metodologia utilizada neste estudo; na seção 4, os resultados são discutidos detalhadamente; finalmente, na seção 5 são apresentadas as considerações finais. 18 2 FUNDAMENTAÇÃO TEÓRICA Neste capítulo são apresentadas as definições e os conceitos teóricos empregados no trabalho. Apresenta-se, inicialmente, o modelo orbital espacial, os elementos orbitais e parâmetros fundamentais para o estudo. Em seguida, introduz-se a obtenção da equação de movimento a partir do problema elíptico restrito de três corpos. Além do modelo dinâmico, apresentam-se as perturbações consideradas, a formulação matemática e o cálculo da contribuição de cada perturbação a partir do modelo da integral da aceleração. Além da parte dinâmica do modelo, é apresentado o conceito de órbita congelada, com o equacio- namento matemático e as características das condições inicias para esse modelo orbital de veículos espaciais ao redor de Vênus. Por fim, apresentam-se as estratégias de controle por empuxo contínuo empregadas para corrigir a trajetória do veículo espacial. 2.1 A ÓRBITA NO ESPAÇO Nos estudos astronômicos e astrodinâmicos, a órbita de um corpo celeste ou de um veículo espacial pode ser definida a partir de seis parâmetros que são chamados de elementos orbitais. Esses elementos definem a orientação da órbita em relação aos eixos do sistema, sua forma e a posição do corpo nessa órbita em um determinado instante de tempo (ROY, 1998). A Figura 2.1 apresenta a órbita no espaço e os elementos orbitais. Figura 2.1 – Elementos que caracterizam uma órbita no espaço: semieixo maior a, excentricidade e, inclinação i, argumento de periapsis ω, longitude do nodo ascendente Ω e anomalia verdadeira θ. Em geral, nos cálculos orbitais ao invés de θ, utiliza-se a anomalia média M . Fonte: elaborado pelo próprio autor. Segundo a literatura, por exemplo Roy (1998) e Murray e Dermott (1999), os elementos orbitais que definem a orientação da órbita de um corpo são: a inclinação (i) a longitude do nodo ascendente (Ω) 19 e o argumento do periapsis (ω). A inclinação se refere ao ângulo formado entre o plano de referência do sistema e o plano orbital do corpo, podendo assumir valores entre 0 ≤ i ≤ 180◦. Segundo Murray e Dermott (1999), quando i < 90◦ diz-se que o movimento é prógrado, enquanto que para i ≥ 90◦ o movimento é dito retrógrado. Por sua vez, a longitude do nodo ascendente compreende o ângulo formado entre o eixo de referência x̂ do sistema e o ponto de interseção do plano equatorial com a órbita (nó ascendente). Trata-se de um parâmetro positivo, que pode tomar valores entre 0 ≤ Ω ≤ 360◦. Por fim, o argumento do periapsis é o ângulo medido entre o nó ascendente e o perigeu da órbita, ou seja, o ponto mais próximo da órbita ao corpo orbitado. Assim como Ω, o argumento do periapsis pode assumir valores entre 0 ≤ ω ≤ 360◦. Em seguida, os elementos que caracterizam o formato e o tamanho da órbita são: a excentricidade (e) e o semieixo maior (a). A excentricidade refere-se a medida relacionada à distância entre o foco e o centro de uma órbita. A partir desse parâmetro pode-se classificar as órbitas como: circular, quando e = 0; elíptica, quando 0 < e < 1; parabólica, quando e = 1; e hiperbólica quando e ≥ 1. Para o caso de órbita elíptica, tem-se que o semieixo maior corresponde a metade da distância entre o perigeu e o apogeu da órbita (ponto mais distante da órbita ao corpo orbitado). Por último, para caracterizar a posição de um corpo em uma órbita em determinado instante de tempo emprega-se a anomalia verdadeira θ, que corresponde ao ângulo formado entre o perigeu e o corpo orbitador. Normalmente, nos cálculos orbitais, em vez de θ, costuma-se empregar a anomalia média (M ), que corresponde a uma função monotonicamente crescente da anomalia verdadeira (CURTIS, 2014). Matematicamente, M pode ser calculada por: M = 2π T t, (2.1) em que T corresponde ao período da órbita e t ao tempo decorrido desde a passagem pelo perigeu. A razão 2π/T é denominada velocidade angular média (n). 2.2 PROBLEMA ELÍPTICO RESTRITO DE TRÊS CORPOS O problema restrito de três corpos é uma simplificação do problema geral de três corpos. Original- mente, o problema geral considera que um dos corpos do sistema possui massa muito menor do que a dos outros dois corpos. Com isso, os corpos de maior massa se movimentam em órbitas Keplerianas circulares sem influência do corpo de massa menor (CURTIS, 2014). São exemplos de situações em que esse tipo de formulação é empregada: o movimento de um veículo espacial em um sistema de um planeta e sua lua, um sistema composto por um planeta, asteroide e o Sol Entretanto, o caso circular nem sempre é representativo dos sistemas celestes reais, portanto, adaptou-se o problema restrito clássico para o caso em que as órbitas dos corpos massivos que se movimentam em torno do centro de gravidade do sistema são elípticas, sendo conhecido como problema elíptico restrito de três corpos (MACEDO; JUNIOR, 2018). A principal diferença entre os dois modelos é que no caso elíptico os vetores de posição dos corpos do sistema não são mais constantes. 20 2.2.1 Formulação matemática No contexto deste trabalho, o problema elíptico restrito de três corpos emprega uma abordagem clássica, como pode ser visto em Casteletti (2013), que resulta nas equações de movimento dos corpos em um sistema de referência inercial, onde tem-se a função perturbadora do terceiro corpo, o Sol. Considerando um sistema de três corpos composto por: Vênus, como corpo central, um veículo espacial que orbita Vênus e o Sol, sendo suas massas MV , m e MSol, respectivamente, tem-se os vetores posição R⃗V , R⃗ e R⃗Sol em relação à origem inercial O do três corpos. Além disso, r⃗ e r⃗Sol representam, respectivamente, os vetores posição do veículo espacial e do Sol em relação ao corpo central. A Figura 2.2 apresenta este sistema de três corpos. Figura 2.2 – Um sistema de referência inercial composto por três corpos. Vênus é o corpo central e possui massa MV , o veículo espacial de massa m orbita Vênus e o Sol com massa MSOL é o corpo perturbador. Os vetores posição R⃗V , R⃗ e R⃗Sol são dados em relação à origem inercial O do três corpos. Os r⃗ e r⃗Sol representam, respectivamente, os vetores posição do veículo espacial e do Sol em relação ao corpo central. Fonte: elaborado pelo próprio autor. O sistema de coordenadas cartesianas tem origem no corpo central e, sendo assim, os vetores posição do veículo espacial e do Sol têm coordenadas r(x, y, z) e rSol(xsol, ySol, zSol), respectivamente. As equações de movimento dos três corpos no sistema de referência inercial são calculadas a partir das equações de movimento e de gravitação geral de Newton: MV ¨⃗ RV = GmMV r⃗ r3 +GMVMSol r⃗Sol r3Sol , (2.2) m ¨⃗ R = GmMSol r⃗Sol − r⃗ |r⃗Sol − r⃗|3 −GmMV r⃗ r3 , (2.3) MSol ¨⃗ RSol = GmMSol r⃗ − r⃗Sol |r⃗ − r⃗Sol|3 −GMVMSol r⃗Sol r3Sol , (2.4) em que G é a constante de gravitação universal, r e rSol representam o módulo dos vetores posição r⃗ e r⃗Sol. 21 O cálculo das acelerações dos corpos secundários em relação ao corpo central é realizado por: ¨⃗r = ¨⃗ R− ¨⃗ RV , (2.5) ¨⃗rSol = ¨⃗ RSol − ¨⃗ RV . (2.6) Dessa forma, as equações de movimento relativos são encontradas a partir da substituição das relações (2.2) a (2.4) nas equações das acelerações relativas, eqs. (2.5) e (2.6). ¨⃗r = −G(MV +m) r⃗ r3 +GMSol ( r⃗Sol − r⃗ |r⃗Sol − r⃗|3 − r⃗Sol r3Sol ) , (2.7) ¨⃗rSol = −G(MV +MSol) r⃗ r3 +Gm ( r⃗ − r⃗Sol |r⃗ − r⃗Sol|3 − r⃗ r3 ) . (2.8) Considerando que m << MV e MSol e, portanto, o veículo espacial não influencia o movimento do corpo central e do Sol, pode-se desconsiderar os termos que envolve m de modo que se obtém as equações de movimento do veículo espacial eq. (2.9) e do Sol eq. (2.10): ¨⃗r = −GMV r⃗ r3 +GMSol ( r⃗Sol − r⃗ |r⃗Sol − r⃗|3 − r⃗Sol r3Sol ) , (2.9) ¨⃗rSol = −G(MV +MSol) r⃗ r3 . (2.10) 2.3 MODELO DINÂMICO Neste trabalho, estamos interessados no estudo do movimento de veículos espaciais ao redor de Vênus. Para isso, emprega-se o modelo planetocêntrico, que descreve um sistema de multicorpos centrado no planeta orbitado. Especificamente neste caso, considera-se um sistema de três corpos composto por: Vênus, como o corpo central, orbitado pelo Sol e pelo veículo espacial.Kuga, Carrara e Rao (2011) apresenta que as principais perturbações atuantes nesse tipo de sistema dinâmico são o potencial gravitacional gerado pela não esfericidade do corpo central, a presença do terceiro corpo (Sol) e a pressão de radiação. Entretanto, Pina (2019) mostrou que para o caso de Vênus, a pressão de radiação não influencia de maneira significativa o comportamento do veículo espacial e, portanto, pode ser desconsiderada para este estudo. Como apresentado em Sanchez, Prado e Yokoyama (2014), tomando como sistema de referência o plano equatorial de Vênus, a equação de movimento do veículo espacial, considerando as perturbações do sistema, é escrita pela eq. (2.11). ¨⃗r = −GMV r3 r⃗ +GMSol ( r⃗Sol − r⃗ |r⃗Sol − r⃗|3 − r⃗Sol r3Sol ) + P⃗G, (2.11) onde P⃗G corresponde à aceleração devido ao potencial gravitacional gerado pela não-esfericidade de Vênus. A órbita descrita pelo Sol em torno do centro do sistema é considerada Kepleriana (SANCHEZ; PRADO; YOKOYAMA, 2014). 22 2.3.1 Aceleração P⃗G A aceleração devido ao potencial gravitacional de Vênus pode ser calculada a partir do modelo potencial recursivo modificado que foi apresentado por Sanchez, Prado e Yokoyama (2014). Pelo fato de Vênus não possuir distribuição de massa homogênea, o movimento do veículo sofrerá alterações, causadas pelo potencial gravitacional do planeta. A aceleração devido ao potencial gravitacional depende da massa e dos valores dos coeficientes dos harmônicos esféricos ligados à geometria do planeta. De acordo com a literatura, por exemplo Kuga, Carrara e Rao (2011), para encontrar a função do potencial gravitacional parte-se do conceito de potencial de força dado por: U = GM r , (2.12) em que M representa a massa que gera o potencial e r é a distância entre as duas massas do sistema. Considerando agora dois pontos dados por suas coordenadas esféricas P (r, θ, λ) e P ′(r′, θ′, λ′), a distância entre eles (∆) será ∆2 = r2 + r′2 − 2rr′ cosψ, (2.13) onde ψ é o ângulo formado entre r e r′. Em seguida, supondo que r′ < r e desenvolvendo a série de potências em relação a r′ r , obtém-se: 1 ∆ = 1 r { P0(cosψ) + r′ r P1(cosψ) + ( r′ r )2 P2(cosψ) + ... } = 1 r ∞∑ n=0 ( r′ r )n Pn(cosψ), (2.14) em que os termos Pn(cosψ) para n = 0, 1, 2, ...,∞ são conhecidos como polinômios de Legendre. Considerando o potencial de força gerado por um elemento infinitesimal de massa e integrando em toda a distribuição de massa do corpo, obtém-se o potencial gravitacional de Vênus a partir da eq. (2.15): U(r, θ, λ) = GMV r + GMV r M∑ n=2 n∑ m=0 ( Re r )n × ( C̄nm cosmλ+ S̄nm sinmλ ) P̄nm(cos θ), (2.15) onde Re é o raio equatorial de Vênus, P̄nm são polinômios associados de Legendre normalizados, C̄nm e S̄nm são os coeficientes dos harmônicos esféricos normalizados, os valores de n e m representam o grau e a ordem dos coeficientes e M é o maior valor de grau e ordem da expansão dos harmônicos esféricos. De acordo com Kuga, Carrara e Rao (2011), os coeficientes harmônicos C̄nm e S̄nm são valores constantes que dependem da geometria do corpo e apresentam as seguintes propriedades: (i) Os termos para n = m = 0 descrevem o potencial gravitacional principal do planeta, pois C00 = 1 e S00 = 0; (ii) Como a origem do sistema de coordenadas coincide com o centro de massa de Vênus, então pode-se mostrar que C10 = C11 = S11 = Sn0 = 0 e; 23 (iii) É comum chamar os harmônicos esféricos Cn0 de coeficientes zonais Jn, onde Cn0 = Jn. Como mostrado em Sanchez, Prado e Yokoyama (2014), para facilitar a manipulação da equação de potencial e possibilitar o cálculo recursivo, é possível rescrever a eq. (2.15) da seguinte forma: U(r, θ, λ) = GMV r + U∗(r, θ, λ), (2.16) onde U∗(r, θ, λ) = GMV r M∑ m=0 {[ M∑ n=2 ( Re r )n C̄nmP̄nm(cos θ) ] cosmλ + [ M∑ n=2 ( Re r )n S̄nmP̄nm(cos θ) ] sinmλ } . (2.17) Definindo P⃗G em termos do gradiente de U∗(r, θ, λ), a expressão para P⃗G pode ser descrita em coordenadas esféricas como: P⃗G = ∇U∗(r, θ, λ) = ∂U∗ ∂r ∂r⃗/∂r |∂r⃗/∂r| + 1 r ∂U∗ ∂θ ∂r⃗/∂θ |∂r⃗/∂θ| + 1 r sin θ ∂U∗ ∂λ ∂r⃗/∂λ |∂r⃗/∂λ| . (2.18) Entretanto, para realizar as integrações numéricas, o modelo dinâmico considera o sistema de coordenadas cartesianas. Assim, as componentes de P⃗G geradas pelos harmônicos esféricos de grau n e ordem m são descritas de acordo com as eqs. (2.19), (2.20) e (2.21): PGm,n,x = ∂U∗ ∂r sin θ cosλ+ 1 r ∂U∗ ∂θ cos θ cosλ− 1 r sin θ ∂U∗ ∂λ sinλ, (2.19) PGm,n,y = ∂U∗ ∂r sin θ sinλ+ 1 r ∂U∗ ∂θ cos θ sinλ+ 1 r sin θ ∂U∗ ∂λ cosλ, (2.20) PGm,n,z = ∂U∗ ∂r cos θ − 1 r ∂U∗ ∂θ sin θ, (2.21) onde ∂U∗ ∂r = −GMV r2 (n+ 1) ( Re r )n [ C̄nmP̄nm(cos θ) cosmλ+ S̄nmP̄nm(cos θ) sinmλ ] , (2.22) ∂U∗ ∂θ = GMV r ( Re r )n [ C̄nmP̄ ′ nm(cos θ) cosmλ+ S̄nmP̄ ′ nm(cos θ) sinmλ ] , (2.23) ∂U∗ ∂λ = −GMV r ( Re r )n m [ C̄nmP̄nm(cos θ) sinmλ− S̄nmP̄nm(cos θ) cosmλ ] , (2.24) em que P̄ ′ nm é a primeira derivada das funções associadas de Legendre normalizadas. Desta forma, a aceleração P⃗G pode ser descrita através dos somatórios das parcelas geradas por PGm,n,x , PGm,n,y e PGm,n,z na direção de cada um dos versores. É válido notar que partindo das expressões obtidas anteriormente, é possível descrever a aceleração individual gerada por cada termo 24 do geopotencial, o que possibilita a realização de uma análise individualizada da contribuição de cada harmônico. A contribuição total da aceleração gerada pela não esfericidade pode ser descrita de acordo com a equação: P⃗G = M∑ m=0 M∑ n=2 PGm,n,x î+ M∑ m=0 M∑ n=2 PGm,n,y ĵ + M∑ m=0 M∑ n=2 PGm,n,z k̂. (2.25) 2.3.2 Coeficientes dos Harmônicos esféricos de Vênus Para o caso de Vênus, são conhecidos os coeficientes dos harmônicos esféricos até grau e ordem 180, conforme apresentado em Konopliv, Banerdt e Sjogren (1999). Esses dados foram obtidos pela missão Magellan, lançada pela Nasa em 1989 com o objetivo de estudar o planeta Vênus quanto à sua topografia, atividade da superfície e mapeamento de seu campo gravitacional. No contexto dos coeficientes dos harmônicos esféricos, é de extrema importância a avaliação coeficientes zonais que, como explicado anteriormente, compreendem os coeficientes dos harmônicos esféricos do tipo Cn0, que são conhecidos por Jn. Isso se deve ao fato de que diversas metodologias simplificam as equações de movimento considerando apenas os coeficientes dos harmônicos zonais, já que, são os que apresentam maior valor numérico. Dentre os planetas rochosos do Sistema Solar, Vênus é o único que apresenta os coeficientes dos harmônicos zonais J2, J3 e J4 de mesma ordem de grandeza. Ao contrário do caso da Terra e de Marte, por exemplo, em que o coeficiente do harmônico zonal J2 é dominante aos demais. Isso implica que, para metodologias simplificadas de análise de Vênus, é interessante empregar os harmônicos zonais até J4. De acordo com Anderson, Macdonald e Yen (2014), os valores desses harmônicos em Vênus são: J2 = 4, 4580× 10−6, J3 = −2, 1082× 10−6 e J4 = −2, 1471× 10−6. 2.4 INTEGRAL DA ACELERAÇÃO O método da integral da aceleração é uma importante ferramenta para avaliar o efeito das perturba- ções sobre o movimento de um veículo espacial. Esse método consiste no cálculo da influência das perturbações, em termos da variação total de velocidade provocadas pela atuação desses perturbadores sobre o veículo (SANCHEZ; PRADO, 2017). Como pode ser visto em Prado (2013), o método da integral da aceleração considera um sistema planetocêntrico e o potencial de perturbação (V ) que atua sobre o veículo espacial. Dessa forma, a integral da magnitude da força em um período orbital do veículo (N ) é dada a partir da integral do módulo do gradiente do potencial de perturbação, que representa a força aplicada no veículo espacial devido ao agente perturbador. PI = ∫ N 0 |∇V |dt (2.26) Prado (2013) mostrou que a partir da eq. (2.26) é possível obter mais dois tipos de integral da aceleração, considerando mudanças em termos da anomalia média e da anomalia excêntrica. Sanchez e Prado (2017) discutem as diferentes aplicações dos três modelos principais da integral da aceleração. 25 O primeiro modelo emprega o valor absoluto da aceleração, ou seja, com esse método é possível encontrar a variação total da velocidade provocada, individualmente, por cada um dos perturbadores, sem que haja a necessidade de desconsiderar as demais perturbações. O segundo modelo, apresenta o efeito direto da variação de energia sobre o veículo espacial, avaliando se ocorreu uma transferência de energia do sistema para o veículo ou se foi o oposto, ou seja, o modelo é uma boa escolha quando se quer avaliar a variação de energia gerada pelo perturbador. Por fim, o terceiro modelo considera a integral devido a cada uma das componentes da aceleração, que se torna uma ótima opção para comparar órbitas que sofrem maior ou menor perturbação dentro de uma região semelhante (SANCHEZ; PRADO, 2017). Neste trabalho emprega-se o primeiro modelo da integral da aceleração, uma vez que, o objetivo de utilizar este método é comparar de maneira absoluta, a perturbação ocasionada por cada um dos perturbadores do sistema. Desta forma, o método da integral da aceleração é utilizado considerando as perturbações geradas pelo potencial gravitacional do Sol e pela não esfericidade de Vênus. Como apresentado por Sanchez e Prado (2017), considerando o tempo final de integração numérica da trajetória do veículo espacial T e o primeiro modelo de integral da aceleração, é possível definir PI = 1 T ∫ T 0 |P⃗j|dt (2.27) onde P⃗j é a aceleração devida à cada perturbação analisada e PI é a variação de velocidade devido a cada perturbação considerada no sistema, medida em metros por segundo por período orbital. 2.5 ÓRBITAS CONGELADAS Como apresentado no Capítulo 1, as órbitas congeladas são aquelas em que a variação de e e ω é praticamente nula ao longo do tempo, resultando na manutenção da altitude do veículo espacial. Segundo a literatura, por exemplo Aorpimai e Palmer (2003) e Liu, Baoyin e Ma (2011), essa característica da órbita congelada se deve ao balanceamento das perturbações seculares dos harmônicos zonais pares com as perturbações de longo período dos harmônicos ímpares. No desenvolvimento matemático das equações de órbita congelada é fundamental a introdução das equações planetárias de Lagrange, que são equações variacionais que descrevem o movimento de veículos espaciais em torno de corpos celestes em termos dos elementos orbitais(KUGA; CARRARA; RAO, 2011). As equações consideram o potencial de perturbação (V ) sobre o veículo. Na abordagem empregada por Wu, Jiang e Li (2015) e Liu, Baoyin e Ma (2011), o potencial de perturbação considera apenas o potencial gravitacional (U ) expandido em termos dos harmônicos zonais J2, J3 e J4, ou seja, neste estudo V = U , sendo representado pela eq. (2.15). Dessa forma, de acordo com a literatura (por exemplo, Kuga, Carrara e Rao (2011) e Hoogland (2015)) as equações de Lagrange são expressas por: da dt = 2 na ∂V ∂M (2.28) , 26 de dt = 1− e2 nea2 ∂V ∂M − √ 1− e2 nea2 ∂V ∂ω , (2.29) di dt = cos i na2 √ 1− e2 sin i ∂V ∂ω − 1 na2 √ 1− e2 sin i ∂V ∂Ω , (2.30) dω dt = √ 1− e2 na2e ∂V ∂e − cos i na2 √ 1− e2 sin i ∂V ∂i , (2.31) dΩ dt = 1 na2 √ 1− e2 sin i ∂V ∂i , (2.32) dM dt = n− 1− e2 na2e ∂V ∂e − 2 na ∂V ∂a , (2.33) onde t é o tempo de evolução orbital e n representa a velocidade angular média que, além de ser expressa por 2π/T , como apresentado no início desser capítulo, também é dada por √ µ/a3, onde µ é o parâmetro gravitacional padrão, ou seja, o produto entre a constante de gravidade G e a massa do corpo, para Vênus µ = 0, 32486× 106 km3/s2 (NASA, 2021b). De acordo com Wu, Jiang e Li (2015) e Liu, Baoyin e Ma (2011), a partir das equações planetárias de Lagrange e, considerando os termos de primeira ordem de longo período, obtém-se as taxas de variação média de e e ω dadas respectivamente pelas eqs. (2.34) e (2.35). ˙̄e = 3nJ3R 3 e sin i 4a3 (1− e2)2 ( 5 2 sin2 i− 2 ) cosω (2.34) ˙̄ω = 3nJ2R 2 e 4a2 (1− e2)2 {( 2− 5 2 sin2 i ) × [ 1 + J3Re √ 2 2J2a (1− e2) ( sin2 i− e cos2 i sin i ) sinω e ] + 3J2R 2 e a2 (1− e2)2 D } (2.35) onde D = ( 4 + 7 12 e2 + 2 √ 1− e2 ) − sin2 i ( 103 12 + 3 8 e2 + 11 2 √ 1− e2 ) + sin4 i ( 215 48 − 15 32 e2 + 15 4 √ 1− e2 ) − 35J4 18J2 2 [( 12 7 + 27 14 e2 ) − sin2(i) ( 93 14 + 27 4 e2 ) + sin4(i) ( 21 4 + 81 16 e2 )] . (2.36) Para a ocorrência de uma órbita congelada a variação de e e ω deve ser praticamente nula, sendo assim, a taxa de variação média desses elementos orbitais também deve ser praticamente nula, ou seja, ˙̄e = 0 e ˙̄ω = 0 (WU; JIANG; LI, 2015). Aplicando essa condição na eq. (2.34) observam-se duas possibilidades: 1. 5 2 sin2 i− 2 = 0; 27 2. 3nJ3R3 e sin i 4a3(1−e2)2 cosω = 0. Caso satisfeita a primeira opção, o resultado será um modelo orbital de inclinação crítica que não, necessariamente, representa um modelo de órbita congelada. Portanto, para os casos em que i é 63,435◦ ou 116,565◦, denominam-se casos de inclinação crítica (LIU; BAOYIN; MA, 2011). Sendo assim, apenas a segunda opção resultará em órbitas tipicamente congeladas. Integrando as equações da taxa de variação média, e levando em consideração as condições para a ocorrência de órbitas congeladas, encontra-se a equação para o cálculo de e de órbitas congeladas ao redor de Vênus (WU; JIANG; LI, 2015). e = − [(J3Re)/(2J2a)]sin(i)sin(ω) 1− (3J2R2 eE)/[a 2(5sin2(i)− 4)] (2.37) onde E = ( 6− 169 12 sin2(i) + 395 48 sin4(i) ) − 35J4 18J2 2 ( 12 7 − 93 14 sin2(i) + 21 4 sin4(i) ) . (2.38) A partir da eq. (2.37) e sabendo que a excentricidade de uma órbita é sempre determinada por valores positivos, o argumento do periapsis de órbitas congelada poderá assumir apenas dois valores, dependendo do valor de i empregado: ω = 90◦ ou 270◦. (2.39) Portanto, selecionando valores de a e i convenientes para o cumprimento dos requisitos da missão, calcula-se o valor de e para as órbitas congeladas e define-se o valor conveniente de ω. 2.6 CONDIÇÕES INICIAIS Coffey, Deprit e Deprit (1994) apresentaram a existência de três famílias de órbitas congeladas ao redor de planetas rochosos quando se considera o problema com todos os harmônicos zonais conhecidos. Resultado semelhante foi encontrado para o caso em que se considera apenas os harmônicos zonais J2, J3 e J4 (WU; JIANG; LI, 2015). Para o caso do planeta Vênus, a partir da eq. (2.37) e considerando valores de semieixo maior de 8000 ≤ a ≤ 13000 km e de inclinação de 0 ≤ i ≤ 180◦, encontra-se o gráfico de condições iniciais das órbitas congeladas ao redor de Vênus, como pode ser visto na Figura 2.3 (FILHO, 2020). Assim como apresentado por Coffey, Deprit e Deprit (1994), Wu, Jiang e Li (2015) e Filho (2020), na Figura 2.3 observa-se a ocorrência de uma família de órbitas instáveis em regiões com descontinuidades próximas aos valores de inclinação crítica (representado na figura pelas linhas tracejadas), para este caso, os valores de e aumentam consideravelmente para todos o espectro de a. Além disso, é possível identificar uma outra família de órbitas que tendem aos valores de órbita de inclinação crítica e podem ser classificadas como estáveis, a terceira e última família se refere às órbitas estáveis que contêm a órbita polar. É importante ressaltar que, com exceção das regiões próximas às inclinações críticas ou para os casos em que o veículo se encontra mais próximo ao planeta, as órbitas congeladas apresentam baixa excentricidade. 28 Figura 2.3 – Condições iniciais de órbitas congeladas em torno de Vênus. Fonte: Elaborado pelo próprio autor. 2.7 ESTRATÉGIAS DE CONTROLE LOW-THRUST No desenvolvimento de missões espaciais, o controle da órbita do veículo é essencial para a manu- tenção das características de uma missão, bem como, para a adaptação e correção de órbitas de veículos acometidos por perturbações externas. Duas estratégias clássicas de controle podem ser empregadas para realizar manobras orbitais, a manobra impulsiva e a manobra de empuxo contínuo (CURTIS, 2014). Na primeira, o veículo espacial realiza disparos propulsivos, que podem ser provenientes de combustão química ou elétrica, de modo que há uma mudança instantânea na magnitude e na direção do vetor velocidade do veículo. Já na segunda, o veículo aplica uma força propulsiva de pequena intensidade por um longo período de tempo, sendo essa a de menor consumo energético e, portanto, classificada como low-thrust. Normalmente, o caso contínuo é empregado em manobras de correção orbital, ao passo que, o caso impulsivo é mais utilizado em manobras de transferência orbital, como o famoso caso da manobra bi-impulsiva de Hohmann (CURTIS, 2014). Figura 2.4 – Direção das acelerações radial (S), tangencial (T ) e transversal (W ). Fonte: Elaborado pelo próprio autor. 29 Nesse trabalho, estudam-se as manobras de empuxo contínuo. Para isso, inicialmente, considera- se uma força propulsiva (Fϵ) composta por três projeções de aceleração perturbante constante nas direções radial (S), tangente (T ) e transversal (W ) à órbita, um modelo esquemático desse sistema é apresentado na Figura 2.4, onde v⃗ representa o vetor velocidade do veículo espacial. Além disso, a partir das equações variacionais de Gauss encontra-se as relações das taxas de variação no tempo para os elementos orbitais na presença das acelerações perturbadoras (ALFRIEND et al., 2009). Portanto, a variação secular desses parâmetros a partir das acelerações S, T e W é apresentado nas eqs. (2.40) a (2.45). ȧ = 2 n √ 1− e2 [e sin(θ)S + (1 + e cos(θ))T ] (2.40) ė = √ 1− e2 na [ sin(θ)S + ( cos(θ) + e+ cos(θ) 1 + e cos(θ) ) T ] (2.41) i̇ = r cos(ω + θ) na2 √ 1− e2 W (2.42) Ω̇ = r sin(ω + θ) na2 √ 1− e2 sin(i) W (2.43) ω̇ = √ 1− e2 nae [ − cos(θ)S + ( 1 + r p ) sin(θ)T ] − Ω̇ cos(i), (2.44) Ṁ = n− 1− e2 nae [ − ( cos(θ)− 2e r p ) S + ( 1 + r p ) sin(θ)T ] (2.45) , onde p é o parâmetro focal, também conhecido por semi-latus rectum, representado por p = a(1− e2). Um método muito útil para lidar com as equações variacionais dos elementos orbitais é a teoria do elemento médio. O objetivo principal desse método é remover a dependência do tempo das equa- ções diferenciais originais, permitindo assim uma simplificação considerável da dinâmica resultante (ALFRIEND et al., 2009). Tomando a taxa de variação média em um período orbital de um elemento arbitrário (σ̇med), o método se dá a partir da eq. (2.46). σ̇med = 1 T0 ∫ T0 0 σ̇dt. (2.46) A equação acima depende do tempo, contudo, as equações de movimento em dinâmica orbital são apresentadas com base na anomalia verdadeira (θ). Logo, é necessário realizar a mudança de variável que é apresentada a seguir (WU; JIANG; LI, 2014): dt = √ p3 µ dθ (1 + e cos θ)2 . (2.47) 30 Portanto, obtém-se a equação para o cálculo do elemento médio em função de θ: σ̇med = 1 2π √ (1− e2)3 ∫ 2π 0 σ̇ (1 + e cos θ)2 dθ. (2.48) Wang et al. (2011), apresentou o resultado da aplicação da teoria do elemento médio às taxas de variação no tempo dos elementos orbitais com as componentes de aceleração constante S, T e W . A Tabela 2.1 apresenta o resultado para cada elemento orbital. Tabela 2.1 – Taxa de variação média dos elementos orbitais por componente da aceleração perturba- dora. S T W ȧmed 0 2aT √ p µ 0 ėmed 0 −3 2 eT √ p µ 0 i̇med 0 0 −3 2 e cos(ω) 1−e2 W √ p µ Ω̇med 0 0 −3 2 e sin(ω) (1−e2) sin(i) W √ p µ ω̇med S √ p µ 0 3 2 e sin(ω) cot(i) 1−e2 W √ p µ Ṁmed −3S √ a µ 0 0 Fonte: adaptado de Wang et al. (2011). A partir das taxas de variação média dos elementos orbitais (Tabela 2.1) projeta-se a estratégia de controle de órbitas congeladas a partir de um modelo de aceleração contínua (empuxo contínuo). Como apresentado anteriormente, em uma órbita congelada a variação de e e ω é praticamente nula ao longo do tempo, contudo, sob ação de forças perturbativas, esses parâmetros se modificam. Portanto, os métodos de controle apresentados a seguir são fundamentados para eliminar a rotação de ω. Esse modelo está sendo desenvolvido nos últimos anos, por exemplo, Wilkins e Kapila (2007) eWang et al. (2011) apresentaram as estratégias de controle de órbitas congeladas para o caso terrestre, com equações desenvolvidas para o harmônico zonal J2. Em seguida, Wu, Jiang e Li (2014) expandiu esse modelo para os harmônicos zonais J2, J3 e J4. Portanto, como em Vênus os harmônicos zonais até J4 são de mesma ordem de grandeza, o modelo utilizado neste trabalho é o de Wu, Jiang e Li (2014). Inicialmente, supõe-se que as taxas de variação média do argumento do periapsis são nulas, de modo que: ω̇ = ω̇med + ˙̄ω = 0. (2.49) . Empregando as componentes da aceleração perturbativa e tomando o resultado prático da equação acima (ω̇med = − ˙̄ω), obtém-se quatro estratégias de controle para as órbitas congeladas: o primeiro método emprega o controle por aceleração radial constante, onde as taxas médias de variação do elemento orbital são tomadas em um período orbital; o segundo método compreende um esquema de controle para a direção da aceleração radial, a partir do módulo da força propulsiva permanecer invariável; o terceiro método emprega um esquema de controle a partir da aceleração tangencial, onde 31 o módulo da força propulsiva é constante e a direção é variável; o quarto método é o chamado de low-thrust, que compreende a ação simultânea da aceleração radial e tangencial (WU; JIANG; LI, 2014). A Figura 3.1 apresenta o modelo esquemático das estratégias de controle propostas. Figura 2.5 – Modelo esquemático dos quatro métodos de controle de órbitas congeladas. Fonte: elaborado pelo autor. 2.7.1 Estratégia de controle: método 1 A primeira estratégia de controle foi proposta por Wilkins e Kapila (2007) e consiste no controle a partir da aceleração radial constante, que emprega a média da taxa de variação dos elementos orbitais, originados de acelerações perturbadoras ao longo de um período orbital, eliminando, dessa formas, o crescimento secular dos elementos orbitais. Para os esquema de controle de órbitas congeladas é fundamental o controle de ω, portanto, a partir das taxas de variação dos elementos orbitais devido as acelerações propulsivas contínuas (Tabela 2.1) e da eq. (2.49), obtém-se a relação S √ p µ = ω̇med = − ˙̄ω. (2.50) Logo, a aceleração de controle radial contínua Fϵ será: Fϵ = S = − √ µ p ˙̄ω. (2.51) Substituindo ˙̄ω pela eq. (2.35), obtém-se: Fϵ = − 3nJ2R 2 e 2a4(1− e2)5/2 {( 2− 5 2 sin2 i )[ 1 + J3Re 2J2p ( sin2 i− e2 cos2 i sin i ) sinω e ] + 3J2R 2 e 2p2 D } (2.52) em que D é a expressão apresentada na eq. (2.36). A variação de velocidade em um período orbital indica o consumo de combustível do veículo espacial para concretizar a manobra. Normalmente, denomina-se a variação de velocidade como 32 velocidade característica (WU; JIANG; LI, 2014) e, portanto, a velocidade característica de um período orbital para o método 1 é ∆vT1 = 2πS n = 2πFϵ n . (2.53) A equação anterior pode ser reescrita como: ∆vT1 = − 3π √ µJ2R 2 e p5/2 {( 2− 5 2 sin2 i )[ 1 + J3Re 2J2p ( sin2 i− e2 cos2 i sin i ) sinω e ] + 3J2R 2 e 2p2 D } . (2.54) 2.7.2 Estratégia de controle: método 2 A segunda estratégia de controle considera apenas a aplicação de uma aceleração radial com módulo constante, igual a S. Dessa forma, a partir das eqs. (2.48) e (2.44), a taxa de variação média de ω será: ω̇med = 1 2π (1− e2) 2 nae ∫ 2π 0 −S cos θ (1 + e cos θ)2 dθ. (2.55) Neste esquema de controle, para o sentido adotado e de acordo com a eq. (2.55), Fϵ terá direção radial se 0 ≤ θ ≤ π/2 ou 3π/2 ≤ θ ≤ 2π/2 e, por consequência, sentido anti-radial se π/2 ≤ θ ≤ 3π/2 (WU; JIANG; LI, 2014). Ademais, o sinal do integrando na eq. (2.55) permanecerá invariável em um período orbital, ao passo que a taxa de de variação média de ω aumentará (WANG et al., 2011). Sendo assim, o resultado será ω̇med = 1 2πe √ p µ ( 4 √ 1− e2 − 8e tan−1 √ 1− e 1 + e + 2πe ) S = CSS, (2.56) onde o termo CS é o parâmetro associado a parcela S da aceleração, dado por CS = 1 2πe √ p µ ( 4 √ 1− e2 − 8e tan−1 √ 1− e 1 + e + 2πe ) . (2.57) Portanto, repetindo o processo empregado no método 1 (Fϵ = S) e considerando a relação expressa na eq. (2.49) Fϵ = − 1 CS 3nJ2R 2 e 2p2 {( 2− 5 2 sin2 i )[ 1 + J3Re 2J2p ( sin2 i− e2 cos2 i sin i ) sinω e ] + 3J2R 2 e 2p2 D } . (2.58) Finalmente, a velocidade característica em um período orbital para o método 2 é ∆vT2 = 2πS n . (2.59) 33 Podendo ser reescrita como ∆vT2 = − 1 CS 3πJ2R 2 e p2 {( 2− 5 2 sin2 i )[ 1 + J3Re 2J2p ( sin2 i− e2 cos2 i sin i ) sinω e ] + 3J2R 2 e 2p2 D } . (2.60) 2.7.3 Estratégia de controle: método 3 A terceira estratégia de controle é proposta a partir da componente tangente da aceleração (T ), ou seja, considera esta aceleração com módulo constante e direção variável ao longo do plano de movimento do veículo. Dessa forma, a partir das eqs. (2.48) e (2.44), a taxa de variação média de ω será: ω̇med = 1 2π (1− e2)2 nae ∫ 2π 0 sin(θ)(2 + e cos θ) (1 + e cos θ)3 Tdθ. (2.61) Neste esquema de controle, para o sentido adotado e de acordo com a eq. (2.61), Fϵ terá direção radial se 0 ≤ θ ≤ π e, por consequência, sentido anti-radial se π ≤ θ ≤ 2π (WU; JIANG; LI, 2014). Sendo assim, o resultado será ω̇med = 2 πe √ a µ (2− e2)T = CTT, (2.62) onde o termo CT é o parâmetro associado a parcela T da aceleração, dado por CT = 2 πe √ a µ (2− e2). (2.63) Logo, a partir da eq. (2.49) e considerando agora a aceleração de controle igual à aceleração transversal (Fϵ = T ), obtém-se: Fϵ = − 1 CT 3nJ2R 2 e 2p2 {( 2− 5 2 sin2 i )[ 1 + J3Re 2J2p ( sin2 i− e2 cos2 i sin i ) sinω e ] + 3J2R 2 e 2p2 D } . (2.64) Finalmente, a velocidade característica em um período orbital para o método 3 é ∆vT3 = 2πT n . (2.65) Podendo ser reescrita como: ∆vT3 = −3πJ2R 2 e 2CTp2 {( 2− 5 2 sin2 i )[ 1 + J3Re 2J2p ( sin2 i− e2 cos2 i sin i ) sinω e ] + 3J2R 2 e 2p2 D } . (2.66) 2.7.4 Estratégia de controle: método 4 O quarto esquema de controle contínuo para órbitas congeladas é o de baixo empuxo (low-thrust), que emprega as acelerações radial e transversal, com módulos fixos e direções variáveis. Nesse sentido, 34 tem-se que a variação média de ω é dada pela soma das parcelas radial (S) e tangente (T ) de modo que, satisfazendo a eq. (2.49), obtém-se ω̇med = CSS + CTT = − ˙̄ω. (2.67) Sendo assim, o módulo da aceleração de controle para esse caso é: Fϵ = √ S2 + T 2. (2.68) A partir das relações obtidas anteriormente, a aceleração mínima de controle necessária é obtida quando considera-se os termos CS , CT e ˙̄ω constantes, o que é uma boa aproximação, já que, trata-se de parâmetros com pequena variação em módulo (WANG et al., 2011). Com isso, aceleração de controle mínima (Fϵ,min) será Fϵ,min = ˙̄ω√ C2 S + C2 T . (2.69) Logo, a velocidade característica em um período orbital para o método 4 é ∆vT4 = 2πFϵ,min n . (2.70) Substituindo os termos de empuxo na equação anterior e a taxa de variação média de ω, obtém-se: ∆vT4 = − 3πJ2R 2 e p2 √ C2 S + C2 T {( 2− 5 2 sin2 i )[ 1 + J3Re 2J2p ( sin2 i− e2 cos2 i sin i ) sinω e ] + 3J2R 2 e 2p2 D } . (2.71) 35 3 METODOLOGIA Nesta seção é apresentada a metodologia de análise do trabalho. Inicialmente, é introduzido o pacote computacional para a integração numérica da equação de movimento para obter a evolução orbital do veículo espacial no tempo. Em seguida, são mostrados o grau e a ordem dos harmônicos esféricos considerados no estudo e as regiões de orbitais analisadas. Por fim, apresenta-se a implementação numérica das equações matemáticas do método de controle e a validação numérica do método. 3.1 INTEGRAÇÃO NUMÉRICA DA EQUAÇÃO DE MOVIMENTO As integrações numéricas das equações de movimento do veículo espacial e do Sol, dadas pelas eqs. (2.9) e (2.10), foram realizadas com o pacote numérico de Sanchez, Prado e Yokoyama (2014). Trata-se de um programa escrito na linguagem FORTRAN, compilado e executado em ambiente Linux e que emprega um integrador do tipo Gauss-Radau para a solução numérica das equações diferenciais. Mais informações sobre o pacote numérico podem ser vistas em Pina (2019) e Filho (2020). Originalmente o pacote numérico foi desenvolvido para o caso de um veículo espacial ao redor da Terra. O pacote numérico teve que ser adaptado para o estudo apresentado aqui. Além de fornecer como resultado as variáveis referentes a evolução orbital do veículo espacial no tempo, informações sobre as parcelas de perturbação sofridas pelo veículo espacial devida a cada perturbação considerada no problema também são obtidas. O cálculo é feito utlizando a eq. (2.26) do método da integral da aceleração definida na seção 2.4. Portanto, para esse estudo, adaptou-se o pacote numérico para o caso de Vênus e a validação das alterações foram parcialmente apresentadas por Pina (2019) e Filho (2020). Outras adaptações necessárias foram realizadas durante o desenvolvimento deste trabalho. 3.2 GRAU E ORDEM DOS COEFICIENTES DOS HARMÔNICOS ESFÉRICOS Na integração numérica da equação de movimento, quanto maior a ordem e o grau dos coeficientes dos harmônicos esféricos, maior será o tempo de máquina necessário para obter a solução numérica do problema. Nesse sentido, Filho (2020) mostrou que para o estudo de órbitas congeladas ao redor de Vênus, a consideração de coeficientes dos harmônicos esféricos de grau e ordem 15 são suficientemente adequados para este tipo de estudo, uma vez que, resultam em soluções cujas ordens de grandeza são muito próximas às soluções com ordem e grau superiores. É importante mencionar que ao conssiderar coeficientes harmônicos esféricos de grau e ordem 15, tem-se uma redução significativa do tempo necessário para a execução das integrações numéricas. Outro ponto que deve ser lembrado, o método de controle estudado aqui leva em consideração a expansão matemática de termos associados aos harmônicos zonais J2, J3 e J4, uma vez que esses termos são os mais significativos do potencial gravitacional de Vênus. Sendo assim, neste trabalho a evolução dos elementos orbitais no tempo e o resultado da integral da aceleração foram obtidos a partir da integração numérica da equação de movimento em duas situações, a primeira emprega harmônicos esféricos de grau e ordem 15 e, a segunda, utiliza apenas os harmônicos esféricos zonais 36 J2, J3 e J4. Isso foi realizado com o objetivo de comparar os resultados com as diferentes quantidades de coeficientes de harmônicos esféricos: somente considerando os mais significativos e considerando até ordem e grau 15. 3.3 PARÂMETROS DE VÊNUS E MARTE Os parâmetros essenciais de Vênus e Marte para este estudo são: Re, os coeficientes dos harmônicos zonais (J2, J3 e J4) e o parâmetro gravitacional (µ). Além disso, outro dado importante é o período orbital (PO), ou seja, seu período de translação ao redor do Sol. A Tabela 3.1 apresenta a síntese dos valores para os parâmetros empregados no estudo de acordo com NASA (2021b) e Anderson, Macdonald e Yen (2014). Alguns dos parâmetros de Marte foram obtidos em Mars Fact Sheet1 Um outro conjunto de dados importantes para o pacote que realiza a integração numérica da equação de movimento, são as efemérides iniciais de Vênus. Isso foi obtido a partir do HORIZONS Web- Interface2, com época de 20 de Janeiro de 2022. Tabela 3.1 – Elementos orbitais e parâmetros dos planetas Vênus e Marte. Parâmetro Valor Vênus Valor Marte Unidade Re 6051, 8 3396, 2 km µ 0, 32486× 106 0, 042828× 106 km3/s2 PO 224, 7 686, 98 dias a 108, 210× 106 227, 956× 106 km e 0, 00677323 0, 09341233 - i 3, 39471 1, 85061 grau Ω 76, 68069 49, 57854 grau ω 131, 53298 336, 04084 grau M 181, 97973 355, 45332 grau J2 4, 4580× 10−6 1, 95545× 10−3 - J3 −2, 1082× 10−6 3, 14498× 10−5 - J4 −2, 1471× 10−6 −1, 53774× 10−5 - Fonte: elaborado pelo próprio autor. 3.4 DEFINIÇÃO DA REGIÃO DE ESTUDO A definição da região de estudo é importante para estabelecer os limites da aplicação do método de controle, uma vez que, ambas as metodologias de cálculo das condições iniciais das órbitas congeladas e de controle por empuxo contínuo, consideram as forças perturbativas do potencial gravitacional do planeta expandida em termos dos coeficientes dos harmônicos zonais. Dessa forma, regiões em que a ação perturbativa do Sol é dominante em relação ao potencial gravitacional, não serão representativas para o método empregado. Nesse sentido, fez-se a integração numérica da equação de movimento do veículo para um intervalo de tempo de 30 períodos orbitais de Vênus para uma órbita na região polar do planeta com condições iniciais: i = 85, 5◦, ω = 90◦, Ω = 0◦ e M = 0◦. Com isso, estabeleceu-se seis casos de análise, onde 1 https://nssdc.gsfc.nasa.gov/planetary/factsheet/marsfact.html 2 https://ssd.jpl.nasa.gov/horizons/app.html/ 37 o valor do semieixo maior (a) foi escolhido de forma crescente e, para cada situação, calculou-se o valor inicial da excentricidade a partir da eq. (2.37). Em seguida, a partir dos resultados obtidos com a integral da aceleração, comparou-se os valores obtidos para a variação de velocidade, em m/s, sofrida pelo corpo em razão da ação perturbativa do potencial gravitacional do planeta (PIg) e da presença do terceiro corpo (PISol). A Tabela 3.2 apresenta os valores de a e e para cada caso de estudo, além dos valores resultantes da integral da aceleração. Os resultados expressos na Tabela 3.2 permitiram observar que o aumento do valor de a reflete na elevação dos efeitos perturbativos devidos à presença do Sol, ao passo que menores valores de a resultam na maior perturbação devido ao potencial gravitacional do planeta. Os resultados também mostraram que para a = 15551, 80km e e = 0, 1055 parece existir uma equivalência em módulo das variações das velocidades devidas às perturbações consideradas (ver órbita 4 da Tabela 3.2). Diante desses resultados, foi possível estabelecer os limites da região de estudo, visto que o método de controle demanda uma dominância da perturbação devido ao potencial gravitacional de Vênus, é possível inferir que o método é melhor empregado em regiões inferiores a a = 15551, 80km, ou seja, regiões mais próximas à superfície do planeta. A fim de ampliar a análise, o mesmo procedimento foi repetido para órbitas de i = 45◦ e i = 110◦, resultando em comportamento semelhante ao obtido para a órbita polar nas mesmas faixas de valores de a. Tabela 3.2 – Variação de velocidade gerada no veículo devido a ação perturbativa do potencial gravita- cional de Vênus (PIg) e da presença do Sol (PISol). Órbita a (km) e PIg (m/s) PISol (m/s) 1 8751,80 0,2776 1033,7800 2,3579 2 12751,80 0,1389 100,5000 33,6300 3 14751,80 0,1131 51,5250 38,8490 4 15551,80 0,1055 40,7870 40,9350 5 16751,80 0,0960 29,5950 44,1030 6 19751,80 0,0786 14,7020 51,9800 Fonte: elaborado pelo próprio autor. 3.5 IMPLEMENTAÇÃO NUMÉRICA DO MÉTODO DE CONTROLE Os métodos de controle apresentados na seção 2.7 foram implementados numericamente no software MATLAB. Como visto nas equações, os resultados dependem das variáveis de entrada (elementos orbitais), e das características do planeta estudado, especificamente, o raio equatorial (Re), os coeficientes dos harmônicos zonais (J2, J3 e J4) e o parâmetro gravitacional (µ). Dessa forma, obtém-se para cada método de controle o valor da aceleração de controle (Fϵ) e a velocidade característica para cada um dos métodos. Com o objetivo de validar a implementação numérica, reproduziram-se os resultados obtidos por Wu, Jiang e Li (2014). Naquele trabalho, os autores estudaram a aplicação dos métodos de controle de órbitas congeladas na Terra e em Marte para a condição inicial dos elementos orbitais: a = 10.000km, e = 0, 2 , i = 30◦ e ω = Ω =M = 0◦. Mais detalhes podem ser encontrados na referência citada. Nas Tabelas 3.3 e 3.4 são mostrados os resultados para o caso da Terra e de Marte, respectivamente. 38 Tabela 3.3 – Comparação entre o resultado obtido por Wu, Jiang e Li (2014) e pelo programa imple- mentado para o caso da Terra considerando J2, J3 e J4. Resultado Wu, Jiang e Li (2014) Resultado obtido Método Fϵ (m/s2) ∆vT (m/s) Fϵ (m/s2) ∆vT (m/s) 1 4, 0163× 10−3 39, 9190 4, 0016× 10−3 39, 9144 2 1, 2310× 10−3 12, 2546 1, 2324× 10−3 12, 2928 3 6, 3076× 10−4 6, 2773 6, 2844× 10−4 6, 2684 4 5, 7326× 10−4 5, 7051 5, 5985× 10−4 5, 5843 Fonte: Elaborado pelo próprio autor. Tabela 3.4 – Comparação entre o resultado obtido por Wu, Jiang e Li (2014) e pelo programa imple- mentado para o caso de Marte considerando J2, J3 e J4. Resultado Wu, Jiang e Li (2014) Resultado obtido Método Fϵ (m/s2) ∆vT (m/s) Fϵ (m/s2) ∆vT (m/s) 1 2, 1997× 10−4 6, 6733 2, 2019× 10−4 6, 6842 2 6, 7441× 10−5 2, 0473 6, 7813× 10−5 2, 0586 3 3, 4546× 10−5 1, 0487 3, 4580× 10−5 1, 0497 4 3, 0747× 10−5 0, 9334 3, 0806× 10−5 0, 9352 Fonte: Elaborado pelo próprio autor. Pode-se notar que, os resultados obtidos têm a mesma ordem de grandeza e valores muito próximos àqueles da literatura considerada. Além disso, é possível observar que os valores da aceleração de controle são menores para o método 4, o que demonstra a economia energética resultante deste método de baixo empuxo. Outra observação interessante é que os valores de controle em Marte são menores em relação àqueles da Terra. Esse fato permite indicar que para as mesmas condições iniciais, a perturbação da órbita sofrida pelo veículo espacial devido ao potencial gravitacional de Marte é menor em relação àquele da Terra. De forma geral, como foi apresentado nos capítulos anteriores, as órbitas congeladas são caracteri- zadas por manter praticamente constante ao longo do tempo ω e e, portanto, o método de controle deve ser aplicado às órbitas que ultrapassem um limite preestabelecido para a variação desses elementos. Nesse sentido, a metodologia de aplicação do método é dividida em três etapas: na primeira parte calcula-se os parâmetros da órbita congelada, conforme método apresentado na seção 2.5; a segunda etapa considera os parâmetros calculados anteriormente como condições iniciais para a integração numérica da equação de movimento; na última etapa avalia-se a variação dos elementos orbitais e, caso eles tenham sido perturbados a ponto de extrapolar o limite estabelecido, aplica-se o método de controle para a restauração das características congeladas da órbita. Resumindo, a descrição dos passos necessários para cumprir as etapas apresentadas para aplicação do método são: 1ª ETAPA: • Definição de a, i, Ω e M para a órbita inicial do veículo; • A partir dos elementos orbitais escolhidos, calcula-se e e ω, respectivamente, pelas eqs. (2.37) e (2.39). De modo a obter a condição incial da órbita congelada. 2ª ETAPA: 39 • Integração numérica da equação de movimento do veículo (2.11), para a condição inicial de óbita congelada calculada na 1ª ETAPA. A integração se dá a partir do pacote numérico de Sanchez, Prado e Yokoyama (2014); • Da integração numérica obtém-se: a evolução dos elementos orbitais em função do tempo e o resultado da integral da aceleração para cada termo perturbativo. 3ª ETAPA: • Define-se um limite de variação máximo para e e ω; • Verificação se na evolução de e e ω no tempo ultrapassaram os limites estabelecidos; • Caso tenha ultrapassado, aplica-se a estratégia de controle, que pode ser os Métodos 1, 2, 3 ou 4, conforme apresentado na seção 2.7; • Da aplicação do método obtem-se: a velocidade característica e a aceleração de controle; • Caso os limites não sejam ultrapassados, a órbita do veículo não precisa ser corrigida. Para facilitar o entendimento do funcionamento do método um fluxograma é apresentado na Figura 3.1. Figura 3.1 – Fluxograma do método para o cálculo da aceleração de controle. Fonte: Elaborado pelo autor. 40 4 RESULTADOS E DISCUSSÕES Neste capítulo são mostrados e discutidos os resultados obtidos a partir da metodologia estabelecida no capítulo anterior. Inicialmente, é mostrada uma análise do efeito dos parâmetros orbitais na obtenção da aceleração de controle e a relação entre esses parâmetros. Em seguida são exploradas as órbitas de inclinação abaixo e acima da condição crítica, bem como, na própria região de inclinação crítica. 4.1 EFEITO DOS ELEMENTOS ORBITAIS NOS MÉTODOS DE CONTROLE Para avaliar a consequência de cada parâmetro orbital no método de controle, tomou-se as condições iniciais (a = 10.000 km, i = 30◦, e = 0, 2 , ω = 0◦, Ω = 0◦ e M = 0◦) vistas na seção 3.5 para elaborar os gráficos de variação da aceleração de controle em função de a, e , i e ω. Vale a pena ressaltar que o período orbital de um veículo espacial ao redor de Vênus com àquelas condições iniciais é ≈ 8,43 horas, enquanto que ao redor de Marte é ≈ 3,06 horas. Com o objetivo de comparar do resultados para o caso Vênus com o caso de um outro planeta do Sistema Solar, repetiu-se o procedimento para o caso de Marte, de modo que foram utilizados os resultados e parâmetros do planeta mostrados na Tabela 3.1. Os resultados obtidos são mostrados na Figura 4.1. Em uma primeira análise, percebe-se que as curvas obtidas para os métodos de controle em Marte e em Vênus apresentam comportamento semelhantes, com ordens de grandezas distintas devido aos parâmetros dos sistemas. Dessa forma, em ambos os planetas, o Método 1 é o que provocou o maior valor em módulo da aceleração de controle necessária para realizar a correção orbital do veículo. Logo, esse não seria um bom método, uma vez que o consumo energético da manobra seria maior em relação aos outros 3 métodos. Nesse sentido, dentre os Métodos 2, 3 e 4, o que resultou nos menores valores em módulo da aceleração de controle foi o Método 4, confirmando, portanto, a característica low-thrust dessa estratégia. Desse modo, nos casos apresentados adiante o foco da análise será o Método 4. Em seguida, com relação a magnitude da resposta, nota-se que em Vênus, a aceleração necessária para a correção resultou em valores bem menores do que para Marte. Lembrando que a metodologia de controle considera apenas o efeito da não esfericidade do planeta e, a região de Vênus em que foi realizada a análise se encontra dentro do limite estabelecido na seção 3.4, a baixa magnitude da resposta para esse planeta pode beneficiar uma missão espacial através do aumento de sua vida útil. A partir da Figura 4.1 também é possível avaliar qualitativamente o efeito dos elementos orbitais para cada um dos métodos na obtenção da aceleração de controle. Devido ao comportamento seme- lhante das curvas as observações sobre o efeito desses parâmetros se estende para ambos os planetas comparados. Nesse contexto, observa-se que a aceleração resultante apresenta maior valor em módulo para os menores valores de a, ou seja, para regiões mais próximas à superfície do corpo orbitado a aceleração necessária para realizar a correção será maior. Quanto à excentricidade, nota-se que a aceleração é maior para as órbitas de elevado e, pois quanto maior o valor de e menor a distância de periapsis do veículo, mais próximo o veículo passará do planeta e mais perturbada a órbita do veículo naquele ponto. Podemos notar também que para 0, 0 < e < 0, 4, o valor dessa aceleração é praticamente constante em todas as estratégias de controle. Para inclinações inferiores a 60◦ e 41 superiores a 120◦, a aceleração de controle sofre grandes variações, enquanto que para as regiões de i em torno da condição polar (i = 90◦), o valor da aceleração tende a um valor constante para os métodos 2, 3 e 4. Por fim, quando analisamos o gráfico da aceleração de controle em função de ω, pode-se ver que para 0◦ < ω < 180◦ a aceleração de controle tem valores constantes para o caso de Marte enquanto que para o caso de Vênus os métodos 3 e 4 tem o mesmo comportamento, sendo praticamento constante para 60◦ < ω < 120◦. Já para os métodos 1 e 2, a aceleração de controle no intervalo considerado, tem um comportamento senoidal no qual o seu valor máximo é em torno de 90◦. Figura 4.1 – Variação da aceleração de controle, para os quatro métodos de análise, em função de a, e e i em (a) Marte ( resultados na coluna à esquerda) e (b) Vênus (resultados na coluna à direita). As condições iniciais da órbita são: a = 10.000 km, i = 30◦, e = 0, 2 , ω = 0◦, Ω = 0◦ e M = 0◦. Continua na próxima página. 42 (a) Marte (b) Vênus Fonte: Elaborado pelo próprio autor. Por fim, para compreender melhor a relação entre os parâmetros orbitais, avaliou-se a inter-relação desses elementos para Vênus considerando apenas o método low-thrust (Método 4). Para tanto fez-se três casos de análise: no primeiro caso fixou-se e = 0, 01 (órbita praticamente circular) e obteve-se a superfície da aceleração de controle para 8000 ≤ a ≤ 16000km e 5 ≤ i ≤ 175◦; no segundo caso fixou-se a = 10.000km e variou-se 0, 01 ≤ e ≤ 0, 8 e 5 ≤ i ≤ 175◦; no último caso fixou-se i = 45◦ e variou-se 0, 01 ≤ e ≤ 0, 8 e 8000 ≤ a ≤ 16000km. As superfícies resultantes são apresentadas na Figura 4.2. Neste caso é importante ressaltar que os motivos da escolha dos dados da análise foram uma adaptação a características de órbitas congeladas e para melhor atender aos limites do método de controle. Nesse sentido, a escolha do valor fixo de e no primeiro caso, foi porque as órbitas congeladas típicas tendem a apresentar excentricidades quase circular. Além disso, em todos os casos o valor de a foi escolhido para contemplar os valores estabelecidos na seção 3.4. Já no terceiro caso, o valor fixo de i foi escolhido com base no trabalho de Filho (2020), onde mostrou-se que órbitas congeladas nessa região de valores de inclinação tendem a condições estáveis dos elementos orbitais. Por último, no intervalo de i desconsiderou-se a região em que sin(i) seria próximo a zero, pois de acordo com as eqs. (2.69) e (2.35) isso resultaria em uma singularidade. 43 Figura 4.2 – Superfície da aceleração de controle em m/s2 para o método low-thrust e função de (a) a e i, (b) e e i, (c) e e a. (a) (b) (c) Fonte: Elaborado pelo próprio autor. Para os gráficos obtidos, a Figura 4.2(a) representa o primeiro caso de análise, onde nota-se a presença de uma grande variação da aceleração de controle em regiões mais próximas à superfície 44 de Vênus, mais especificamente para valores de a entre 8.000 e 11.000km, esse fato ocorre em todo o espectro de inclinação. Já o segundo caso é apresentado na Figura 4.2(b), nela vê-se que ocorre uma elevada influência da excentricidade quando e atinge valores mais próximos a 1. Tal efeito é acentuado nas regiões limites de inclinação. Por último, o terceiro caso é mostrado na Figura 4.2(c), nesta situação os fatores que resultam nos maiores valores para a aceleração de controle são as regiões próximas à superfície de Vênus e de alta excentricidade. Em suma, das observações acima pode-se dizer que para o método low-thrust, visando a redução do gasto energético, as regiões que seriam mais adequadas para a aplicação da correção são as de excentricidade no máximo 0,4, inclinações mais distantes de 0◦ e 180◦ e semieixo maior mais próximo da região onde ocorre o equilíbrio das perturbações devido a não esfericidade de Vênus e à presença do Sol. 4.2 AVALIAÇÃO DAS ÓRBITAS CONGELADAS NAS REGIÕES DE INTERESSE AO REDOR DE VÊNUS Anteriormente foi apresentado o comportamento do método de controle para uma órbita arbitrária, com o objetivo de discutir o comportamento geral dessas situações. Nesta seção, o intuito é investigar órbitas cujas condições iniciais são de uma órbita congelada. Para tanto, respeitando os limites de análise estabelecidos na metodologia, selecionou-se 20 órbitas com valores de semieixo maior 9.000 < a < 15.000 km (ver seção 2.6, Figura 2.3) para órbitas com: • Inclinação abaixo da crítica: são as órbitas que tende à inclinação crítica. Os valores de inclinação dessa região, no primeiro quadrante, variam entre 0 e 60◦. Desse modo, selecionou-se dois casos de estudo, sendo i = 35◦ e i = 45◦. • Região de inclinação crítica: são as órbitas congeladas que tendem a uma maior instabilidade. Considerando o primeiro quadrante, escolheu-se o caso em que i = 62◦. • Região de inclinação polar: é a família de órbitas que compreende o entorno de i = 90◦. Para essa região selecionou-se os casos de i = 75◦ e i = 85◦. Para cada um dos casos tomaram-se quatro valores de a e com as eqs. (2.37) e (2.39) calcularam-se os elementos orbitais da órbita inicial do veículo. O resultado de e e ω para cada valor de i e a é apresentado na Tabela 4.1. Os demais elementos foram considerados nulos (Ω =M = 0◦). Portanto, com as considerações acima, na Tabela 4.1 as órbitas de 1 a 8 são casos compreendidos na região de abaixo da inclinação crítica. Já as órbita de 9 a 12 referem-se à região de inclinação critica. As demais pertencem à região de inclinação polar. A partir das condições inciais do veículo espacial, fez-se a integração numérica da equação de movimento do veículo por um tempo de 30 períodos orbitais de Vênus considerando, inicialmente, os coeficientes dos harmônicos esféricos de grau e ordem 15 e, em seguida, somente os coeficientes dos harmônicos zonais J2, J3 e J4. Nesse contexto, as seções a seguir apresentam a evolução orbital de ω e e, o módulo da aceleração de controle e da variação de velocidade para cada órbita. Para isso, empregou-se a estratégia low-thrust. 45 Tabela 4.1 – Condições iniciais das órbitas congeladas selecionadas para o estudo. Órbita i (grau) a (km) e ω (grau) 1 35 9000 0,0867 90 2 11000 0,0721 90 3 13000 0,0616 90 4 15000 0,0537 90 5 45 9000 0,1498 90 6 11000 0,1104 90 7 13000 0,0884 90 8 15000 0,0741 90 9 62 9000 0,0583 270 10 11000 0,0897 270 11 13000 0,1535 270 12 15000 0,3714 270 13 75 9000 0,184 90 14 11000 0,1413 90 15 13000 0,1155 90 16 15000 0,098 90 17 85 9000 0,258 90 18 11000 0,1748 90 19 13000 0,1346 90 20 15000 0,1104 90 Fonte: Elaborado pelo próprio autor. 4.2.1 Região de inclinação abaixo da crítica Nesta seção são apresentados os resultados obtidos para a região de inclinação abaixo da crítica, ou seja, as órbitas de 1 a 8 da Tabela 4.1. Para realizar a aplicação do controle orbital adotou-se duas estratégias: a primeira é considerar um limite máximo na variação de ω, já a segunda o limite é imposto sob a variação de e. Dessa forma, se a evolução orbital do elemento analisado for superior a essa variação, aplica-se a correção orbital. A seguir são discriminados os limites dessas variações na região de estudo. 4.2.1.1 Limitando ω A evolução orbital de ω em 30 períodos orbitais de Vênus para as órbitas de i = 35◦ é apresentada na Figura 4.3 e para i = 45◦ na Figura 4.4. Nessas Figuras, a coluna da esquerda representa a evolução do elemento orbital para o potencial gravitacional de Vênus até coeficientes dos harmônicos esféricos de grau e ordem 15, ao passo que na coluna da direita é a evolução orbital considerando somente os coeficientes dos harmônicos zonais J2, J3 e J4. Inicialmente, nota-se que para i = 35◦ o comportamento das curvas apresentam semelhanças quanto às perturbações, ou seja, na Figura 4.3(a) e (b) ocorreu uma menor perturbação de ω com a elevação do valor de a. Entretanto, as curvas para o caso de coeficientes dos harmônicos esféricos até grau e ordem 15 apresentaram perturbações de curto período não presentes no caso de harmônicos zonais. Discutir a origem dessas perturbações e demais características não são objetivos deste trabalho, portanto, tal comportamento não será aprofundado. 46 Por sua vez, para i = 45◦ o comportamento das curvas foi semelhante nas duas situações. Como no caso de 35◦, percebe-se que o aumento de a reduz o efeito perturbativo em ω. Figura 4.3 – Evolução orbital de ω em 30 períodos orbitais de Vênus e i = 35◦ considerando (a) os coeficientes dos harmônicos esféricos até ordem e grau 15 e (b) os coeficientes dos harmônicos zonais J2, J3 e J4. (a) (b) Fonte: elaborado pelo próprio autor. Figura 4.4 – Evolução orbital de ω em 30 períodos orbitais de Vênus e i = 45◦ para (a) os coeficientes dos harmônicos esféricos até ordem e grau 15 e (b) os coeficientes dos harmônicos zonais J2, J3 e J4. (a) (b) Fonte: Elaborado pelo próprio autor. A partir dos resultados obtidos das integrações numéricas da equação de movimento do veículo avaliou-se a aplicação do método de controle e o tempo em períodos orbitais de Vênus (PO) que é necessário aplicar a correção da órbita. Para tanto, visando a redução do consumo energético, estabeleceu-se um limite de variação de pequenos ângulos para ω, de modo que escolheu-se como limite ∆ω = 3◦ e ∆ω = 5◦. Portanto, aplicando o controle quando a órbita atingia essas restrições de variação obteve-se os resultados expressos na Tabela 4.2. Uma observação importante é que, na tabela 47 a seguir e das próximas seções, os traços indicam as órbitas que não atingiram o limite de variação imposto. Tabela 4.2 – Aplicação do método low-thrust para as órbitas da região de inclinação abaixo da crítica, considerando ∆ω = 3◦ e ∆ω = 5◦. ∆ω = 3◦ ∆ω = 5◦ Órbita Coeficientes dos Harmônicos esféricos Tempo (PO) ∆vT4 (m/s) Fϵ,min (m/s2) Tempo (PO) ∆vT4 (m/s) Fϵ,min (m/s2) 1 Até 15 29,94 4, 689976× 10−4 4, 982756× 10−8 - - - J2, J3 e J4 10,25 2, 277138× 10−4 2, 162191× 10−8 21,00 6, 790086× 10−4 7, 213975× 10−8 2 Até 15 - - - - - - J2, J3 e J4 - - - - - - 3 Até 15 - - - - - - J2, J3 e J4 - - - - - - 4 Até 15 - - - - - - J2, J3 e J4 - - - - - - 5 Até 15 3,84 2, 009248× 10−4 2, 134711× 10−8 6,06 1, 652500× 10−4 1, 755659× 10−8 J2, J3 e J4 5,29 1, 213943× 10−4 1, 289677× 10−8 9,11 3, 539702× 10−5 3, 760613× 10−9 6 Até 15 12,47 4, 312733× 10−5 3, 390944× 10−9 21,11 1, 412856× 10−5 1, 110886× 10−9 J2, J3 e J4 13,15 1, 958010× 10−5 1, 539524× 10−9 22,09 3, 297344× 10−5 2, 592613× 10−9 7 Até 15 21,04 9, 328959× 10−6 5, 709313× 10−10 - - - J2, J3 e J4 19,09 1, 008762× 10−5 6, 173598× 10−10 - - - 8 Até 15 22,04 8, 624988× 10−6 4, 258820× 10−10 - - - J2, J3 e J4 20,11 1, 131726× 10−5 5, 588189× 10−10 - - - Fonte: Elaborado pelo próprio autor. Da Tabela 4.2 nota-se que os casos de i = 35◦ (órbitas 1 a 4), em geral, não sofreram variação que ultrapassasse o limite de controle, portanto, não foi necessário aplicar nenhuma correção. Entretanto, há uma exceção para a órbita 1, esse caso resultou na necessidade de aplicação de controle para um tempo de 29,94 períodos orbitais de Vênus, quando considera-se os coeficientes dos harmônicos esféricos até grau e ordem 15. Ao passo que, para apenas os coeficientes dos harmônicos zonais J2, J3 e J4, a aplicação do método foi necessária para os dois limites. Ademais, para as órbitas de i = 45◦ (5 a 8) o método foi aplicado em todos os casos quando ∆ω = 3◦ e necessário apenas para os casos mais próximos à superfície (órbitas 5 e 6) quando o limite de variação era de 5◦. Para esses casos não houve uma grande diferença no tempo de aplicação da correção. Diante desses dados, pode-se verificar que o valor da menor inclinação aparenta promover uma maior estabilidade da órbita do veículo e, portanto, reduz-se a necessidade da aplicação do método de correção. Além disso, para esses casos em que a órbita é originalmente congelada, e que somente a perutrabação da não esfericidade do planeta é levada em consideração, a variação de velocidade necessária para retomar condições congeladas é menor. 4.2.1.2 Limitando e As mesmas órbitas analisadas para a variação de ω são novamente estudadas, só que agora o elemento orbital considerado para a aplicação do método de controle é a excentricidade. Nesse sentido, a Figura 4.5 apresenta a evolução orbital de e para i = 35◦ e a Figura 4.6 para i = 45◦. Nessas Figuras (a) é a evolução considerando os coeficientes dos harmônicos esféricos até grau e ordem 15 e (b) considerando somente os coeficientes dos harmônicos zonais J2, J3 e J4. 48 Na Figura 4.5 nota-se que e permanece praticamente constante para os casos mais distantes da superfície de Vênus, ao passo que, a maior variação de e ocorre para a órbita de menor a. Tal efeito pode ser visto nos dois casos de potencial gravitacional. Além disso, assim como foi observado na evolução orbital de ω, o caso de e para harmônicos até grau e ordem 15 é marcado pela presença de perturbações de curto período, as quais não serão aprofundadas. Já na Figura 4.6, é possível observar que a evolução orbital de e em ambos os casos é basicamente o mesmo, não apresentando grandes diferenças qualitativas. Figura 4.5 – Evolução orbital de e em 30 períodos orbitais de Vênus e i = 35◦ considerando (a) os coeficientes dos harmônicos esféricos até ordem e grau 15 e (b) coeficientes dos harmônicos zonais J2, J3 e J4. (a) (b) Fonte: Elaborado pelo próprio autor. Figura 4.6 – Evolução orbital de e em 30 períodos orbitais de Vênus e i = 45◦ considerando (a) coeficientes dos harmônicos esféricos até ordem e grau 15 e (b) os coeficientes dos harmônicos zonais J2, J3 e J4. (a) (b) Fonte: Elaborado pelo próprio autor. Com a evolução orbital de e no tempo, avaliou-se a aplicação do controle orbital para dois casos limites de variação de e. Considerando o comportamento regular da excentricidade de órbitas 49 congeladas definiu-se pequenas variações neste parâmetro estabelecendo, portanto, como limites ∆e = 0, 001 e ∆e = 0, 01. A Tabela 4.3 apresenta os dados de aceleração de controle, variação de velocidade e tempo de evolução em que deve ocorrer a manobra, para quando foi necessário aplicar a estratégia low-thrust. Tabela 4.3 – Aplicação do método low-thrust para as órbitas da região de inclinação abaixo da crítica, considerando ∆e = 0, 001 e ∆e = 0, 01. ∆e = 0, 001 ∆e = 0, 01 Órbita Coeficientes dos harmônicos esféricos Tempo (PO) ∆vT4 (m/s) Fϵ,min (m/s2) Tempo (PO) ∆vT4 (m/s) Fϵ,min (m/s2) 1 Até 15 0,14 1, 511228× 10−5 1, 605559× 10−9 - - - J2, J3 e J4 9,72 2, 401376× 10−4 2, 551320× 10−8 - - - 2 Até 15 0,28 3, 063350× 10−5 2, 408610× 10−9 - - - J2, J3 e J4 - - - - - - 3 Até 15 0,42 2, 716436× 10−5 1, 662443× 10−9 - - - J2, J3 e J4 8,29 2, 009047× 10−5 1, 229525× 10−9 - - - 4 Até 15 5,15 2, 005223× 10−5 9, 901302× 10−10 - - - J2, J3 e J4 5,34 1, 864137× 10−5 9, 204612× 10−10 - - - 5 Até 15 0,18 1, 778131× 10−4 1, 889080× 10−8 16,19 6, 418413× 10−4 6, 818728× 10−8 J2, J3 e J4 7,74 4, 998744× 10−5 5, 310861× 10−9 21,49 5, 630641× 10−4 5, 981900× 10−8 6 Até 15 0,42 4, 912737× 10−5 3, 862675× 10−9 - - - J2, J3 e J4 21,78 5, 203003× 10−5 4, 091035× 10−9 - - - 7 Até 15 3,42 2, 345032× 10−5 1, 435141× 10−9 - - - J2, J3 e J4 5,34 2, 248733× 10−5 1, 376207× 10−9 - - - 8 Até 15 2,38 1, 561940× 10−5 7, 712431× 10−10 - - - J2, J3 e J4 4,26 1, 558504× 10−5 7, 695451× 10−10 - - - Fonte: Elaborado pelo próprio autor. A partir da Tabela 4.3, pode-se verificar que para os casos de i = 35◦ (órbitas 1 a 4) houve necessidade de aplicação do método em quase todos os casos para ∆e = 0, 001, a exceção é a órbita 2 para harmônicos zonais. Ademais, em nenhuma das situações em que ∆e = 0, 01 foi necessário corrigir a trajetória. Das órbitas que necessitaram correção em ambas as situações de potencial gravitacional, o tempo para a aplicação da estratégia naquelas em que considerou-se os coeficientes dos harmônicos esféricos até grau e ordem 15 foi bem menor do que para o caso dos coeficientes dos harmônicos zonais J2, J3 e J4. Quanto às órbitas de i = 45◦ (5 a 8), todas necessitaram da aplicação da estratégia de controle quando ∆e = 0, 001 e apenas a órbita 5 atingiu o limite de ∆e = 0, 01. Assim como no caso de i = 35◦, o tempo para a aplicação da manobra foi menor quando levou-se em consideração mais harmônicos do potencial gravitacional. Por fim, os valores de variação de velocidade necessários para a concretização da correção apresenta pequeno módulo, mas ao contrário do que é observado na Tabela 4.2, neste caso a ordem de grandeza foi a mesma na maior parte das situações analisadas. 4.2.2 Região de inclinação crítica A segunda região de estudo é a de inclinação crítica, cujas condições iniciais são apresentadas na Tabela 4.1 pela órbitas de 9 a 12. Será retomado o mesmo procedimento de análise da região anterior. A seguir serão apresentados os resultados para a limitação de ω e e. 50 4.2.2.1 Limitando ω A evolução orbital de ω em 30 períodos orbitais de Vênus para a órbita de i = 60◦ é apresentada na Figura 4.7. A Figura 4.7(a) mostra a evolução para coeficientes dos harmônicos esféricos até grau e ordem 15, nela pode-se observar que a variação de ω da órbita 9 é mais intensa do que as demais, além disso, diferentemente do que foi visto na região de inclinação abaixo da crítica, nesta região de inclinação crítica as órbitas não tendem a uma maior estabilidade com o aumento da distância em relação à superfície do planeta, já que, com exceção da órbita 9, elas estão concentradas em uma mesma faixa de variação de valores. Por sua vez a Figura 4.7(b) apresenta a evolução para os coeficientes dos harmônicos zonais J2, J3 e J4. Neste caso, diferentemente de tudo do que foi observado até agora, as órbitas mais próximas à superfície tenderam a uma maior estabilidade no comportamento do elemento orbital. Ademais, a variação dos dados na própria curva não é tão intensa como para o caso de (a). Figura 4.7 – Evolução orbital de ω em 30 períodos orbitais de Vênus e i = 62◦ para (a) harmônicos até ordem e grau 15 e (b) harmônicos zonais. (a) (b) Fonte: Elaborado pelo próprio autor. A partir da evolução orbital de ω em 30 PO e considerando, novamente, o limite ∆ω = 3◦ e ∆ω = 5◦, na Tabela 4.4 são apresentados os resultados da aplicação do método de controle para as órbitas da região crítica que ultrapassaram o limite estabelecido. Tabela 4.4 – Aplicação do método low-thrust para a órbita da região de inclinação crítica, considerando ∆ω = 3◦ e ∆ω = 5◦. ∆ω = 3◦ ∆ω = 5◦ Órbita Coeficientes dos harmônicos esféricos Tempo (PO) ∆vT4 (m/s) Fϵ,min (m/s2) Tempo (PO) ∆vT4 (m/s) Fϵ,min (m/s2) 9 Até 15 5,36 2, 369547× 10−5 2, 517511× 10−9 8,62 2, 229989× 10−5 2, 369220× 10−9 J2, J3 e J4 10,84 1, 915082× 10−5 2, 034654× 10−9 17,83 3, 193945× 10−5 3, 393358× 10−9 10 Até 15 8,56 1, 268983× 10−5 9, 977891× 10−10 12,91 2, 369012× 10−5 1, 862706× 10−9 J2, J3 e J4 8,96 5, 640435× 10−6 4, 439530× 10−10 14,95 1, 496646× 10−5 1, 176787× 10−9 11 Até 15 6,66 1, 833980× 10−5 1, 122408× 10−9 11,84 1, 354948× 10−7 8, 292250× 10−12 J2, J3 e J4 7,58 1, 893572× 10−5 1, 158891× 10−9 12,07 1, 037997× 10−5 6, 352544× 10−10 12 Até 15 6,58 3, 418728× 10−4 1, 688134× 10−8 10,88 3, 039305× 10−4 1, 500727× 10−8 J2, J3 e J4 6,73 3, 376529× 10−4 1, 667244× 10−8 11,15 3, 180592× 10−4 1, 570534× 10−8 Fonte: Elaborado pelo próprio autor. 51 Na Tabela acima pode-se verificar que todas as órbitas estudadas necessitaram da aplicação do controle orbital para manter condições congeladas. Além disso, nota-se que, em geral, ambos os casos de configuração do potencial gravitacional resultaram em tempo semelhante de aplicação da estratégia. A única exceção é a órbita 9, na qual a evolução para apenas os coeficientes dos harmônicos zonais J2, J3 e J4 beneficiou a missão, pois aumentou consideravelmente o tempo para a aplicação do método. Por fim, como foi observado nas órbitas da região de baixa inclinação, o módulo da variação de velocidade e da aceleração de controle apresentaram valores de pequena ordem de grandeza. 4.2.2.2 Limitando e Por fim, para avaliar esta região quanto à evolução de e, elaborou-se as curvas de evolução desse parâmetro em 30 PO, conforme é apresentado na Figura 4.8. Figura 4.8 – Evolução orbital de e em 30 períodos orbitais de Vênus e i = 62◦ considerando (a) coeficientes dos harmônicos esféricos até ordem e grau 15 e (b) os coeficientes dos harmônicos zonais J2, J3 e J4. (a) (b) Fonte: Elaborado pelo próprio autor. A Figura 4.8(a) mostra a evolução considerando os coeficientes dos harmônicos esféricos até grau e ordem 15 e, por sua vez, a Figura 4.8(b) apresenta a evolução considerando os coeficientes dos harmônicos zonais J2, J3 e J4. Em ambas observa-se que a variação de e da órbita 12 é mais intensa do que das demais e, além disso, a perturbação de e é praticamente a mesma para as duas situações. Tabela 4.5 – Aplicação do método low-thrust para as órbitas da região de inclinação crítica, conside- rando ∆e = 0, 001 e ∆e = 0, 01. ∆e = 0, 001 ∆e = 0, 01 Órbita Coeficientes dos harmônicos esféricos Tempo (PO) ∆vT4 (m/s) Fϵ,min (m/s2) Tempo (PO) ∆vT4 (m/s) Fϵ,min (m/s2) 9 Até 15 6,17 3, 239173× 10−5 3, 441425× 10−9 - - - J2, J3 e J4 - - - - - - 10 Até 15 20,19 2, 251615× 10−5 1, 770423× 10−9 - - - J2, J3 e J4 28,27 2, 942318× 10−5 2, 313531× 10−9 - - - 11 Até 15 18,36 9, 037216× 10−6 5, 530892× 10−10 - - - J2, J3 e J4 20,28 7, 466439× 10−6 4, 569571× 10−10 - - - 12 Até 15 2,65 3, 544317× 10−4 1, 750122× 10−8 29,19 1, 880099× 10−4 9, 283904× 10−9 J2, J3 e J4 2,64 3, 652599× 10−4 1, 803596× 10−8 - - - Fonte: Elaborado pelo próprio autor. 52 A partir da evolução orbital de e em 30 POV e considerando o limite ∆e = 0, 001 e ∆e = 0, 01, na Tabela 4.5 são apresentados os resultados da aplicação do método de controle para as órbitas da região crítica que ultrapassaram o limite estabelecido. Nessa Tabela vê-se que para a órbita de maior valor de a (caso 12) resultou no menor tempo para a ocorrência do controle. Além