Acessibilidade / Reportar erro

Movimento Vertical de Minifoguetes: Equações de Trajetórias e Análises Gráficas

Vertical Model Rocket Movement: Trajectory Equations and Graphical Analysis

Resumos

Neste trabalho são descritas as equações diferenciais de movimento e as respectivas soluções para a velocidade [v(t)] e altitude [y(t)] em função do tempo, relativas a um foguete que segue uma trajetória vertical em um movimento de subida. São considerados quatro modelos: no primeiro é desprezada a força de arrasto do ar; no segundo é considerada uma força de arrasto linear com a velocidade; no terceiro e quarto modelos é considerado que a força de arrasto é proporcional ao quadrado da velocidade. O terceiro modelo se distingue por considerar a massa total de um foguete constante. No quarto modelo, apresenta-se uma solução completa, considerando-se que a massa decresce linearmente com o tempo durante a exaustão de gases. Para testar estes modelos foi utilizado o voo experimental do minifoguete Epsilon-8 do Grupo de Foguetes Carl Sagan da UFPR. Foram utilizados parâmetros experimentais relacionados à curva de empuxo em função do tempo (teste estático) do Epsilon-8 que permitiram estimar a velocidade de exaustão dos gases expelidos do motor. As curvas teóricas para v(t) e y(t) do terceiro e quarto modelos, por serem mais completos, foram utilizados no ajuste das respectivas curvas experimentais do Epsilon-8. O terceiro modelo, embora tenha uma solução aproximada, é razoável na previsão do apogeu, com um erro de 8,6%. O quarto modelo, por ser o mais completo, prevê o apogeu com um erro de apenas 1,5%.

Palavras-chave
Ensino; Simulação computacional; Minifoguete; Força de arrasto; Equações de movimento


This work describes the differential equations of motion and the respective solutions for speed v(t) and altitude y(t) as a function of time, relative to a rocket that follows a vertical trajectory in an upward motion were described. Four models are considered: in the first, the drag force of the air is neglected; in the second it is considered a linear drag force with speed; in the third and fourth models it is considered that the drag force is proportional to the square of the speed. The third model is distinguished by considering the total mass of a rocket constant. In the fourth model, a complete solution is presented, considering that the mass decreases linearly with time. To test these models, the experimental flight of the Epsilon-8 model rocket from the UFPR Carl Sagan Rocket Group was used. Experimental parameters related to the thrust curve as a function of the Epsilon-8 time (static test) were used, in which from the model rocket mass variation rate it was possible to obtain the rate of exhaust of the gases expelled in the engine. The third and fourth models, since they are more realistic, were used to adjust the experimental curve for v(t) and y(t). The third model, although it has an approximate solution, proves to be reasonable in the prediction of the peak, with an error of 8.6%. The fourth model, being the most complete, predicts the peak with an error of only 1.5% precisely in the time of the experimental peak.

Keywords
Teaching; Computer simulation; Model Rocket; Drag force; Motion equations


1. Introdução

As atividades de ensino utilizando projetos de foguetes, realizadas nas escolas, universidades, feiras de ciências e em competições, tornam o ambiente de aprendizado mais dinâmico e atrativo e promovem a integração entre alunos de diferentes turmas e/ou instituições. Estas atividades tornam-se desencadeadoras de ações interdisciplinares em instituições de ensino dos mais variados níveis [11. G.V.V. Gabriel e C. Rafael, em Anais do XXIV Congresso deIniciação Científica da Universidade Federal de Pelotas (Pelotas,2015), disponível em:https://wp.ufpel. edu.br/pibidfisica/files/2013/03/CE_01696_foguetes.pdf.
https://wp.ufpel. edu.br/pibidfisica/fil...
, 22. R.R. Cuzinatto, A.M. D’Ambrosio, H.F. Andrade, B.R. Duarte, V.C. Lorencetti,S.A. Maéstri, R.D. Martins e M.F.T. Filho, Revista Ciência em Extensão 11, 40 (2015).]. Durante a construção de um foguete didático, os estudantes aprendem a manusear ferramentas e a realizar medições diversas. No lançamento, muitos estudantes ficam surpresos e encantados com o voo do foguete, pela elevada aceleração e a altitude que atingem, resultados de seus esforços de trabalhos em equipe. Como toda atividade experimental, o ambiente se torna mais favorável para que os estudantes interajam mais com o professor, deixando de ser apenas espectadores, o que aprimora várias questões associadas ao aprendizado [22. R.R. Cuzinatto, A.M. D’Ambrosio, H.F. Andrade, B.R. Duarte, V.C. Lorencetti,S.A. Maéstri, R.D. Martins e M.F.T. Filho, Revista Ciência em Extensão 11, 40 (2015)., 33. L. Hennemann, Journal of Aeronautical Sciences 7, 1 (2016).].

A abordagem do foguete em sala permite explorar assuntos de diversas áreas do conhecimento. A equação diferencial do movimento é descrita com base na segunda lei de Newton. O empuxo produzido pelo motor do foguete está associado à pressão interna no mesmo que, dependendo do modelo, é ocasionada por reações químicas dos componentes do combustível [44. A.L. Alves, A.N. Paneto, K.A. Littike, S.S. Bento e C.H. Marchi, Revista Brasileira de Ensino de Física 42, e20200390 (2020).]. Já na análise de voo, diversos modelos físicos/matemáticos podem ser aplicados na equação do movimento para a obtenção de uma equação de trajetória do foguete [55. H. Rodrigues, M.O. Pinho, D. Portes Jr. e A. Santiago, International Journal of Mathematical Education in Science and Technology 40, 523 (2009)., 66. A.V. Kraff, G.S. Vázquez, R.R. Mijangos e J.A. Heredia-Cancino, Revista Mexicana de Física 61, 6 (2015).]. Desta forma, o uso da análise gráfica com o computador, no chamado ajuste da equação teórica de movimento aos dados experimentais é de extrema importância, pois possibilita aos estudantes a obtenção de parâmetros associados ao movimento do foguete. Parâmetros como o coeficiente de arrasto do ar (b), a taxa de ejeção de massa (R = dm/dt), a velocidade com que os gases são ejetados (ve) do motor do foguete, bem como dados técnicos do foguete e do tipo de combustível podem ser obtidos utilizando programas computacionais [44. A.L. Alves, A.N. Paneto, K.A. Littike, S.S. Bento e C.H. Marchi, Revista Brasileira de Ensino de Física 42, e20200390 (2020)., 77. L.K. Araki e C.H. Marchi, Sciences and Applications 1, 1 (2008).]. Quanto mais variáveis físicas são empregadas na equação diferencial de movimento e quanto mais sofisticado é o modelo que descreve esta equação, mais a solução para a equação de trajetória se aproxima dos dados experimentais do voo. Neste trabalho, desenvolve-se uma atividade teórica deste nível utilizando como foguete didático, o minifoguete.

Um minifoguete é uma versão em miniatura de um foguete real que utiliza como propulsor um motor a propelente sólido inflamável. Mostramos em um trabalho recente que o minifoguete didático pode ser construído com materiais de baixo custo adquiridos em papelaria e lojas de construção [44. A.L. Alves, A.N. Paneto, K.A. Littike, S.S. Bento e C.H. Marchi, Revista Brasileira de Ensino de Física 42, e20200390 (2020).]. A motivação de se trabalhar com o minifoguete é que este, por retratar um foguete real, pode aumentar o interesse dos estudantes pelas áreas de ciências, tecnologias e outras; e também porque ele atinge altitudes muito mais elevadas do que os que utilizam a garrafa PET [88. http://www.foguete.ufpr.br/festival/2019/festival-2019.htm, acessado em 19/10/2020.
http://www.foguete.ufpr.br/festival/2019...
], o que impressiona os estudantes. Neste trabalho, a obtenção dos dados para a construção dos gráficos experimentais de v(t) e y(t) foram adquiridas por meio de um altímetro a bordo de um minifoguete fabricado na UFPR em 2015, e que utiliza um motor-foguete comercial, conforme detalhado na seção 2.1 2.1. O minifoguete Um minifoguete, assim como um foguete real, é um veículo que se desloca expelindo, em sua parte traseira, um fluxo de gás à alta velocidade devido à queima do propelente [11]. De acordo com a terceira lei de Newton, esta ação resulta em uma reação de mesma intensidade, porém sentido oposto, deslocando o foguete para cima a partir do solo [4]. Minifoguetes são utilizados no ensino e em competições. Já outros, mais sofisticados do tipo sondagem, são aplicados às pesquisas aeroespaciais cujas altitudes são restritas à troposfera [12]. Na Figura 1(a) está ilustrada a seção transversal de um minifoguete, utilizado em atividades de ensino, com as principais partes que o compõe. Na Figura 1(b) é mostrada uma fotografia do minifoguete utilizado neste trabalho. Figura 1 (a) Ilustração do corte de seção transversal de um minifoguete. (b) Minifoguete Epsilon-8. Como se observa na Figura 1(a) a estrutura do minifoguete é delimitada pela linha em preto e inclui a coifa (nariz), a parte traseira que sustenta as empenas (ou aletas) e o motor demarcado pela linha azul. O motor contém o propelente na cor amarela que está no interior de um cilindro por onde é feito um pequeno furo em sua extremidade direita [13]. A queima do propelente é feita por um dispositivo de ignição elétrica que contém um ignitor e dois condutores. Estes ao serem ligados a uma pequena bateria transmitem energia que inflama o ignitor e consequentemente o propelente. Um dispositivo de ejeção do motor (não ilustrado aqui) desacopla o conjunto constituído pela coifa e compartimento do altímetro, o que automaticamente aciona o paraquedas. Um elástico segura todo o conjunto do foguete para um pouso seguro. O minifoguete utilizado (Figura 1(b)) para a obtenção dos dados experimentais foi o Epsilon-8. A ele foi acoplado um motor comercial de classificação E7 segundo a nomenclatura da NAR (National Association of Rocketry) [14]. Nesta nomenclatura a letra especifica que o impulso total produzido pelo motor está na faixa de 20,00 a 40,00 N.s, e o número em seguida especifica o empuxo médio produzido de 7 N. De fabricação do Grupo de Foguetes Carl Sagan (GFCS) da UFPR [8], o Epsilon-8 possui diâmetro máximo de 2,03 cm, comprimento de 50,1 cm e massa de 131,2 g [15]. Ele é o atual recordista brasileiro de apogeu para a classe E, atingindo em 2015 um apogeu de 723 m. Os dados deste voo foram obtidos com um altímetro micro-peak da Altus Metrum [9] e foram utilizados, neste trabalho, para os ajustes das equações de movimento em dados experimentais reais de v(t) e y(t). Os dados do voo foram obtidos em 08/08/2015 pelo Grupo de Foguetes Carl Sagan da Universidade Federal do Paraná (UFPR), que é liderado por Carlos Henrique Marchi, um dos autores do presente trabalho. . O altímetro é um dispositivo eletrônico que relaciona as variações de pressão com a altitude, sendo muito utilizado em competições [99. AltusMetrum, Introducing MicroPeak by AltusMetrum, disponível em:https://altusmetrum.org/MicroPeak/, acessado em 19/10/2020.
https://altusmetrum.org/MicroPeak/...
].

Embora seja prático e didático trabalhar com minifoguetes que carregam altímetros, na literatura brasileira de ensino de física este tópico é escasso. Apenas um artigo de nossa autoria foi publicado recentemente, abordando o uso do minifoguete no ensino de física [44. A.L. Alves, A.N. Paneto, K.A. Littike, S.S. Bento e C.H. Marchi, Revista Brasileira de Ensino de Física 42, e20200390 (2020).]. Por outro lado, os artigos sobre este tópico enfocam uma área técnica específica, associadas a cursos de engenharias. Em livros didáticos de física a nível de graduação, aborda-se um caso simples para a dinâmica de um foguete que geralmente despreza a força de arrasto do ar [1010. H.D. Young e R.A. Freedman, Física I(Pearson, São Paulo, 2003), v. 1.].

Neste trabalho utilizam-se quatro modelos teóricos para descrever o movimento vertical ascendente de um minifoguete: no modelo I, a equação da trajetória é obtida desprezando-se a resistência do ar; no modelo II supõe-se haver uma dependência linear da força de resistência do ar com a velocidade do foguete; no modelo III esta dependência é do tipo quadrática, entretanto, considera-se a massa do foguete como um todo constante; no modelo IV apresenta-se uma solução analítica completa para v(t) e numérica para y(t) durante a queima do propelente. Faz-se aqui uma apresentação das equações obtidas durante a queima do propelente, para os quatro modelos em estudo. Adicionalmente, apresentam-se também as soluções para as equações diferenciais de movimento, após a queima do propelente. Com isto torna-se possível a obtenção dos gráficos v(t) e y(t) desde o lançamento até o apogeu de um minifoguete. Para abordar as diferenças entre os modelos, realiza-se uma ampla discussão por meio da análise gráfica para v(t) e y(t). Em nenhum dos modelos supracitados, consideraram-se as variações da gravidade e da densidade do ar com a altitude.

Utiliza-se o voo real de um minifoguete com apogeu de 723 m, onde as soluções teóricas para as equações de movimento de cada um dos modelos podem ser confrontadas com os gráficos experimentais de v(t) e y(t). Os valores previamente calculados de b, R e ve, que compõem as equações teóricas da trajetória, podem ser utilizados como parâmetros iniciais durante as simulações computacionais. Em todos os modelos aqui estudados, o impulso do minifoguete foi considerado constante e por isso os parâmetros R e ve não foram alterados durante as simulações. Por outro lado, modificando-se os valores de b foi possível aproximar (ou ajustar) uma equação teórica, de um dado modelo, aos dados experimentais de v(t) e y(t) do minifoguete. Finalmente realiza-se uma discussão comparativa no qual conclui-se que o modelo IV se ajusta melhor aos dados experimentais.

2. Fundamentação Teórica

2.1. O minifoguete

Um minifoguete, assim como um foguete real, é um veículo que se desloca expelindo, em sua parte traseira, um fluxo de gás à alta velocidade devido à queima do propelente [1112. http://ftp.demec.ufpr.br/foguete/apostila/, acessado em 26/08/2020.
http://ftp.demec.ufpr.br/foguete/apostil...
]. De acordo com a terceira lei de Newton, esta ação resulta em uma reação de mesma intensidade, porém sentido oposto, deslocando o foguete para cima a partir do solo [44. A.L. Alves, A.N. Paneto, K.A. Littike, S.S. Bento e C.H. Marchi, Revista Brasileira de Ensino de Física 42, e20200390 (2020).]. Minifoguetes são utilizados no ensino e em competições. Já outros, mais sofisticados do tipo sondagem, são aplicados às pesquisas aeroespaciais cujas altitudes são restritas à troposfera [1213. NASA, Space, disponível em:https://www1.grc.nasa. gov/space/?ref=spaceeducation/rocket/rktcompare.html, acessado em 26/08/2019.
https://www1.grc.nasa. gov/space/?ref=sp...
]. Na Figura 1(a) está ilustrada a seção transversal de um minifoguete, utilizado em atividades de ensino, com as principais partes que o compõe. Na Figura 1(b) é mostrada uma fotografia do minifoguete utilizado neste trabalho.

Figura 1
(a) Ilustração do corte de seção transversal de um minifoguete. (b) Minifoguete Epsilon-8.

Como se observa na Figura 1(a) a estrutura do minifoguete é delimitada pela linha em preto e inclui a coifa (nariz), a parte traseira que sustenta as empenas (ou aletas) e o motor demarcado pela linha azul. O motor contém o propelente na cor amarela que está no interior de um cilindro por onde é feito um pequeno furo em sua extremidade direita [1313. NASA, Model Solid Rocket Engine, disponível em:https://www.grc.nasa.gov/WWW/K-12/rocket/rktengine.html, acessado em 26/08/2020.
https://www.grc.nasa.gov/WWW/K-12/rocket...
]. A queima do propelente é feita por um dispositivo de ignição elétrica que contém um ignitor e dois condutores. Estes ao serem ligados a uma pequena bateria transmitem energia que inflama o ignitor e consequentemente o propelente. Um dispositivo de ejeção do motor (não ilustrado aqui) desacopla o conjunto constituído pela coifa e compartimento do altímetro, o que automaticamente aciona o paraquedas. Um elástico segura todo o conjunto do foguete para um pouso seguro.

O minifoguete utilizado (Figura 1(b)) para a obtenção dos dados experimentais foi o Epsilon-8. A ele foi acoplado um motor comercial de classificação E7 segundo a nomenclatura da NAR (National Association of Rocketry) [1414. Clube de Astronomia, Espaçomodelismo, disponível em:https://www.clubedeastronomia.com.br/modelismo.php, acessado em 19/10/2020.
https://www.clubedeastronomia.com.br/mod...
]. Nesta nomenclatura a letra especifica que o impulso total produzido pelo motor está na faixa de 20,00 a 40,00 N.s, e o número em seguida especifica o empuxo médio produzido de 7 N. De fabricação do Grupo de Foguetes Carl Sagan (GFCS) da UFPR [88. http://www.foguete.ufpr.br/festival/2019/festival-2019.htm, acessado em 19/10/2020.
http://www.foguete.ufpr.br/festival/2019...
], o Epsilon-8 possui diâmetro máximo de 2,03 cm, comprimento de 50,1 cm e massa de 131,2 g [1515. Foguete UFPR, 723 metros: novo recorde de apogeu para minifoguete commotor da classe E, disponível emhttp://fogueteufpr.blogspot.com/2015/12/723-metros-novo-recorde-de-apogeu-para.html,acessado em 19/10/2020.
http://fogueteufpr.blogspot.com/2015/12/...
]. Ele é o atual recordista brasileiro de apogeu para a classe E, atingindo em 2015 um apogeu de 723 m. Os dados deste voo foram obtidos com um altímetro micro-peak da Altus Metrum [99. AltusMetrum, Introducing MicroPeak by AltusMetrum, disponível em:https://altusmetrum.org/MicroPeak/, acessado em 19/10/2020.
https://altusmetrum.org/MicroPeak/...
] e foram utilizados, neste trabalho, para os ajustes das equações de movimento em dados experimentais reais de v(t) e y(t). Os dados do voo foram obtidos em 08/08/2015 pelo Grupo de Foguetes Carl Sagan da Universidade Federal do Paraná (UFPR), que é liderado por Carlos Henrique Marchi, um dos autores do presente trabalho.

2.2. A Equação diferencial de movimentoe o empuxo

Considere um sistema S constituído unicamente por um foguete, e um sistema S’ cuja fronteira separa o foguete de seu meio exterior (Figura 2(a)). Este meio contribui com as forças externas de atração gravitacional que a Terra exerce sobre o foguete (P = mg) e a força de arrasto do ar (Fr). A força de empuxo T, que acelera o foguete para cima (valores de y crescente) é uma força interna ao sistema S’ e surge devido ao fato de o foguete estar exaurindo gases a altas velocidades, proveniente da queima do propelente em seu motor (Figura 2(b)). Esta é uma força de reação devido à 3a lei de Newton. Na Figura 2(c) está representado um diagrama de corpo livre para o sistema S, que inclui tanto as forças externas quanto a interna. Por simplicidade, são descritas nos próximos parágrafos as grandezas físicas vetoriais apenas em termos de suas componentes na direção y.

Figura 2
(a) Foguete a uma velocidade v em um instante de tempo t. (b) Ao ejetar uma quantidade de massa de propelente dm, sua velocidade no instante t+dt será v + dv. (c) Diagrama de corpo livre.

Na descrição da equação de movimento, considera-se o foguete subindo a uma velocidade v, de tal forma que a magnitude do momento linear inicial do foguete no instante t é Po = mv (Figura 2(a)). No instante t+dt o foguete expele uma massa infinitesimal dm (negativa) de tal forma que a massa e a velocidade do foguete são, respectivamente, m+dm e v+dv. A velocidade em que dm é exaurida é dada em relação a um referencial inercial em repouso na Terra por v + ve, no qual ve é a velocidade de exaustão dos gases em relação ao foguete. Desta forma, o momento linear total do sistema S’ no instante t+dt é: Pf = (m + dm)(v + dv) + (−dm)(vve). A segunda lei de Newton aplicada ao sistema S’ possibilitará obter uma expressão para a equação do foguete:

(1) F e x t = d P d t .

Na Equação (1) dP representa a variação infinitesimal do momento linear entre os instantes t e t+dt. Conforme descrito nos parágrafos anteriores, esta variação é dada por: Pf-Po = (m + dm)(v + dv) + (−dm)(vve)−mv. Desenvolvendo estes termos, desprezando o produto dm.dv e lembrando que para o sistema S’, ∑Fext = −mg + Fr (com Fr negativo) a Equação (1) se torna:

(2) - m g + F r = m d v + v e d m d t ,

que é equivalente a:

(3) m d v d t = - v e d m d t - m g + F r .

O membro esquerdo da Equação (3) representa a massa remanescente, m = m(t) do sistema S, vezes a aceleração do mesmo [44. A.L. Alves, A.N. Paneto, K.A. Littike, S.S. Bento e C.H. Marchi, Revista Brasileira de Ensino de Física 42, e20200390 (2020).]. Esse termo é interpretado como o somatório de todas as forças que atuam no sistema do foguete. Já o primeiro termo no membro direito desta equação é definido como a magnitude do empuxo, T. Tendo em vista que dm/dt é negativo, o empuxo é positivo na Equação (3).

Na prática, T é função do tempo e um gráfico T(t) pode ser obtido experimentalmente a partir de testes estáticos do motor do foguete. Entretanto, para algumas estimativas teóricas para a solução da Equação (3), T pode ser considerado constante, o que é válido para um caso teórico no regime estacionário do motor [44. A.L. Alves, A.N. Paneto, K.A. Littike, S.S. Bento e C.H. Marchi, Revista Brasileira de Ensino de Física 42, e20200390 (2020).]. Neste caso, tanto ve quanto dm/dt assumem valores constantes durante o voo do foguete. Apresentam-se na próxima seção as soluções da equação (3) para a obtenção de v(t) e y(t) em quatro modelos abordados, considerando-se que ve e dm/dt são constantes.

2.3. Soluções Analíticas para a Equaçãode Movimento de um Foguete

Nas seções 2.3.1 2.3.1. Modelo I: Movimento sem a forçade arrasto do ar Quando a força de arrasto do ar é desprezada, a Equação (3) se torna: (4) m ⁢ d ⁢ v d ⁢ t = - v e ⁢ d ⁢ m d ⁢ t - m ⁢ g . Assumindo que dm/dt é uma constante nesta equação, então a massa diminui linearmente com o tempo: (5) m ⁢ ( t ) = m o - R ⁢ t , na qual R = −dm/dt (lembrando que dm/dt é negativo) e mo=mp + mf é a massa inicial do foguete, com mp denotando a massa total do propelente e mf a massa final do foguete [5, 6]. Substituindo a Equação (5) na (4), obtém-se: (6) d ⁢ v d ⁢ t = R ⁢ v e m o - R ⁢ t - g . As soluções da Equação (6) para v(t) e y(t), são dadas por: (7) v I ⁢ A ⁢ ( t ) = - g ⁢ t - v e ⁢ l ⁢ n ⁢ [ 1 - R m o ⁢ t ] , (8) y I ⁢ A ⁢ ( t ) = v e ⁢ t - 1 2 ⁢ g ⁢ t 2 + m o - R ⁢ t R ⁢ v e ⁢ l ⁢ n ⁢ [ 1 - R m o ⁢ t ] . Na etapa B, após a queima de todo o propelente (t≥tq), a Equação (3) se torna: (9) d ⁢ v d ⁢ t = - g , cujas soluções para v(t)e y(t) são, respectivamente: (10) v I ⁢ B = v q - g ⁢ ( t - t q ) , (11) y I ⁢ B = y q + v q ⁢ ( t - t q ) - g 2 ⁢ ( t - t q ) 2 , nas quais vq é a velocidade do foguete no tempo tq; yq é altura do foguete no tempo tq. De outra forma: vq = vIA(tq) e yq = yIA(tq). a 2.3.4 2.3.4. Modelo IV: Solução completa com forçade arrasto quadrática Uma solução analítica completa da Equação (24) para v(t), durante a queima do propelente, foi obtida por Rodrigues et al. (2014) [5]. O autor demonstrou, utilizando técnicas de cálculo diferencial, que a Equação (3) pode ser reescrita em termos de uma equação diferencial cujas soluções são funções de Bessel. Veja inicialmente que a Equação (24), agora com m = mo−Rt, pode ser escrita como: (31) d ⁢ v d ⁢ t + b m ⁢ v 2 - R ⁢ v e m + g = 0 , Durante as etapas para a solução da Equação (31), é possível realizar uma mudança da variável v, para uma nova variável χ [5]. Neste caso, esta equação se transforma em: (32) ξ 2 ⁢ d 2 ⁢ χ d ⁢ t 2 + ξ ⁢ d ⁢ χ d ⁢ t + ( ξ 2 - ν 2 ) ⁢ χ = 0 , que é uma Equação de Bessel de ordem ν [5, 19]. As constantes ξ e ν são definidas como: (33) ξ = ξ ⁢ ( m ) = 2 R ⁢ b ⁢ m ⁢ g , (34) ν = 2 ⁢ b ⁢ v e R . A solução geral da Equação (32) é uma combinação linear do tipo: (35) χ ⁢ ( ξ ) = C 1 ⁢ J ν ⁢ ( ξ ) + C 2 ⁢ Y ν ⁢ ( ξ ) , no qual Jν(ξ) e Yν(ξ) são funções de Bessel de primeira e segunda ordens respectivamente [5, 19, 20], C1 e C2 são constantes reais. Após algumas álgebras utilizando propriedades convenientes das funções de Bessel [5], é possível escrever a solução geral para a velocidade como: (36) v I ⁢ V ⁢ A ⁢ ( ξ ) = R 2 ⁢ b ⁢ [ ξ ⁢ C ⁢ Y ν + 1 ⁢ ( ξ ) + J ν + 1 ⁢ ( ξ ) C ⁢ Y ν ⁢ ( ξ ) + J ν ⁢ ( ξ ) - v o ] , (37) c ⁢ o ⁢ m ⁢ C = ( 2 ⁢ b ⁢ v o + R ⁢ ν ) ⁢ J ν ⁢ ( ξ o ) - R ⁢ ξ o ⁢ J ν + 1 ⁢ ( ξ o ) R ⁢ ξ o ⁢ Y ν + 1 ⁢ ( ξ o ) - ( 2 ⁢ b ⁢ v o + R ⁢ ν ) ⁢ Y ν ⁢ ( ξ ) na qual vo é a velocidade inicial do foguete, que pode ser definida como zero e ξo=2R⁢b⁢mo⁢g. Tendo em vista que ξ = ξ(m), Equação (36), e que m = m(t), Equação (5), é possível por meio da Equação (36) construir um gráfico da velocidade em função do tempo para (0≤t≤tq), etapa A, conhecendo-se os parâmetros mo,g,R,b e ve, previamente obtidos de cálculos teóricos/experimentais. A construção do gráfico y(t) pode ser obtida a partir de v(t) utilizando programas de computador adequados. Para este modelo, foi utilizado o programa GNU Octave [21] para simular as funções v(t) e y(t); a partir da velocidade foi obtida a altura através de um processo de integração numérica denominado quadratura de Gauss. No Apêndice I estão descritas as funções utilizadas no GNU Octave, bem como os scriptsimplementados para simulação do modelo apresentado nesta seção. No Apêndice II foi deixada uma explicação sucinta sobre a quadratura de Gauss. Para a etapa B neste modelo as Equações (28) e (29) podem ser aplicadas. Nesta condição deve-se considerar que m = m(tq). apresentam-se as soluções para a equação de movimento de um foguete durante o tempo de queima (tq) do propelente, fase propulsada, que aqui denomina-se de etapa A. Após este tempo, na fase balística, denominam-se as soluções como etapa B. Nesta Etapa, após o tempo tq, o motor deixa de funcionar e a magnitude do empuxo T, na Equação (3), se torna nula. Nesta etapa, o foguete segue um voo no qual apenas as forças de resistência do ar e o peso estão presentes. Utiliza-se a notação W(t)nA,B para designar uma equação de trajetória (soluções da Equação (3)). Nesta notação, W(t) representa v(t) ou y(t), n representa o modelo (I, II …) e A ou B, uma das etapas de voo. No modelo I a seguir, coloca-se a solução para o movimento do foguete, desconsiderando a força de arrasto do ar, com a intenção de que o estudante aprecie as soluções e compare com modelos que levam em consideração esta força. Este movimento representaria o movimento de um foguete lançado da superfície de um astro sem atmosfera, como a Lua.

Os estudantes são encorajados a verificar as deduções e soluções nas referências deixadas, para cada uma das equações de trajetória aqui apresentadas. Para todas as equações aqui colocadas, considerou-se que as velocidades e alturas iniciais (t = 0s) possuem, ambas, valores nulos.

2.3.1. Modelo I: Movimento sem a forçade arrasto do ar

Quando a força de arrasto do ar é desprezada, a Equação (3) se torna:

(4) m d v d t = - v e d m d t - m g .

Assumindo que dm/dt é uma constante nesta equação, então a massa diminui linearmente com o tempo:

(5) m ( t ) = m o - R t ,

na qual R = −dm/dt (lembrando que dm/dt é negativo) e mo=mp + mf é a massa inicial do foguete, com mp denotando a massa total do propelente e mf a massa final do foguete [55. H. Rodrigues, M.O. Pinho, D. Portes Jr. e A. Santiago, International Journal of Mathematical Education in Science and Technology 40, 523 (2009)., 66. A.V. Kraff, G.S. Vázquez, R.R. Mijangos e J.A. Heredia-Cancino, Revista Mexicana de Física 61, 6 (2015).]. Substituindo a Equação (5) na (4), obtém-se:

(6) d v d t = R v e m o - R t - g .

As soluções da Equação (6) para v(t) e y(t), são dadas por:

(7) v I A ( t ) = - g t - v e l n [ 1 - R m o t ] ,
(8) y I A ( t ) = v e t - 1 2 g t 2 + m o - R t R v e l n [ 1 - R m o t ] .

Na etapa B, após a queima de todo o propelente (ttq), a Equação (3) se torna:

(9) d v d t = - g ,

cujas soluções para v(t)e y(t) são, respectivamente:

(10) v I B = v q - g ( t - t q ) ,
(11) y I B = y q + v q ( t - t q ) - g 2 ( t - t q ) 2 ,

nas quais vq é a velocidade do foguete no tempo tq; yq é altura do foguete no tempo tq. De outra forma: vq = vIA(tq) e yq = yIA(tq).

2.3.2. Modelo II: Movimento com força de arrasto linear

A força de arrasto do ar depende da velocidade v de um objeto e geralmente é descrita em termos de potência n da velocidade, ou seja,

(12) F r = - b v n ,

na qual b é uma constante positiva que especifica o poder da força de arrasto [1818. Mecânica Newtoniana: Partícula Única,disponível em:http://jararaca.ufsm.br/websites/petfisica/download/Arquivos/Capitulo%202.pdf, acessado em 19/10/2020.
http://jararaca.ufsm.br/websites/petfisi...
]. Existem dois regimes diferentes no qual essa dependência pode ser do tipo linear ou quadrática em v.O parâmetro que diferencia estes regimes é o número de Reynolds (nRe) [44. A.L. Alves, A.N. Paneto, K.A. Littike, S.S. Bento e C.H. Marchi, Revista Brasileira de Ensino de Física 42, e20200390 (2020).]. Experimentalmente, para um objeto relativamente pequeno movendo-se no ar, n=1 na Equação (12) para velocidades de até 24 m/s. Para velocidades mais elevadas, a força de arrasto é aproximadamente proporcional ao quadrado da velocidade [1818. Mecânica Newtoniana: Partícula Única,disponível em:http://jararaca.ufsm.br/websites/petfisica/download/Arquivos/Capitulo%202.pdf, acessado em 19/10/2020.
http://jararaca.ufsm.br/websites/petfisi...
]. As constantes b’s na Equação (12) se diferem por ordem de grandeza e unidade de medição. Nesta seção frisa-se mais na obtenção das equações de movimento e, por simplicidade, utilizou-se a mesma letra, sem subíndice, para representar esta variável. Na seção 3 3. Resultados Teóricos Com base nas Equações desenvolvidas na seção 2.3 foram construídos gráficos teóricos para os modelos de I a IV desta mesma seção. Foram utilizados os parâmetros descritos na Tabela 1 e aqueles obtidos a partir dela (ve e R). Nas Figuras 4(a)–(c) à esquerda estão apresentados gráficos de v(t) para os modelos de I a IV e, nas Figuras 4(d)–(f) à direita os respectivos gráficos de y(t). Os gráficos se distinguem por diferentes constantes de arrasto, b, que estão discriminados por diferentes cores. Quando a força de arrasto tiver uma dependência linear, a constante b será designada por bL e quando tiver uma dependência quadrática, será designada por bQ. Estas constantes se diferem pela unidade de medição: bL é dada em kg/s e bQ em kg/m. Figura 4 (a)–(c) Curvas de v(t) para os modelos de I a IV. O modelo I está representado pela curva em preto. Nas partes (d)–(f), à direita, os respectivos gráficos de y(t). Os valores destas constantes foram escolhidos para que as curvas de velocidade e altitude pudessem ser comparadas entre os diferentes modelos. Para efeitos de comparações, nas Figuras 4(a) e 4(d), na cor em preto, foram colocados os gráficos para v(t) e y(t) com b=0 relativo ao modelo I, juntamente com os gráficos referentes ao modelo II. Todos os gráficos foram plotados até seus respectivos tempos de apogeu, ta. Nos gráficos de v(t) este tempo se refere à velocidade nula para t > 0 de cada curva. Já nos gráficos de y(t), para os valores de t = ta, as curvas atingem a altitude máxima (apogeu). Todos os gráficos de velocidade são caracterizados por duas partes: a primeira nas quais as velocidades aumentam com o tempo, deste t = 0s ate t = tq = 5,358 s e a segunda nas quais as velocidades decrescem instantaneamente a partir deste tempo, tornando-se nulas em seus respectivos ta′⁢s. Na primeira parte do modelo I, por exemplo, as curvas v(t) e y(t) são descritas pelas Equações (7) e (8) para vIA(t) e yIA(t), respectivamente. Após este tempo as curvas são descritas pelas Equações (10) e (11) para vIB(t) e yIB(t), respectivamente. A mesma sistemática é válida para as demais curvas. Os picos nos gráficos de v(t) podem ser explicados como se segue. Antes de tq o termo de empuxo T=-ve⁢d⁢md⁢t supera os da soma (mg + |Fr|) na Equação (3), o que confere uma aceleração positiva ao foguete. Teoricamente após este tempo o termo T torna-se instantaneamente nulo e apenas as forças retardatárias mg e Fr prevalecem. Isto confere uma mudança rápida no sentido da aceleração do foguete e assim a velocidade passa a decrescer instantaneamente. Vale lembrar que em tq a força Fr atinge seu valor máximo (pois a velocidade é máxima) o que pode contribuir ainda mais com os picos característicos nos gráficos de v(t). Nos gráficos y(t) surge um ponto de inflexão em t = tq no qual a velocidade de um foguete passa a decrescer. Observa-se nos gráficos da Figura 1(a) e 1(d) que, como esperado, os valores de velocidades e altitudes para b=0 (modelo I) são superiores ao de qualquer um dos modelos que apresentam b≠ 0. Ainda, como esperado para os modelos de II a IV, um aumento nos valores de b acarreta em decréscimos nas curvas de v(t) e y(t). As curvas com bL = 0,001 kg/s (modelo II) e bQ = 0,00005 kg/m (modelo III), que apresentam menores valores desta grandeza, são as que mais se aproximam da curva do modelo I. A velocidade para o modelo I em tq é da ordem de 271 m/s e supera a do modelo II com bL = 0,001 kg/s em aproximadamente 6 m/s. Após este tempo, a ausência do arrasto no modelo I confere o maior valor de ta e assim maior altitude (≈4.463,2 m) excedendo a do modelo II com bL = 0,001 kg/s em aproximadamente 801,5 m. Ainda, as curvas de v(t) para este modelo, para tempos até tq, apresentam diferentes formas conforme os valores de b aumentam. Para valores de bL ≈ 0,001 kg/s as derivadas em cada ponto da curva aumentam com o tempo, indicando um aumento da aceleração do foguete com o tempo. Conforme os valores de bL se tornam maiores as curvas mudam de tal forma que as derivadas em cada ponto passam a diminuir, indicando um decréscimo da aceleração do foguete com o tempo. Estes são os casos das curvas para bL = 0,01, 0,03 e 0,05 kg/s. Para que os gráficos de velocidade dos modelos III (Figura 4(b)) e IV (Figura 4(c)) estejam na faixa de 0–300 m/s (para se realizar uma comparação com os modelos I e II) é necessário que os valores de bQ sejam inferiores a uma potência de pelo menos 10−2. Por representarem modelos mais realísticos para movimento de foguetes, estes modelos puderam ser simulados com a mesma faixa de valores da constante bQ. Inicialmente, aparenta-se não haver distinção entre as formas das curvas, considerando um mesmo valor de b dentre esses modelos. Observa-se, entretanto, que os valores de velocidade nos tempos tq e ta são ligeiramente diferentes, o que acarreta em apogeus distintos para estes modelos (Figuras 4(e) e (f)). Para efeitos de comparação, são apresentadas na Figura 5(a) as curvas v(t) para os modelos III e IV com os valores de bQ = 0,00005 kg/m e 0,0002 kg/m, e na Figura 5(b), as respectivas curvas y(t). Observa-se que, apesar da semelhança marcante entre as curvas (para ambos os valores de bQ), os valores de velocidade para o modelo III são superiores aos do modelo IV em todos os instantes de tempo e, como resultado disto, as altitudes para este modelo são superiores às do modelo IV. A mesma análise é válida para as curvas com bQ = 0,0001 kg/m e 0,0005 kg/m. Figura 5 Comparação para as curvas (a) v(t) e (b) y(t) para os modelos III e IV. Foram utilizados valores de bQ = 0,00005 kg/m e 0,0002 kg/m. Os modelos II, III e IV supracitados parecem todos serem úteis para descrever o movimento vertical ascendente de foguetes. Entretanto, sabe-se que o modelo II não se adequa ao movimento de foguete, pois este modelo se aplica a objetos que se deslocam em meios viscosos e adaptável para situações nas quais as velocidades são inferiores a 24 m/s [18]. O minifoguete Epsilon-8 atinge uma velocidade ≈128,3 m/s até o instante tq e como será visto, o modelo I não é muito ideal para simular o movimento deste minifoguete. São feitas, na próxima seção, análises comparativas com um voo real (experimental) do minifoguete Epsilon-8, abordado nas seções 2.1 e 2.4. , faz-se melhor a distinção das constantes b’s dos modelos II, III e IV.

Assumindo para o modelo II uma dependência linear para Fr, a Equação (3) se torna:

(13) m d v d t = - v e d m d t - m g - b v .

Substituindo a Equação (5) na Equação (13) e reagrupando os termos, obtém-se:

(14) d v d t + b m o - R t v = v e m o - R t R - g .

As soluções analíticas para v(t) e y(t) são dadas pelas Equações (15) e (16) [66. A.V. Kraff, G.S. Vázquez, R.R. Mijangos e J.A. Heredia-Cancino, Revista Mexicana de Física 61, 6 (2015).]:

(15) v I I A ( t ) = v e R b + g R - b ( m o - R t ) + c 1 ( m o - R t ) b R ,

em que:

(16) c 1 = - v e R ( R - b ) - g m o b R ( R - b ) m o b R ;
(17) y I I A ( t ) = v e R b t - g 2 R ( R - b ) ( m o - R t ) 2 - c 1 R + b ( m o - R t ) 1 + b R + c 2 ,

em que:

(18) c 2 = g m o 2 2 R ( R - b ) + c 1 m o 1 + b R R + b .

Após este tempo, na etapa B, a Equação (3) se torna:

(19) m f d v d t = - m f g - b v ,

na qual mf é a massa final do foguete após a queima de todo o propelente no instante tq. As soluções da Equação (19) para v(t) e y(t) podem ser resolvidas por uma integração direta. As soluções são, respectivamente,

(20) v I I B ( t ) = [ ( m f g b + v q ) e - b m f ( t - t q ) - m f g b ]
(21) y I I B ( t ) = y q - m f b g ( t - t q ) + m f 2 b 2 g ( 1 + b m f g v q ) y I I B ( t ) × ( 1 - e - b m f ( t - t q ) ) ,

nas quais vq = vIIA(tq) e yq = yIIA(tq).

2.3.3. Modelo III: Solução aproximada commassa constante e força de arrastoquadrática

Um modelo mais realístico do movimento de um foguete, durante sua subida, considera o termo n=2 na Equação (12). O Epsilon-8, por exemplo, atinge velocidades ≈100 m/s em um intervalo de tempo relativamente curto ≈3 s. Neste caso, nRe≥103 e a dependência quadrática deve ser levada em consideração [44. A.L. Alves, A.N. Paneto, K.A. Littike, S.S. Bento e C.H. Marchi, Revista Brasileira de Ensino de Física 42, e20200390 (2020).]. Como aproximação considera-se um valor médio para a massa do foguete (m), constante durante a queima do propelente [44. A.L. Alves, A.N. Paneto, K.A. Littike, S.S. Bento e C.H. Marchi, Revista Brasileira de Ensino de Física 42, e20200390 (2020)., 1616. M. Abate, E. Anandapadmanaban, L. Bao, S. Challani, J. Gaughan, A. Jiang, A.Lingineni, A. Vora, C. Yang e D. Zhao, Correlation Between Simulated, Calculated, and Measured Model Rocket Flight, disponível em:http://www.drew.edu/govschool/wp-content/uploads/sites/99/T1-Final-Paper.pdf, acessado em 26/08/2020.
http://www.drew.edu/govschool/wp-content...
, 1717. J. Waters, Undergraduate Journal of Mathematical: Modeling onde + Two 6, 1 (2014).]. Assim, adota-se também m constante durante todo o voo do foguete. Esta aproximação é válida se o tempo de queima é pequeno quando comparado ao tempo para o foguete atingir o apogeu. Neste caso, m é dado pela média:

(22) m = m o + m f 2 ,

na qual mo e mf são as respectivas massa inicial e final do minifoguete. A equação (3) escrita na forma:

(23) m d v d t = - v e d m d t - m g - b v 2 ,

pode ser reescrita como:

(24) m β - b v 2 d v = d t ,

na qual: β = -vedmdt-mg. As soluções da Equação (24) para v(t) e y(t), na etapa A, descritas por Waters (2014) [1717. J. Waters, Undergraduate Journal of Mathematical: Modeling onde + Two 6, 1 (2014).] e relatadas por Alves et al. [44. A.L. Alves, A.N. Paneto, K.A. Littike, S.S. Bento e C.H. Marchi, Revista Brasileira de Ensino de Física 42, e20200390 (2020).], são dadas respectivamente por:

(25) v I I I A ( t ) = β b t g h [ ( b β m ) t ] ,
(26) y I I I A ( t ) = m b l n [ cosh ( b β m t ) ] .

As Equações (25) e (26) são válidas para β> 0. Esta condição é válida para o Epsilon-8 (veja a seção 2.4 2.4. Parâmetros teóricos/experimentaisrelacionados ao minifoguete Epsilon-8 Como descrito na seção 2.1 foi utilizado o minifoguete Epsilon-8 com motor a propelente sólido classe E7. Parâmetros de importância sobre as características do motor classe E7 podem ser obtidos experimentalmente através de testes estáticos [22]. Nestes testes, o motor é acoplado a um anteparo que ligado a um sistema eletrônico, fornece a força/empuxo em função do tempo [T(t)], quando o motor é ignitado. A Figura 3 representa um dos testes estáticos realizado, dentre 6 obtidos pelo GFCS da UFPR. O valor do empuxo médio, representado em vermelho, foi de 6,58 N. Figura 3 Curva experimental de empuxo em função do tempo para o motor E7 utilizado no minifoguete. No Gráfico da Figura 3, com os pontos experimentais representados por círculos, a magnitude do impulso total (Itot) é definida de acordo com a Equação (36) como a integral de [T(t)], durante o tempo de funcionamento do motor, ou seja, até o tempo final de queima (tq) do propelente: (38) I t ⁢ o ⁢ t = ∫ 0 t q T ⁢ ( t ) ⁢ d t . A magnitude do empuxo médio (Tm) que o motor exerce sobre o minifoguete é calculado pela Equação: (39) T m = I t ⁢ o ⁢ t t q . As equações de movimento descritas nas seções 2.3.1 a 2.3.4 consideram um empuxo constante, igual ao valor médio do empuxo real calculado através da Equação (39) e representado pela linha vermelha na Figura 3. Esta aproximação é razoável no chamado regime estacionário e quando os regimes transientes (início e término do funcionamento do motor) são breves [23]. Assim, a magnitude do impulso total pode ser calculada como: (40) I t ⁢ o ⁢ t = v e ⁢ | d ⁢ m d ⁢ t | ⁢ t q . Outro parâmetro relevante é o impulso específico (Isp), obtido pela divisão do impulso total pelo peso do propelente (mpg) de acordo com a equação: (41) I s ⁢ p = I t ⁢ o ⁢ t m p ⁢ g , na qual g é a magnitude da aceleração da gravidade ao nível do mar. Com esta grandeza mede-se a eficiência do propelente, já que, quanto maior for o impulso específico, maior será a variação da quantidade de movimento, por unidade de peso do propelente. Admite-se neste trabalho que o empuxo T(t) é constante, [4]. Desta forma, considerando-se que ve e dm/dt são constantes, é possível estimar a velocidade de exaustão dos gases de acordo com a equação: (42) v e = I t m p = I s ⁢ p ⁢ g . Mede-se experimentalmente tq e mp sendo possível estimar R = dm/dt, como: (43) R = m p / t q . A Tabela 1 fornece a massa do propelente utilizado no minifoguete Epsilon-8 e os parâmetros médios relevantes obtidos dos testes estáticos do motor. Nesta tabela também se encontram outros parâmetros experimentais relacionados ao minifoguete e ao voo do mesmo. Outros parâmetros como velocidade máxima e altitude máxima atingida serão discutidos na análise gráfica da seção 3. Tabela 1 Parâmetros do motor do minifoguete Epsilon-8. A grandeza mp é a massa do propelente, It é o valor médio do impulso total, obtido do teste estático de 6 motores idênticos e tq é o tempo de queima do propelente; moé a massa inicial do minifoguete (antes da queima do propelente), mfé a massa final do mesmo e taé o tempo necessário para atingir o apogeu. mp(g) It(médio)(N.s) tq (s) mo (g) mf(g) ta(s) 44,70 ± 0,01 34,90 ± 0,01 5,358 ± 0,001 131,20 ± 0,01 86,50 ± 0,01 11,7 ± 0,1 Utilizando a Equação (42) e os valores de It (médio) e mp da Tabela 1 é possível estimar ve = (780,8 ± 0,4) m/s. Adicionalmente, com os valores de mp e tq da Tabela 1 é possível estimar, utilizando a Equação (43), R = (8,343 ± 0,003) g/s. ). Para a etapa B, o empuxo é nulo. Logo, a Equação (23) torna-se:

(27) m d v d t = - m g - b v 2 ,

cujas soluções para v(t) e y(t) são, respectivamente:

(28) v I I I B ( t ) = m g b t g [ b g m ( t q - t ) + t g - 1 ( b m g v q ) ] ,
(29) y I I I B ( t ) = y q + m b l n { c o s [ b g m ( t a - t ) ] c o s [ b g m ( t a - t q ) ] } .

Na Equação (28) tem-se que vq = vIIIA(tq) e, na Equação (29), yq = yIIIA(tq), em que ta é o instante no qual o foguete atinge o apogeu. Este instante é dado pela Equação:

(30) t a = t q + m b g t g - 1 ( b m g v q ) .

2.3.4. Modelo IV: Solução completa com forçade arrasto quadrática

Uma solução analítica completa da Equação (24) para v(t), durante a queima do propelente, foi obtida por Rodrigues et al. (2014) [55. H. Rodrigues, M.O. Pinho, D. Portes Jr. e A. Santiago, International Journal of Mathematical Education in Science and Technology 40, 523 (2009).]. O autor demonstrou, utilizando técnicas de cálculo diferencial, que a Equação (3) pode ser reescrita em termos de uma equação diferencial cujas soluções são funções de Bessel. Veja inicialmente que a Equação (24), agora com m = moRt, pode ser escrita como:

(31) d v d t + b m v 2 - R v e m + g = 0 ,

Durante as etapas para a solução da Equação (31), é possível realizar uma mudança da variável v, para uma nova variável χ [55. H. Rodrigues, M.O. Pinho, D. Portes Jr. e A. Santiago, International Journal of Mathematical Education in Science and Technology 40, 523 (2009).]. Neste caso, esta equação se transforma em:

(32) ξ 2 d 2 χ d t 2 + ξ d χ d t + ( ξ 2 - ν 2 ) χ = 0 ,

que é uma Equação de Bessel de ordem ν [55. H. Rodrigues, M.O. Pinho, D. Portes Jr. e A. Santiago, International Journal of Mathematical Education in Science and Technology 40, 523 (2009)., 1919. R. Teixeira, H. Medeiros e M.L. Menezes, Séries de Fourier eEquações Diferenciais Parciais, disponível em:http://www.professores.uff.br/freddyh/wp-content/uploads/sites/105/2017/08/notasRalph.pdf,acessado em 06/11/2020.
http://www.professores.uff.br/freddyh/wp...
]. As constantes ξ e ν são definidas como:

(33) ξ = ξ ( m ) = 2 R b m g ,
(34) ν = 2 b v e R .

A solução geral da Equação (32) é uma combinação linear do tipo:

(35) χ ( ξ ) = C 1 J ν ( ξ ) + C 2 Y ν ( ξ ) ,

no qual Jν(ξ) e Yν(ξ) são funções de Bessel de primeira e segunda ordens respectivamente [55. H. Rodrigues, M.O. Pinho, D. Portes Jr. e A. Santiago, International Journal of Mathematical Education in Science and Technology 40, 523 (2009)., 1919. R. Teixeira, H. Medeiros e M.L. Menezes, Séries de Fourier eEquações Diferenciais Parciais, disponível em:http://www.professores.uff.br/freddyh/wp-content/uploads/sites/105/2017/08/notasRalph.pdf,acessado em 06/11/2020.
http://www.professores.uff.br/freddyh/wp...
, 2020. E. Brietzke, Funções de Bessel, disponível em:www.mat.ufrgs.br/~brietzke/bes1/bes1.html, acessado em 06/11/2020.
www.mat.ufrgs.br/~brietzke/bes1/bes1.htm...
], C1 e C2 são constantes reais. Após algumas álgebras utilizando propriedades convenientes das funções de Bessel [55. H. Rodrigues, M.O. Pinho, D. Portes Jr. e A. Santiago, International Journal of Mathematical Education in Science and Technology 40, 523 (2009).], é possível escrever a solução geral para a velocidade como:

(36) v I V A ( ξ ) = R 2 b [ ξ C Y ν + 1 ( ξ ) + J ν + 1 ( ξ ) C Y ν ( ξ ) + J ν ( ξ ) - v o ] ,
(37) c o m C = ( 2 b v o + R ν ) J ν ( ξ o ) - R ξ o J ν + 1 ( ξ o ) R ξ o Y ν + 1 ( ξ o ) - ( 2 b v o + R ν ) Y ν ( ξ )

na qual vo é a velocidade inicial do foguete, que pode ser definida como zero e ξo=2Rbmog.

Tendo em vista que ξ = ξ(m), Equação (36), e que m = m(t), Equação (5), é possível por meio da Equação (36) construir um gráfico da velocidade em função do tempo para (0≤ttq), etapa A, conhecendo-se os parâmetros mo,g,R,b e ve, previamente obtidos de cálculos teóricos/experimentais. A construção do gráfico y(t) pode ser obtida a partir de v(t) utilizando programas de computador adequados. Para este modelo, foi utilizado o programa GNU Octave [2121. J.W. Eaton, D. Bateman, S. Hauberg e R. Wehbring, GNU Octave a high-level interactive language for numerical computations (2020), v. 5.2.0.] para simular as funções v(t) e y(t); a partir da velocidade foi obtida a altura através de um processo de integração numérica denominado quadratura de Gauss. No Apêndice I estão descritas as funções utilizadas no GNU Octave, bem como os scriptsimplementados para simulação do modelo apresentado nesta seção. No Apêndice II foi deixada uma explicação sucinta sobre a quadratura de Gauss. Para a etapa B neste modelo as Equações (28) e (29) podem ser aplicadas. Nesta condição deve-se considerar que m = m(tq).

2.4. Parâmetros teóricos/experimentaisrelacionados ao minifoguete Epsilon-8

Como descrito na seção 2.1 2.1. O minifoguete Um minifoguete, assim como um foguete real, é um veículo que se desloca expelindo, em sua parte traseira, um fluxo de gás à alta velocidade devido à queima do propelente [11]. De acordo com a terceira lei de Newton, esta ação resulta em uma reação de mesma intensidade, porém sentido oposto, deslocando o foguete para cima a partir do solo [4]. Minifoguetes são utilizados no ensino e em competições. Já outros, mais sofisticados do tipo sondagem, são aplicados às pesquisas aeroespaciais cujas altitudes são restritas à troposfera [12]. Na Figura 1(a) está ilustrada a seção transversal de um minifoguete, utilizado em atividades de ensino, com as principais partes que o compõe. Na Figura 1(b) é mostrada uma fotografia do minifoguete utilizado neste trabalho. Figura 1 (a) Ilustração do corte de seção transversal de um minifoguete. (b) Minifoguete Epsilon-8. Como se observa na Figura 1(a) a estrutura do minifoguete é delimitada pela linha em preto e inclui a coifa (nariz), a parte traseira que sustenta as empenas (ou aletas) e o motor demarcado pela linha azul. O motor contém o propelente na cor amarela que está no interior de um cilindro por onde é feito um pequeno furo em sua extremidade direita [13]. A queima do propelente é feita por um dispositivo de ignição elétrica que contém um ignitor e dois condutores. Estes ao serem ligados a uma pequena bateria transmitem energia que inflama o ignitor e consequentemente o propelente. Um dispositivo de ejeção do motor (não ilustrado aqui) desacopla o conjunto constituído pela coifa e compartimento do altímetro, o que automaticamente aciona o paraquedas. Um elástico segura todo o conjunto do foguete para um pouso seguro. O minifoguete utilizado (Figura 1(b)) para a obtenção dos dados experimentais foi o Epsilon-8. A ele foi acoplado um motor comercial de classificação E7 segundo a nomenclatura da NAR (National Association of Rocketry) [14]. Nesta nomenclatura a letra especifica que o impulso total produzido pelo motor está na faixa de 20,00 a 40,00 N.s, e o número em seguida especifica o empuxo médio produzido de 7 N. De fabricação do Grupo de Foguetes Carl Sagan (GFCS) da UFPR [8], o Epsilon-8 possui diâmetro máximo de 2,03 cm, comprimento de 50,1 cm e massa de 131,2 g [15]. Ele é o atual recordista brasileiro de apogeu para a classe E, atingindo em 2015 um apogeu de 723 m. Os dados deste voo foram obtidos com um altímetro micro-peak da Altus Metrum [9] e foram utilizados, neste trabalho, para os ajustes das equações de movimento em dados experimentais reais de v(t) e y(t). Os dados do voo foram obtidos em 08/08/2015 pelo Grupo de Foguetes Carl Sagan da Universidade Federal do Paraná (UFPR), que é liderado por Carlos Henrique Marchi, um dos autores do presente trabalho. foi utilizado o minifoguete Epsilon-8 com motor a propelente sólido classe E7. Parâmetros de importância sobre as características do motor classe E7 podem ser obtidos experimentalmente através de testes estáticos [2222. R.F.V. Marcus, Metodologia de Projeto e Validação de MotoresFoguete a Propelente Sólido. Dissertação de Mestrado, Universidade deSão Paulo, São Paulo (2013).]. Nestes testes, o motor é acoplado a um anteparo que ligado a um sistema eletrônico, fornece a força/empuxo em função do tempo [T(t)], quando o motor é ignitado. A Figura 3 representa um dos testes estáticos realizado, dentre 6 obtidos pelo GFCS da UFPR. O valor do empuxo médio, representado em vermelho, foi de 6,58 N.

Figura 3
Curva experimental de empuxo em função do tempo para o motor E7 utilizado no minifoguete.

No Gráfico da Figura 3, com os pontos experimentais representados por círculos, a magnitude do impulso total (Itot) é definida de acordo com a Equação (36) como a integral de [T(t)], durante o tempo de funcionamento do motor, ou seja, até o tempo final de queima (tq) do propelente:

(38) I t o t = 0 t q T ( t ) d t .

A magnitude do empuxo médio (Tm) que o motor exerce sobre o minifoguete é calculado pela Equação:

(39) T m = I t o t t q .

As equações de movimento descritas nas seções 2.3.1 2.3.1. Modelo I: Movimento sem a forçade arrasto do ar Quando a força de arrasto do ar é desprezada, a Equação (3) se torna: (4) m ⁢ d ⁢ v d ⁢ t = - v e ⁢ d ⁢ m d ⁢ t - m ⁢ g . Assumindo que dm/dt é uma constante nesta equação, então a massa diminui linearmente com o tempo: (5) m ⁢ ( t ) = m o - R ⁢ t , na qual R = −dm/dt (lembrando que dm/dt é negativo) e mo=mp + mf é a massa inicial do foguete, com mp denotando a massa total do propelente e mf a massa final do foguete [5, 6]. Substituindo a Equação (5) na (4), obtém-se: (6) d ⁢ v d ⁢ t = R ⁢ v e m o - R ⁢ t - g . As soluções da Equação (6) para v(t) e y(t), são dadas por: (7) v I ⁢ A ⁢ ( t ) = - g ⁢ t - v e ⁢ l ⁢ n ⁢ [ 1 - R m o ⁢ t ] , (8) y I ⁢ A ⁢ ( t ) = v e ⁢ t - 1 2 ⁢ g ⁢ t 2 + m o - R ⁢ t R ⁢ v e ⁢ l ⁢ n ⁢ [ 1 - R m o ⁢ t ] . Na etapa B, após a queima de todo o propelente (t≥tq), a Equação (3) se torna: (9) d ⁢ v d ⁢ t = - g , cujas soluções para v(t)e y(t) são, respectivamente: (10) v I ⁢ B = v q - g ⁢ ( t - t q ) , (11) y I ⁢ B = y q + v q ⁢ ( t - t q ) - g 2 ⁢ ( t - t q ) 2 , nas quais vq é a velocidade do foguete no tempo tq; yq é altura do foguete no tempo tq. De outra forma: vq = vIA(tq) e yq = yIA(tq). a 2.3.4 2.3.4. Modelo IV: Solução completa com forçade arrasto quadrática Uma solução analítica completa da Equação (24) para v(t), durante a queima do propelente, foi obtida por Rodrigues et al. (2014) [5]. O autor demonstrou, utilizando técnicas de cálculo diferencial, que a Equação (3) pode ser reescrita em termos de uma equação diferencial cujas soluções são funções de Bessel. Veja inicialmente que a Equação (24), agora com m = mo−Rt, pode ser escrita como: (31) d ⁢ v d ⁢ t + b m ⁢ v 2 - R ⁢ v e m + g = 0 , Durante as etapas para a solução da Equação (31), é possível realizar uma mudança da variável v, para uma nova variável χ [5]. Neste caso, esta equação se transforma em: (32) ξ 2 ⁢ d 2 ⁢ χ d ⁢ t 2 + ξ ⁢ d ⁢ χ d ⁢ t + ( ξ 2 - ν 2 ) ⁢ χ = 0 , que é uma Equação de Bessel de ordem ν [5, 19]. As constantes ξ e ν são definidas como: (33) ξ = ξ ⁢ ( m ) = 2 R ⁢ b ⁢ m ⁢ g , (34) ν = 2 ⁢ b ⁢ v e R . A solução geral da Equação (32) é uma combinação linear do tipo: (35) χ ⁢ ( ξ ) = C 1 ⁢ J ν ⁢ ( ξ ) + C 2 ⁢ Y ν ⁢ ( ξ ) , no qual Jν(ξ) e Yν(ξ) são funções de Bessel de primeira e segunda ordens respectivamente [5, 19, 20], C1 e C2 são constantes reais. Após algumas álgebras utilizando propriedades convenientes das funções de Bessel [5], é possível escrever a solução geral para a velocidade como: (36) v I ⁢ V ⁢ A ⁢ ( ξ ) = R 2 ⁢ b ⁢ [ ξ ⁢ C ⁢ Y ν + 1 ⁢ ( ξ ) + J ν + 1 ⁢ ( ξ ) C ⁢ Y ν ⁢ ( ξ ) + J ν ⁢ ( ξ ) - v o ] , (37) c ⁢ o ⁢ m ⁢ C = ( 2 ⁢ b ⁢ v o + R ⁢ ν ) ⁢ J ν ⁢ ( ξ o ) - R ⁢ ξ o ⁢ J ν + 1 ⁢ ( ξ o ) R ⁢ ξ o ⁢ Y ν + 1 ⁢ ( ξ o ) - ( 2 ⁢ b ⁢ v o + R ⁢ ν ) ⁢ Y ν ⁢ ( ξ ) na qual vo é a velocidade inicial do foguete, que pode ser definida como zero e ξo=2R⁢b⁢mo⁢g. Tendo em vista que ξ = ξ(m), Equação (36), e que m = m(t), Equação (5), é possível por meio da Equação (36) construir um gráfico da velocidade em função do tempo para (0≤t≤tq), etapa A, conhecendo-se os parâmetros mo,g,R,b e ve, previamente obtidos de cálculos teóricos/experimentais. A construção do gráfico y(t) pode ser obtida a partir de v(t) utilizando programas de computador adequados. Para este modelo, foi utilizado o programa GNU Octave [21] para simular as funções v(t) e y(t); a partir da velocidade foi obtida a altura através de um processo de integração numérica denominado quadratura de Gauss. No Apêndice I estão descritas as funções utilizadas no GNU Octave, bem como os scriptsimplementados para simulação do modelo apresentado nesta seção. No Apêndice II foi deixada uma explicação sucinta sobre a quadratura de Gauss. Para a etapa B neste modelo as Equações (28) e (29) podem ser aplicadas. Nesta condição deve-se considerar que m = m(tq). consideram um empuxo constante, igual ao valor médio do empuxo real calculado através da Equação (39) e representado pela linha vermelha na Figura 3. Esta aproximação é razoável no chamado regime estacionário e quando os regimes transientes (início e término do funcionamento do motor) são breves [2323. Richard Nakka’s Experimental Rocketry Web Site, disponível emhttps://www.nakka-rocketry.net/articles/nakka_theory_pages.pdf,acessado em 19/10/2020.
https://www.nakka-rocketry.net/articles/...
]. Assim, a magnitude do impulso total pode ser calculada como:

(40) I t o t = v e | d m d t | t q .

Outro parâmetro relevante é o impulso específico (Isp), obtido pela divisão do impulso total pelo peso do propelente (mpg) de acordo com a equação:

(41) I s p = I t o t m p g ,

na qual g é a magnitude da aceleração da gravidade ao nível do mar. Com esta grandeza mede-se a eficiência do propelente, já que, quanto maior for o impulso específico, maior será a variação da quantidade de movimento, por unidade de peso do propelente. Admite-se neste trabalho que o empuxo T(t) é constante, [44. A.L. Alves, A.N. Paneto, K.A. Littike, S.S. Bento e C.H. Marchi, Revista Brasileira de Ensino de Física 42, e20200390 (2020).]. Desta forma, considerando-se que ve e dm/dt são constantes, é possível estimar a velocidade de exaustão dos gases de acordo com a equação:

(42) v e = I t m p = I s p g .

Mede-se experimentalmente tq e mp sendo possível estimar R = dm/dt, como:

(43) R = m p / t q .

A Tabela 1 fornece a massa do propelente utilizado no minifoguete Epsilon-8 e os parâmetros médios relevantes obtidos dos testes estáticos do motor. Nesta tabela também se encontram outros parâmetros experimentais relacionados ao minifoguete e ao voo do mesmo. Outros parâmetros como velocidade máxima e altitude máxima atingida serão discutidos na análise gráfica da seção 3 3. Resultados Teóricos Com base nas Equações desenvolvidas na seção 2.3 foram construídos gráficos teóricos para os modelos de I a IV desta mesma seção. Foram utilizados os parâmetros descritos na Tabela 1 e aqueles obtidos a partir dela (ve e R). Nas Figuras 4(a)–(c) à esquerda estão apresentados gráficos de v(t) para os modelos de I a IV e, nas Figuras 4(d)–(f) à direita os respectivos gráficos de y(t). Os gráficos se distinguem por diferentes constantes de arrasto, b, que estão discriminados por diferentes cores. Quando a força de arrasto tiver uma dependência linear, a constante b será designada por bL e quando tiver uma dependência quadrática, será designada por bQ. Estas constantes se diferem pela unidade de medição: bL é dada em kg/s e bQ em kg/m. Figura 4 (a)–(c) Curvas de v(t) para os modelos de I a IV. O modelo I está representado pela curva em preto. Nas partes (d)–(f), à direita, os respectivos gráficos de y(t). Os valores destas constantes foram escolhidos para que as curvas de velocidade e altitude pudessem ser comparadas entre os diferentes modelos. Para efeitos de comparações, nas Figuras 4(a) e 4(d), na cor em preto, foram colocados os gráficos para v(t) e y(t) com b=0 relativo ao modelo I, juntamente com os gráficos referentes ao modelo II. Todos os gráficos foram plotados até seus respectivos tempos de apogeu, ta. Nos gráficos de v(t) este tempo se refere à velocidade nula para t > 0 de cada curva. Já nos gráficos de y(t), para os valores de t = ta, as curvas atingem a altitude máxima (apogeu). Todos os gráficos de velocidade são caracterizados por duas partes: a primeira nas quais as velocidades aumentam com o tempo, deste t = 0s ate t = tq = 5,358 s e a segunda nas quais as velocidades decrescem instantaneamente a partir deste tempo, tornando-se nulas em seus respectivos ta′⁢s. Na primeira parte do modelo I, por exemplo, as curvas v(t) e y(t) são descritas pelas Equações (7) e (8) para vIA(t) e yIA(t), respectivamente. Após este tempo as curvas são descritas pelas Equações (10) e (11) para vIB(t) e yIB(t), respectivamente. A mesma sistemática é válida para as demais curvas. Os picos nos gráficos de v(t) podem ser explicados como se segue. Antes de tq o termo de empuxo T=-ve⁢d⁢md⁢t supera os da soma (mg + |Fr|) na Equação (3), o que confere uma aceleração positiva ao foguete. Teoricamente após este tempo o termo T torna-se instantaneamente nulo e apenas as forças retardatárias mg e Fr prevalecem. Isto confere uma mudança rápida no sentido da aceleração do foguete e assim a velocidade passa a decrescer instantaneamente. Vale lembrar que em tq a força Fr atinge seu valor máximo (pois a velocidade é máxima) o que pode contribuir ainda mais com os picos característicos nos gráficos de v(t). Nos gráficos y(t) surge um ponto de inflexão em t = tq no qual a velocidade de um foguete passa a decrescer. Observa-se nos gráficos da Figura 1(a) e 1(d) que, como esperado, os valores de velocidades e altitudes para b=0 (modelo I) são superiores ao de qualquer um dos modelos que apresentam b≠ 0. Ainda, como esperado para os modelos de II a IV, um aumento nos valores de b acarreta em decréscimos nas curvas de v(t) e y(t). As curvas com bL = 0,001 kg/s (modelo II) e bQ = 0,00005 kg/m (modelo III), que apresentam menores valores desta grandeza, são as que mais se aproximam da curva do modelo I. A velocidade para o modelo I em tq é da ordem de 271 m/s e supera a do modelo II com bL = 0,001 kg/s em aproximadamente 6 m/s. Após este tempo, a ausência do arrasto no modelo I confere o maior valor de ta e assim maior altitude (≈4.463,2 m) excedendo a do modelo II com bL = 0,001 kg/s em aproximadamente 801,5 m. Ainda, as curvas de v(t) para este modelo, para tempos até tq, apresentam diferentes formas conforme os valores de b aumentam. Para valores de bL ≈ 0,001 kg/s as derivadas em cada ponto da curva aumentam com o tempo, indicando um aumento da aceleração do foguete com o tempo. Conforme os valores de bL se tornam maiores as curvas mudam de tal forma que as derivadas em cada ponto passam a diminuir, indicando um decréscimo da aceleração do foguete com o tempo. Estes são os casos das curvas para bL = 0,01, 0,03 e 0,05 kg/s. Para que os gráficos de velocidade dos modelos III (Figura 4(b)) e IV (Figura 4(c)) estejam na faixa de 0–300 m/s (para se realizar uma comparação com os modelos I e II) é necessário que os valores de bQ sejam inferiores a uma potência de pelo menos 10−2. Por representarem modelos mais realísticos para movimento de foguetes, estes modelos puderam ser simulados com a mesma faixa de valores da constante bQ. Inicialmente, aparenta-se não haver distinção entre as formas das curvas, considerando um mesmo valor de b dentre esses modelos. Observa-se, entretanto, que os valores de velocidade nos tempos tq e ta são ligeiramente diferentes, o que acarreta em apogeus distintos para estes modelos (Figuras 4(e) e (f)). Para efeitos de comparação, são apresentadas na Figura 5(a) as curvas v(t) para os modelos III e IV com os valores de bQ = 0,00005 kg/m e 0,0002 kg/m, e na Figura 5(b), as respectivas curvas y(t). Observa-se que, apesar da semelhança marcante entre as curvas (para ambos os valores de bQ), os valores de velocidade para o modelo III são superiores aos do modelo IV em todos os instantes de tempo e, como resultado disto, as altitudes para este modelo são superiores às do modelo IV. A mesma análise é válida para as curvas com bQ = 0,0001 kg/m e 0,0005 kg/m. Figura 5 Comparação para as curvas (a) v(t) e (b) y(t) para os modelos III e IV. Foram utilizados valores de bQ = 0,00005 kg/m e 0,0002 kg/m. Os modelos II, III e IV supracitados parecem todos serem úteis para descrever o movimento vertical ascendente de foguetes. Entretanto, sabe-se que o modelo II não se adequa ao movimento de foguete, pois este modelo se aplica a objetos que se deslocam em meios viscosos e adaptável para situações nas quais as velocidades são inferiores a 24 m/s [18]. O minifoguete Epsilon-8 atinge uma velocidade ≈128,3 m/s até o instante tq e como será visto, o modelo I não é muito ideal para simular o movimento deste minifoguete. São feitas, na próxima seção, análises comparativas com um voo real (experimental) do minifoguete Epsilon-8, abordado nas seções 2.1 e 2.4. .

Tabela 1
Parâmetros do motor do minifoguete Epsilon-8. A grandeza mp é a massa do propelente, It é o valor médio do impulso total, obtido do teste estático de 6 motores idênticos e tq é o tempo de queima do propelente; moé a massa inicial do minifoguete (antes da queima do propelente), mfé a massa final do mesmo e taé o tempo necessário para atingir o apogeu.

Utilizando a Equação (42) e os valores de It (médio) e mp da Tabela 1 é possível estimar ve = (780,8 ± 0,4) m/s. Adicionalmente, com os valores de mp e tq da Tabela 1 é possível estimar, utilizando a Equação (43), R = (8,343 ± 0,003) g/s.

3. Resultados Teóricos

Com base nas Equações desenvolvidas na seção 2.3 2.3. Soluções Analíticas para a Equaçãode Movimento de um Foguete Nas seções 2.3.1 a 2.3.4 apresentam-se as soluções para a equação de movimento de um foguete durante o tempo de queima (tq) do propelente, fase propulsada, que aqui denomina-se de etapa A. Após este tempo, na fase balística, denominam-se as soluções como etapa B. Nesta Etapa, após o tempo tq, o motor deixa de funcionar e a magnitude do empuxo T, na Equação (3), se torna nula. Nesta etapa, o foguete segue um voo no qual apenas as forças de resistência do ar e o peso estão presentes. Utiliza-se a notação W(t)nA,B para designar uma equação de trajetória (soluções da Equação (3)). Nesta notação, W(t) representa v(t) ou y(t), n representa o modelo (I, II …) e A ou B, uma das etapas de voo. No modelo I a seguir, coloca-se a solução para o movimento do foguete, desconsiderando a força de arrasto do ar, com a intenção de que o estudante aprecie as soluções e compare com modelos que levam em consideração esta força. Este movimento representaria o movimento de um foguete lançado da superfície de um astro sem atmosfera, como a Lua. Os estudantes são encorajados a verificar as deduções e soluções nas referências deixadas, para cada uma das equações de trajetória aqui apresentadas. Para todas as equações aqui colocadas, considerou-se que as velocidades e alturas iniciais (t = 0s) possuem, ambas, valores nulos. 2.3.1. Modelo I: Movimento sem a forçade arrasto do ar Quando a força de arrasto do ar é desprezada, a Equação (3) se torna: (4) m ⁢ d ⁢ v d ⁢ t = - v e ⁢ d ⁢ m d ⁢ t - m ⁢ g . Assumindo que dm/dt é uma constante nesta equação, então a massa diminui linearmente com o tempo: (5) m ⁢ ( t ) = m o - R ⁢ t , na qual R = −dm/dt (lembrando que dm/dt é negativo) e mo=mp + mf é a massa inicial do foguete, com mp denotando a massa total do propelente e mf a massa final do foguete [5, 6]. Substituindo a Equação (5) na (4), obtém-se: (6) d ⁢ v d ⁢ t = R ⁢ v e m o - R ⁢ t - g . As soluções da Equação (6) para v(t) e y(t), são dadas por: (7) v I ⁢ A ⁢ ( t ) = - g ⁢ t - v e ⁢ l ⁢ n ⁢ [ 1 - R m o ⁢ t ] , (8) y I ⁢ A ⁢ ( t ) = v e ⁢ t - 1 2 ⁢ g ⁢ t 2 + m o - R ⁢ t R ⁢ v e ⁢ l ⁢ n ⁢ [ 1 - R m o ⁢ t ] . Na etapa B, após a queima de todo o propelente (t≥tq), a Equação (3) se torna: (9) d ⁢ v d ⁢ t = - g , cujas soluções para v(t)e y(t) são, respectivamente: (10) v I ⁢ B = v q - g ⁢ ( t - t q ) , (11) y I ⁢ B = y q + v q ⁢ ( t - t q ) - g 2 ⁢ ( t - t q ) 2 , nas quais vq é a velocidade do foguete no tempo tq; yq é altura do foguete no tempo tq. De outra forma: vq = vIA(tq) e yq = yIA(tq). 2.3.2. Modelo II: Movimento com força de arrasto linear A força de arrasto do ar depende da velocidade v de um objeto e geralmente é descrita em termos de potência n da velocidade, ou seja, (12) F r = - b ⁢ v n , na qual b é uma constante positiva que especifica o poder da força de arrasto [18]. Existem dois regimes diferentes no qual essa dependência pode ser do tipo linear ou quadrática em v.O parâmetro que diferencia estes regimes é o número de Reynolds (nRe) [4]. Experimentalmente, para um objeto relativamente pequeno movendo-se no ar, n=1 na Equação (12) para velocidades de até 24 m/s. Para velocidades mais elevadas, a força de arrasto é aproximadamente proporcional ao quadrado da velocidade [18]. As constantes b’s na Equação (12) se diferem por ordem de grandeza e unidade de medição. Nesta seção frisa-se mais na obtenção das equações de movimento e, por simplicidade, utilizou-se a mesma letra, sem subíndice, para representar esta variável. Na seção 3, faz-se melhor a distinção das constantes b’s dos modelos II, III e IV. Assumindo para o modelo II uma dependência linear para Fr, a Equação (3) se torna: (13) m ⁢ d ⁢ v d ⁢ t = - v e ⁢ d ⁢ m d ⁢ t - m ⁢ g - b ⁢ v . Substituindo a Equação (5) na Equação (13) e reagrupando os termos, obtém-se: (14) d ⁢ v d ⁢ t + b m o - R ⁢ t ⁢ v = v e m o - R ⁢ t ⁢ R - g . As soluções analíticas para v(t) e y(t) são dadas pelas Equações (15) e (16) [6]: (15) v I ⁢ I ⁢ A ⁢ ( t ) = v e ⁢ R b + g R - b ⁢ ( m o - R ⁢ t ) + c 1 ⁢ ( m o - R ⁢ t ) b R , em que: (16) c 1 = - v e ⁢ R ⁢ ( R - b ) - g ⁢ m o ⁢ b R ⁢ ( R - b ) ⁢ m o b R ; (17) y I ⁢ I ⁢ A ⁢ ( t ) = v e ⁢ R b ⁢ t - g 2 ⁢ R ⁢ ( R - b ) ⁢ ( m o - R ⁢ t ) 2 - c 1 R + b ⁢ ( m o - R ⁢ t ) 1 + b R + c 2 , em que: (18) c 2 = g ⁢ m o 2 2 ⁢ R ⁢ ( R - b ) + c 1 ⁢ m o 1 + b R R + b . Após este tempo, na etapa B, a Equação (3) se torna: (19) m f ⁢ d ⁢ v d ⁢ t = - m f ⁢ g - b ⁢ v , na qual mf é a massa final do foguete após a queima de todo o propelente no instante tq. As soluções da Equação (19) para v(t) e y(t) podem ser resolvidas por uma integração direta. As soluções são, respectivamente, (20) v I ⁢ I ⁢ B ⁢ ( t ) = [ ( m f ⁢ g b + v q ) ⁢ e - b m f ⁢ ( t - t q ) - m f ⁢ g b ] (21) y I ⁢ I ⁢ B ⁢ ( t ) = y q - m f b ⁢ g ⁢ ( t - t q ) + m f 2 b 2 ⁢ g ⁢ ( 1 + b m f ⁢ g ⁢ v q ) y I ⁢ I ⁢ B ⁢ ( t ) × ( 1 - e - b m f ⁢ ( t - t q ) ) , nas quais vq = vIIA(tq) e yq = yIIA(tq). 2.3.3. Modelo III: Solução aproximada commassa constante e força de arrastoquadrática Um modelo mais realístico do movimento de um foguete, durante sua subida, considera o termo n=2 na Equação (12). O Epsilon-8, por exemplo, atinge velocidades ≈100 m/s em um intervalo de tempo relativamente curto ≈3 s. Neste caso, nRe≥103 e a dependência quadrática deve ser levada em consideração [4]. Como aproximação considera-se um valor médio para a massa do foguete (m), constante durante a queima do propelente [4, 16, 17]. Assim, adota-se também m constante durante todo o voo do foguete. Esta aproximação é válida se o tempo de queima é pequeno quando comparado ao tempo para o foguete atingir o apogeu. Neste caso, m é dado pela média: (22) m = m o + m f 2 , na qual mo e mf são as respectivas massa inicial e final do minifoguete. A equação (3) escrita na forma: (23) m ⁢ d ⁢ v d ⁢ t = - v e ⁢ d ⁢ m d ⁢ t - m ⁢ g - b ⁢ v 2 , pode ser reescrita como: (24) m β - b ⁢ v 2 ⁢ d ⁢ v = d ⁢ t , na qual: β = -ve⁢d⁢md⁢t-m⁢g. As soluções da Equação (24) para v(t) e y(t), na etapa A, descritas por Waters (2014) [17] e relatadas por Alves et al. [4], são dadas respectivamente por: (25) v I ⁢ I ⁢ I ⁢ A ⁢ ( t ) = β b ⁢ t ⁢ g ⁢ h ⁢ [ ( b ⁢ β m ) ⁢ t ] , (26) y I ⁢ I ⁢ I ⁢ A ⁢ ( t ) = m b ⁢ l ⁢ n ⁢ [ cosh ⁢ ( b ⁢ β m ⁢ t ) ] . As Equações (25) e (26) são válidas para β> 0. Esta condição é válida para o Epsilon-8 (veja a seção 2.4). Para a etapa B, o empuxo é nulo. Logo, a Equação (23) torna-se: (27) m ⁢ d ⁢ v d ⁢ t = - m ⁢ g - b ⁢ v 2 , cujas soluções para v(t) e y(t) são, respectivamente: (28) v I ⁢ I ⁢ I ⁢ B ⁢ ( t ) = m ⁢ g b ⁢ t ⁢ g ⁢ [ b ⁢ g m ⁢ ( t q - t ) + t ⁢ g - 1 ⁢ ( b m ⁢ g ⁢ v q ) ] , (29) y I ⁢ I ⁢ I ⁢ B ⁢ ( t ) = y q + m b ⁢ l ⁢ n ⁢ { c ⁢ o ⁢ s ⁢ [ b ⁢ g m ⁢ ( t a - t ) ] c ⁢ o ⁢ s ⁢ [ b ⁢ g m ⁢ ( t a - t q ) ] } . Na Equação (28) tem-se que vq = vIIIA(tq) e, na Equação (29), yq = yIIIA(tq), em que ta é o instante no qual o foguete atinge o apogeu. Este instante é dado pela Equação: (30) t a = t q + m b ⁢ g ⁢ t ⁢ g - 1 ⁢ ( b m ⁢ g ⁢ v q ) . 2.3.4. Modelo IV: Solução completa com forçade arrasto quadrática Uma solução analítica completa da Equação (24) para v(t), durante a queima do propelente, foi obtida por Rodrigues et al. (2014) [5]. O autor demonstrou, utilizando técnicas de cálculo diferencial, que a Equação (3) pode ser reescrita em termos de uma equação diferencial cujas soluções são funções de Bessel. Veja inicialmente que a Equação (24), agora com m = mo−Rt, pode ser escrita como: (31) d ⁢ v d ⁢ t + b m ⁢ v 2 - R ⁢ v e m + g = 0 , Durante as etapas para a solução da Equação (31), é possível realizar uma mudança da variável v, para uma nova variável χ [5]. Neste caso, esta equação se transforma em: (32) ξ 2 ⁢ d 2 ⁢ χ d ⁢ t 2 + ξ ⁢ d ⁢ χ d ⁢ t + ( ξ 2 - ν 2 ) ⁢ χ = 0 , que é uma Equação de Bessel de ordem ν [5, 19]. As constantes ξ e ν são definidas como: (33) ξ = ξ ⁢ ( m ) = 2 R ⁢ b ⁢ m ⁢ g , (34) ν = 2 ⁢ b ⁢ v e R . A solução geral da Equação (32) é uma combinação linear do tipo: (35) χ ⁢ ( ξ ) = C 1 ⁢ J ν ⁢ ( ξ ) + C 2 ⁢ Y ν ⁢ ( ξ ) , no qual Jν(ξ) e Yν(ξ) são funções de Bessel de primeira e segunda ordens respectivamente [5, 19, 20], C1 e C2 são constantes reais. Após algumas álgebras utilizando propriedades convenientes das funções de Bessel [5], é possível escrever a solução geral para a velocidade como: (36) v I ⁢ V ⁢ A ⁢ ( ξ ) = R 2 ⁢ b ⁢ [ ξ ⁢ C ⁢ Y ν + 1 ⁢ ( ξ ) + J ν + 1 ⁢ ( ξ ) C ⁢ Y ν ⁢ ( ξ ) + J ν ⁢ ( ξ ) - v o ] , (37) c ⁢ o ⁢ m ⁢ C = ( 2 ⁢ b ⁢ v o + R ⁢ ν ) ⁢ J ν ⁢ ( ξ o ) - R ⁢ ξ o ⁢ J ν + 1 ⁢ ( ξ o ) R ⁢ ξ o ⁢ Y ν + 1 ⁢ ( ξ o ) - ( 2 ⁢ b ⁢ v o + R ⁢ ν ) ⁢ Y ν ⁢ ( ξ ) na qual vo é a velocidade inicial do foguete, que pode ser definida como zero e ξo=2R⁢b⁢mo⁢g. Tendo em vista que ξ = ξ(m), Equação (36), e que m = m(t), Equação (5), é possível por meio da Equação (36) construir um gráfico da velocidade em função do tempo para (0≤t≤tq), etapa A, conhecendo-se os parâmetros mo,g,R,b e ve, previamente obtidos de cálculos teóricos/experimentais. A construção do gráfico y(t) pode ser obtida a partir de v(t) utilizando programas de computador adequados. Para este modelo, foi utilizado o programa GNU Octave [21] para simular as funções v(t) e y(t); a partir da velocidade foi obtida a altura através de um processo de integração numérica denominado quadratura de Gauss. No Apêndice I estão descritas as funções utilizadas no GNU Octave, bem como os scriptsimplementados para simulação do modelo apresentado nesta seção. No Apêndice II foi deixada uma explicação sucinta sobre a quadratura de Gauss. Para a etapa B neste modelo as Equações (28) e (29) podem ser aplicadas. Nesta condição deve-se considerar que m = m(tq). foram construídos gráficos teóricos para os modelos de I a IV desta mesma seção. Foram utilizados os parâmetros descritos na Tabela 1 e aqueles obtidos a partir dela (ve e R). Nas Figuras 4(a)–(c) à esquerda estão apresentados gráficos de v(t) para os modelos de I a IV e, nas Figuras 4(d)–(f) à direita os respectivos gráficos de y(t). Os gráficos se distinguem por diferentes constantes de arrasto, b, que estão discriminados por diferentes cores. Quando a força de arrasto tiver uma dependência linear, a constante b será designada por bL e quando tiver uma dependência quadrática, será designada por bQ. Estas constantes se diferem pela unidade de medição: bL é dada em kg/s e bQ em kg/m.

Figura 4
(a)–(c) Curvas de v(t) para os modelos de I a IV. O modelo I está representado pela curva em preto. Nas partes (d)–(f), à direita, os respectivos gráficos de y(t).

Os valores destas constantes foram escolhidos para que as curvas de velocidade e altitude pudessem ser comparadas entre os diferentes modelos. Para efeitos de comparações, nas Figuras 4(a) e 4(d), na cor em preto, foram colocados os gráficos para v(t) e y(t) com b=0 relativo ao modelo I, juntamente com os gráficos referentes ao modelo II. Todos os gráficos foram plotados até seus respectivos tempos de apogeu, ta. Nos gráficos de v(t) este tempo se refere à velocidade nula para t > 0 de cada curva. Já nos gráficos de y(t), para os valores de t = ta, as curvas atingem a altitude máxima (apogeu).

Todos os gráficos de velocidade são caracterizados por duas partes: a primeira nas quais as velocidades aumentam com o tempo, deste t = 0s ate t = tq = 5,358 s e a segunda nas quais as velocidades decrescem instantaneamente a partir deste tempo, tornando-se nulas em seus respectivos tas. Na primeira parte do modelo I, por exemplo, as curvas v(t) e y(t) são descritas pelas Equações (7) e (8) para vIA(t) e yIA(t), respectivamente. Após este tempo as curvas são descritas pelas Equações (10) e (11) para vIB(t) e yIB(t), respectivamente. A mesma sistemática é válida para as demais curvas.

Os picos nos gráficos de v(t) podem ser explicados como se segue. Antes de tq o termo de empuxo T=-vedmdt supera os da soma (mg + |Fr|) na Equação (3), o que confere uma aceleração positiva ao foguete. Teoricamente após este tempo o termo T torna-se instantaneamente nulo e apenas as forças retardatárias mg e Fr prevalecem. Isto confere uma mudança rápida no sentido da aceleração do foguete e assim a velocidade passa a decrescer instantaneamente. Vale lembrar que em tq a força Fr atinge seu valor máximo (pois a velocidade é máxima) o que pode contribuir ainda mais com os picos característicos nos gráficos de v(t). Nos gráficos y(t) surge um ponto de inflexão em t = tq no qual a velocidade de um foguete passa a decrescer.

Observa-se nos gráficos da Figura 1(a) e 1(d) que, como esperado, os valores de velocidades e altitudes para b=0 (modelo I) são superiores ao de qualquer um dos modelos que apresentam b≠ 0. Ainda, como esperado para os modelos de II a IV, um aumento nos valores de b acarreta em decréscimos nas curvas de v(t) e y(t). As curvas com bL = 0,001 kg/s (modelo II) e bQ = 0,00005 kg/m (modelo III), que apresentam menores valores desta grandeza, são as que mais se aproximam da curva do modelo I. A velocidade para o modelo I em tq é da ordem de 271 m/s e supera a do modelo II com bL = 0,001 kg/s em aproximadamente 6 m/s. Após este tempo, a ausência do arrasto no modelo I confere o maior valor de ta e assim maior altitude (≈4.463,2 m) excedendo a do modelo II com bL = 0,001 kg/s em aproximadamente 801,5 m. Ainda, as curvas de v(t) para este modelo, para tempos até tq, apresentam diferentes formas conforme os valores de b aumentam. Para valores de bL ≈ 0,001 kg/s as derivadas em cada ponto da curva aumentam com o tempo, indicando um aumento da aceleração do foguete com o tempo. Conforme os valores de bL se tornam maiores as curvas mudam de tal forma que as derivadas em cada ponto passam a diminuir, indicando um decréscimo da aceleração do foguete com o tempo. Estes são os casos das curvas para bL = 0,01, 0,03 e 0,05 kg/s.

Para que os gráficos de velocidade dos modelos III (Figura 4(b)) e IV (Figura 4(c)) estejam na faixa de 0–300 m/s (para se realizar uma comparação com os modelos I e II) é necessário que os valores de bQ sejam inferiores a uma potência de pelo menos 10−2. Por representarem modelos mais realísticos para movimento de foguetes, estes modelos puderam ser simulados com a mesma faixa de valores da constante bQ. Inicialmente, aparenta-se não haver distinção entre as formas das curvas, considerando um mesmo valor de b dentre esses modelos. Observa-se, entretanto, que os valores de velocidade nos tempos tq e ta são ligeiramente diferentes, o que acarreta em apogeus distintos para estes modelos (Figuras 4(e) e (f)).

Para efeitos de comparação, são apresentadas na Figura 5(a) as curvas v(t) para os modelos III e IV com os valores de bQ = 0,00005 kg/m e 0,0002 kg/m, e na Figura 5(b), as respectivas curvas y(t). Observa-se que, apesar da semelhança marcante entre as curvas (para ambos os valores de bQ), os valores de velocidade para o modelo III são superiores aos do modelo IV em todos os instantes de tempo e, como resultado disto, as altitudes para este modelo são superiores às do modelo IV. A mesma análise é válida para as curvas com bQ = 0,0001 kg/m e 0,0005 kg/m.

Figura 5
Comparação para as curvas (a) v(t) e (b) y(t) para os modelos III e IV. Foram utilizados valores de bQ = 0,00005 kg/m e 0,0002 kg/m.

Os modelos II, III e IV supracitados parecem todos serem úteis para descrever o movimento vertical ascendente de foguetes. Entretanto, sabe-se que o modelo II não se adequa ao movimento de foguete, pois este modelo se aplica a objetos que se deslocam em meios viscosos e adaptável para situações nas quais as velocidades são inferiores a 24 m/s [1818. Mecânica Newtoniana: Partícula Única,disponível em:http://jararaca.ufsm.br/websites/petfisica/download/Arquivos/Capitulo%202.pdf, acessado em 19/10/2020.
http://jararaca.ufsm.br/websites/petfisi...
]. O minifoguete Epsilon-8 atinge uma velocidade ≈128,3 m/s até o instante tq e como será visto, o modelo I não é muito ideal para simular o movimento deste minifoguete. São feitas, na próxima seção, análises comparativas com um voo real (experimental) do minifoguete Epsilon-8, abordado nas seções 2.1 2.1. O minifoguete Um minifoguete, assim como um foguete real, é um veículo que se desloca expelindo, em sua parte traseira, um fluxo de gás à alta velocidade devido à queima do propelente [11]. De acordo com a terceira lei de Newton, esta ação resulta em uma reação de mesma intensidade, porém sentido oposto, deslocando o foguete para cima a partir do solo [4]. Minifoguetes são utilizados no ensino e em competições. Já outros, mais sofisticados do tipo sondagem, são aplicados às pesquisas aeroespaciais cujas altitudes são restritas à troposfera [12]. Na Figura 1(a) está ilustrada a seção transversal de um minifoguete, utilizado em atividades de ensino, com as principais partes que o compõe. Na Figura 1(b) é mostrada uma fotografia do minifoguete utilizado neste trabalho. Figura 1 (a) Ilustração do corte de seção transversal de um minifoguete. (b) Minifoguete Epsilon-8. Como se observa na Figura 1(a) a estrutura do minifoguete é delimitada pela linha em preto e inclui a coifa (nariz), a parte traseira que sustenta as empenas (ou aletas) e o motor demarcado pela linha azul. O motor contém o propelente na cor amarela que está no interior de um cilindro por onde é feito um pequeno furo em sua extremidade direita [13]. A queima do propelente é feita por um dispositivo de ignição elétrica que contém um ignitor e dois condutores. Estes ao serem ligados a uma pequena bateria transmitem energia que inflama o ignitor e consequentemente o propelente. Um dispositivo de ejeção do motor (não ilustrado aqui) desacopla o conjunto constituído pela coifa e compartimento do altímetro, o que automaticamente aciona o paraquedas. Um elástico segura todo o conjunto do foguete para um pouso seguro. O minifoguete utilizado (Figura 1(b)) para a obtenção dos dados experimentais foi o Epsilon-8. A ele foi acoplado um motor comercial de classificação E7 segundo a nomenclatura da NAR (National Association of Rocketry) [14]. Nesta nomenclatura a letra especifica que o impulso total produzido pelo motor está na faixa de 20,00 a 40,00 N.s, e o número em seguida especifica o empuxo médio produzido de 7 N. De fabricação do Grupo de Foguetes Carl Sagan (GFCS) da UFPR [8], o Epsilon-8 possui diâmetro máximo de 2,03 cm, comprimento de 50,1 cm e massa de 131,2 g [15]. Ele é o atual recordista brasileiro de apogeu para a classe E, atingindo em 2015 um apogeu de 723 m. Os dados deste voo foram obtidos com um altímetro micro-peak da Altus Metrum [9] e foram utilizados, neste trabalho, para os ajustes das equações de movimento em dados experimentais reais de v(t) e y(t). Os dados do voo foram obtidos em 08/08/2015 pelo Grupo de Foguetes Carl Sagan da Universidade Federal do Paraná (UFPR), que é liderado por Carlos Henrique Marchi, um dos autores do presente trabalho. e 2.4 2.4. Parâmetros teóricos/experimentaisrelacionados ao minifoguete Epsilon-8 Como descrito na seção 2.1 foi utilizado o minifoguete Epsilon-8 com motor a propelente sólido classe E7. Parâmetros de importância sobre as características do motor classe E7 podem ser obtidos experimentalmente através de testes estáticos [22]. Nestes testes, o motor é acoplado a um anteparo que ligado a um sistema eletrônico, fornece a força/empuxo em função do tempo [T(t)], quando o motor é ignitado. A Figura 3 representa um dos testes estáticos realizado, dentre 6 obtidos pelo GFCS da UFPR. O valor do empuxo médio, representado em vermelho, foi de 6,58 N. Figura 3 Curva experimental de empuxo em função do tempo para o motor E7 utilizado no minifoguete. No Gráfico da Figura 3, com os pontos experimentais representados por círculos, a magnitude do impulso total (Itot) é definida de acordo com a Equação (36) como a integral de [T(t)], durante o tempo de funcionamento do motor, ou seja, até o tempo final de queima (tq) do propelente: (38) I t ⁢ o ⁢ t = ∫ 0 t q T ⁢ ( t ) ⁢ d t . A magnitude do empuxo médio (Tm) que o motor exerce sobre o minifoguete é calculado pela Equação: (39) T m = I t ⁢ o ⁢ t t q . As equações de movimento descritas nas seções 2.3.1 a 2.3.4 consideram um empuxo constante, igual ao valor médio do empuxo real calculado através da Equação (39) e representado pela linha vermelha na Figura 3. Esta aproximação é razoável no chamado regime estacionário e quando os regimes transientes (início e término do funcionamento do motor) são breves [23]. Assim, a magnitude do impulso total pode ser calculada como: (40) I t ⁢ o ⁢ t = v e ⁢ | d ⁢ m d ⁢ t | ⁢ t q . Outro parâmetro relevante é o impulso específico (Isp), obtido pela divisão do impulso total pelo peso do propelente (mpg) de acordo com a equação: (41) I s ⁢ p = I t ⁢ o ⁢ t m p ⁢ g , na qual g é a magnitude da aceleração da gravidade ao nível do mar. Com esta grandeza mede-se a eficiência do propelente, já que, quanto maior for o impulso específico, maior será a variação da quantidade de movimento, por unidade de peso do propelente. Admite-se neste trabalho que o empuxo T(t) é constante, [4]. Desta forma, considerando-se que ve e dm/dt são constantes, é possível estimar a velocidade de exaustão dos gases de acordo com a equação: (42) v e = I t m p = I s ⁢ p ⁢ g . Mede-se experimentalmente tq e mp sendo possível estimar R = dm/dt, como: (43) R = m p / t q . A Tabela 1 fornece a massa do propelente utilizado no minifoguete Epsilon-8 e os parâmetros médios relevantes obtidos dos testes estáticos do motor. Nesta tabela também se encontram outros parâmetros experimentais relacionados ao minifoguete e ao voo do mesmo. Outros parâmetros como velocidade máxima e altitude máxima atingida serão discutidos na análise gráfica da seção 3. Tabela 1 Parâmetros do motor do minifoguete Epsilon-8. A grandeza mp é a massa do propelente, It é o valor médio do impulso total, obtido do teste estático de 6 motores idênticos e tq é o tempo de queima do propelente; moé a massa inicial do minifoguete (antes da queima do propelente), mfé a massa final do mesmo e taé o tempo necessário para atingir o apogeu. mp(g) It(médio)(N.s) tq (s) mo (g) mf(g) ta(s) 44,70 ± 0,01 34,90 ± 0,01 5,358 ± 0,001 131,20 ± 0,01 86,50 ± 0,01 11,7 ± 0,1 Utilizando a Equação (42) e os valores de It (médio) e mp da Tabela 1 é possível estimar ve = (780,8 ± 0,4) m/s. Adicionalmente, com os valores de mp e tq da Tabela 1 é possível estimar, utilizando a Equação (43), R = (8,343 ± 0,003) g/s. .

4. Aplicação Real: Simulações dos gráficos v(t) e y(t) para o voo do minifoguete Epsilon-8

Nas Figuras 6(a)–(c) são apresentados pontos experimentais (círculos abertos) e as curvas simuladas (linhas contínuas) para v(t) e y(t) do voo do minifoguete Epsilon-8. Os pontos experimentais, obtidas com o uso do altímetro micro-peak, aparecem em todos os gráficos. Todas as curvas, assim como as da seção 3 3. Resultados Teóricos Com base nas Equações desenvolvidas na seção 2.3 foram construídos gráficos teóricos para os modelos de I a IV desta mesma seção. Foram utilizados os parâmetros descritos na Tabela 1 e aqueles obtidos a partir dela (ve e R). Nas Figuras 4(a)–(c) à esquerda estão apresentados gráficos de v(t) para os modelos de I a IV e, nas Figuras 4(d)–(f) à direita os respectivos gráficos de y(t). Os gráficos se distinguem por diferentes constantes de arrasto, b, que estão discriminados por diferentes cores. Quando a força de arrasto tiver uma dependência linear, a constante b será designada por bL e quando tiver uma dependência quadrática, será designada por bQ. Estas constantes se diferem pela unidade de medição: bL é dada em kg/s e bQ em kg/m. Figura 4 (a)–(c) Curvas de v(t) para os modelos de I a IV. O modelo I está representado pela curva em preto. Nas partes (d)–(f), à direita, os respectivos gráficos de y(t). Os valores destas constantes foram escolhidos para que as curvas de velocidade e altitude pudessem ser comparadas entre os diferentes modelos. Para efeitos de comparações, nas Figuras 4(a) e 4(d), na cor em preto, foram colocados os gráficos para v(t) e y(t) com b=0 relativo ao modelo I, juntamente com os gráficos referentes ao modelo II. Todos os gráficos foram plotados até seus respectivos tempos de apogeu, ta. Nos gráficos de v(t) este tempo se refere à velocidade nula para t > 0 de cada curva. Já nos gráficos de y(t), para os valores de t = ta, as curvas atingem a altitude máxima (apogeu). Todos os gráficos de velocidade são caracterizados por duas partes: a primeira nas quais as velocidades aumentam com o tempo, deste t = 0s ate t = tq = 5,358 s e a segunda nas quais as velocidades decrescem instantaneamente a partir deste tempo, tornando-se nulas em seus respectivos ta′⁢s. Na primeira parte do modelo I, por exemplo, as curvas v(t) e y(t) são descritas pelas Equações (7) e (8) para vIA(t) e yIA(t), respectivamente. Após este tempo as curvas são descritas pelas Equações (10) e (11) para vIB(t) e yIB(t), respectivamente. A mesma sistemática é válida para as demais curvas. Os picos nos gráficos de v(t) podem ser explicados como se segue. Antes de tq o termo de empuxo T=-ve⁢d⁢md⁢t supera os da soma (mg + |Fr|) na Equação (3), o que confere uma aceleração positiva ao foguete. Teoricamente após este tempo o termo T torna-se instantaneamente nulo e apenas as forças retardatárias mg e Fr prevalecem. Isto confere uma mudança rápida no sentido da aceleração do foguete e assim a velocidade passa a decrescer instantaneamente. Vale lembrar que em tq a força Fr atinge seu valor máximo (pois a velocidade é máxima) o que pode contribuir ainda mais com os picos característicos nos gráficos de v(t). Nos gráficos y(t) surge um ponto de inflexão em t = tq no qual a velocidade de um foguete passa a decrescer. Observa-se nos gráficos da Figura 1(a) e 1(d) que, como esperado, os valores de velocidades e altitudes para b=0 (modelo I) são superiores ao de qualquer um dos modelos que apresentam b≠ 0. Ainda, como esperado para os modelos de II a IV, um aumento nos valores de b acarreta em decréscimos nas curvas de v(t) e y(t). As curvas com bL = 0,001 kg/s (modelo II) e bQ = 0,00005 kg/m (modelo III), que apresentam menores valores desta grandeza, são as que mais se aproximam da curva do modelo I. A velocidade para o modelo I em tq é da ordem de 271 m/s e supera a do modelo II com bL = 0,001 kg/s em aproximadamente 6 m/s. Após este tempo, a ausência do arrasto no modelo I confere o maior valor de ta e assim maior altitude (≈4.463,2 m) excedendo a do modelo II com bL = 0,001 kg/s em aproximadamente 801,5 m. Ainda, as curvas de v(t) para este modelo, para tempos até tq, apresentam diferentes formas conforme os valores de b aumentam. Para valores de bL ≈ 0,001 kg/s as derivadas em cada ponto da curva aumentam com o tempo, indicando um aumento da aceleração do foguete com o tempo. Conforme os valores de bL se tornam maiores as curvas mudam de tal forma que as derivadas em cada ponto passam a diminuir, indicando um decréscimo da aceleração do foguete com o tempo. Estes são os casos das curvas para bL = 0,01, 0,03 e 0,05 kg/s. Para que os gráficos de velocidade dos modelos III (Figura 4(b)) e IV (Figura 4(c)) estejam na faixa de 0–300 m/s (para se realizar uma comparação com os modelos I e II) é necessário que os valores de bQ sejam inferiores a uma potência de pelo menos 10−2. Por representarem modelos mais realísticos para movimento de foguetes, estes modelos puderam ser simulados com a mesma faixa de valores da constante bQ. Inicialmente, aparenta-se não haver distinção entre as formas das curvas, considerando um mesmo valor de b dentre esses modelos. Observa-se, entretanto, que os valores de velocidade nos tempos tq e ta são ligeiramente diferentes, o que acarreta em apogeus distintos para estes modelos (Figuras 4(e) e (f)). Para efeitos de comparação, são apresentadas na Figura 5(a) as curvas v(t) para os modelos III e IV com os valores de bQ = 0,00005 kg/m e 0,0002 kg/m, e na Figura 5(b), as respectivas curvas y(t). Observa-se que, apesar da semelhança marcante entre as curvas (para ambos os valores de bQ), os valores de velocidade para o modelo III são superiores aos do modelo IV em todos os instantes de tempo e, como resultado disto, as altitudes para este modelo são superiores às do modelo IV. A mesma análise é válida para as curvas com bQ = 0,0001 kg/m e 0,0005 kg/m. Figura 5 Comparação para as curvas (a) v(t) e (b) y(t) para os modelos III e IV. Foram utilizados valores de bQ = 0,00005 kg/m e 0,0002 kg/m. Os modelos II, III e IV supracitados parecem todos serem úteis para descrever o movimento vertical ascendente de foguetes. Entretanto, sabe-se que o modelo II não se adequa ao movimento de foguete, pois este modelo se aplica a objetos que se deslocam em meios viscosos e adaptável para situações nas quais as velocidades são inferiores a 24 m/s [18]. O minifoguete Epsilon-8 atinge uma velocidade ≈128,3 m/s até o instante tq e como será visto, o modelo I não é muito ideal para simular o movimento deste minifoguete. São feitas, na próxima seção, análises comparativas com um voo real (experimental) do minifoguete Epsilon-8, abordado nas seções 2.1 e 2.4. , foram plotadas até seus respectivos tempos de apogeu. As curvas contínuas nos gráficos da Figura 6(a) representam a melhor simulação para v(t) e y(t) utilizando o modelo II, enquanto as Figuras 6(b) e (c) representam as melhores simulações para os modelos III e IV respectivamente. No método de ajuste adotado, para os modelos II a IV, os parâmetros ve,R,tq e mo foram mantidos fixos e nos valores descritos na seção 2.4 2.4. Parâmetros teóricos/experimentaisrelacionados ao minifoguete Epsilon-8 Como descrito na seção 2.1 foi utilizado o minifoguete Epsilon-8 com motor a propelente sólido classe E7. Parâmetros de importância sobre as características do motor classe E7 podem ser obtidos experimentalmente através de testes estáticos [22]. Nestes testes, o motor é acoplado a um anteparo que ligado a um sistema eletrônico, fornece a força/empuxo em função do tempo [T(t)], quando o motor é ignitado. A Figura 3 representa um dos testes estáticos realizado, dentre 6 obtidos pelo GFCS da UFPR. O valor do empuxo médio, representado em vermelho, foi de 6,58 N. Figura 3 Curva experimental de empuxo em função do tempo para o motor E7 utilizado no minifoguete. No Gráfico da Figura 3, com os pontos experimentais representados por círculos, a magnitude do impulso total (Itot) é definida de acordo com a Equação (36) como a integral de [T(t)], durante o tempo de funcionamento do motor, ou seja, até o tempo final de queima (tq) do propelente: (38) I t ⁢ o ⁢ t = ∫ 0 t q T ⁢ ( t ) ⁢ d t . A magnitude do empuxo médio (Tm) que o motor exerce sobre o minifoguete é calculado pela Equação: (39) T m = I t ⁢ o ⁢ t t q . As equações de movimento descritas nas seções 2.3.1 a 2.3.4 consideram um empuxo constante, igual ao valor médio do empuxo real calculado através da Equação (39) e representado pela linha vermelha na Figura 3. Esta aproximação é razoável no chamado regime estacionário e quando os regimes transientes (início e término do funcionamento do motor) são breves [23]. Assim, a magnitude do impulso total pode ser calculada como: (40) I t ⁢ o ⁢ t = v e ⁢ | d ⁢ m d ⁢ t | ⁢ t q . Outro parâmetro relevante é o impulso específico (Isp), obtido pela divisão do impulso total pelo peso do propelente (mpg) de acordo com a equação: (41) I s ⁢ p = I t ⁢ o ⁢ t m p ⁢ g , na qual g é a magnitude da aceleração da gravidade ao nível do mar. Com esta grandeza mede-se a eficiência do propelente, já que, quanto maior for o impulso específico, maior será a variação da quantidade de movimento, por unidade de peso do propelente. Admite-se neste trabalho que o empuxo T(t) é constante, [4]. Desta forma, considerando-se que ve e dm/dt são constantes, é possível estimar a velocidade de exaustão dos gases de acordo com a equação: (42) v e = I t m p = I s ⁢ p ⁢ g . Mede-se experimentalmente tq e mp sendo possível estimar R = dm/dt, como: (43) R = m p / t q . A Tabela 1 fornece a massa do propelente utilizado no minifoguete Epsilon-8 e os parâmetros médios relevantes obtidos dos testes estáticos do motor. Nesta tabela também se encontram outros parâmetros experimentais relacionados ao minifoguete e ao voo do mesmo. Outros parâmetros como velocidade máxima e altitude máxima atingida serão discutidos na análise gráfica da seção 3. Tabela 1 Parâmetros do motor do minifoguete Epsilon-8. A grandeza mp é a massa do propelente, It é o valor médio do impulso total, obtido do teste estático de 6 motores idênticos e tq é o tempo de queima do propelente; moé a massa inicial do minifoguete (antes da queima do propelente), mfé a massa final do mesmo e taé o tempo necessário para atingir o apogeu. mp(g) It(médio)(N.s) tq (s) mo (g) mf(g) ta(s) 44,70 ± 0,01 34,90 ± 0,01 5,358 ± 0,001 131,20 ± 0,01 86,50 ± 0,01 11,7 ± 0,1 Utilizando a Equação (42) e os valores de It (médio) e mp da Tabela 1 é possível estimar ve = (780,8 ± 0,4) m/s. Adicionalmente, com os valores de mp e tq da Tabela 1 é possível estimar, utilizando a Equação (43), R = (8,343 ± 0,003) g/s. . Com isto, são testados estes modelos na condição aproximada de impulso constante e a validade das Equações (42) e (43) na descrição das variáveis ve e R respectivamente. As constantes de arrasto bL e bQ foram variadas para que as funções v(t) e y(t) teóricas se ajustassem aos pontos experimentais destas respectivas funções. Durante o método de ajuste adotado, as constantes b’s foram modificadas em cada um dos modelos para que no tempo tq = 5,358 s a velocidade de queima teórica, vqteor, fosse da ordem do valor experimental, vqexp=128,3 m/s. Feito isto, foi observado se no instante do apogeu experimental taexp=11,7 s o valor máximo da curva teórica ymáx se aproximava do valor do apogeu do minifoguete ≈723 m (veja a seção 2.1 2.1. O minifoguete Um minifoguete, assim como um foguete real, é um veículo que se desloca expelindo, em sua parte traseira, um fluxo de gás à alta velocidade devido à queima do propelente [11]. De acordo com a terceira lei de Newton, esta ação resulta em uma reação de mesma intensidade, porém sentido oposto, deslocando o foguete para cima a partir do solo [4]. Minifoguetes são utilizados no ensino e em competições. Já outros, mais sofisticados do tipo sondagem, são aplicados às pesquisas aeroespaciais cujas altitudes são restritas à troposfera [12]. Na Figura 1(a) está ilustrada a seção transversal de um minifoguete, utilizado em atividades de ensino, com as principais partes que o compõe. Na Figura 1(b) é mostrada uma fotografia do minifoguete utilizado neste trabalho. Figura 1 (a) Ilustração do corte de seção transversal de um minifoguete. (b) Minifoguete Epsilon-8. Como se observa na Figura 1(a) a estrutura do minifoguete é delimitada pela linha em preto e inclui a coifa (nariz), a parte traseira que sustenta as empenas (ou aletas) e o motor demarcado pela linha azul. O motor contém o propelente na cor amarela que está no interior de um cilindro por onde é feito um pequeno furo em sua extremidade direita [13]. A queima do propelente é feita por um dispositivo de ignição elétrica que contém um ignitor e dois condutores. Estes ao serem ligados a uma pequena bateria transmitem energia que inflama o ignitor e consequentemente o propelente. Um dispositivo de ejeção do motor (não ilustrado aqui) desacopla o conjunto constituído pela coifa e compartimento do altímetro, o que automaticamente aciona o paraquedas. Um elástico segura todo o conjunto do foguete para um pouso seguro. O minifoguete utilizado (Figura 1(b)) para a obtenção dos dados experimentais foi o Epsilon-8. A ele foi acoplado um motor comercial de classificação E7 segundo a nomenclatura da NAR (National Association of Rocketry) [14]. Nesta nomenclatura a letra especifica que o impulso total produzido pelo motor está na faixa de 20,00 a 40,00 N.s, e o número em seguida especifica o empuxo médio produzido de 7 N. De fabricação do Grupo de Foguetes Carl Sagan (GFCS) da UFPR [8], o Epsilon-8 possui diâmetro máximo de 2,03 cm, comprimento de 50,1 cm e massa de 131,2 g [15]. Ele é o atual recordista brasileiro de apogeu para a classe E, atingindo em 2015 um apogeu de 723 m. Os dados deste voo foram obtidos com um altímetro micro-peak da Altus Metrum [9] e foram utilizados, neste trabalho, para os ajustes das equações de movimento em dados experimentais reais de v(t) e y(t). Os dados do voo foram obtidos em 08/08/2015 pelo Grupo de Foguetes Carl Sagan da Universidade Federal do Paraná (UFPR), que é liderado por Carlos Henrique Marchi, um dos autores do presente trabalho. ).

Figura 6
Representação dos pontos experimentais (círculos abertos) e as curvas simuladas (linhas contínuas) para v(t) e y(t) do voo do minifoguete Epsilon-8. Em (a) simulação para o modelo II, (b) modelo III e (c) modelo IV.

Seguindo o método supracitado, ao se escolher bL = 0,0358 kg/s para o modelo II (Figura 6(a)), observa-se que a curva teórica para v(t) está, na média, abaixo dos pontos experimentais para praticamente todos os instantes no tempo. O mesmo ocorre com a curva para y(t) que prevê um apogeu de ymáx = 633,5 m em t = 9,4 s. Este apogeu representa uma diferença percentual de 12,4% relativo ao apogeu experimental. Embora este erro não seja “grande”, os ajustes indicam que a aproximação utilizando uma força de arrasto linear com a velocidade durante todo o voo, não se adequa bem aos dados experimentais do minifoguete e assim outros modelos precisam ser trabalhados.

Seguindo a mesma sistemática, foi escolhido bQ = 0,000306 kg/m para o modelo III. Observa-se uma maior aproximação da curva teórica para v(t) aos dados experimentais para o intervalo de tempo 0≤t≤1,5 s. Entretanto, para o intervalo de 1,5≤ttq, a curva teórica está acima dos dados experimentais. O mesmo ocorre para ttq. Como resultado, observa-se um bom ajuste da curva teórica para y(t) aos dados experimentais para o intervalo 0≤t≤5,0 s. Acima deste intervalo, inicia-se uma pequena divergência dos valores teóricos aos dados experimentais, de tal forma que o apogeu teórico previsto é ymáx = 785,3 m, que coincide como instante experimental ta = 11,7 s. Este apogeu representa uma diferença percentual de 8,6% relativo ao apogeu experimental do minifoguete. Assim, o modelo III, que considera a massa constante durante todo o voo de subida, é razoável para simular o movimento do Epsilon-8.

Para efeitos de comparação, o ajuste do modelo IV foi realizado com o mesmo valor de bQ utilizado no modelo III. Com exceção do intervalo 2,6 s ≤ttq a curva teórica para v(t) é dentre os modelos estudados, a que melhor se ajusta aos dados experimentais. Observa-se também na curva y(t), melhor ajuste da curva aos pontos experimentais. O apogeu teórico, previsto para o instante ta = 11,7 s foi de 731,3 m, o que representa uma diferença percentual de apenas 1,5% relativo ao apogeu experimental. O modelo IV, como relatado na seção 2.3.4 2.3.4. Modelo IV: Solução completa com forçade arrasto quadrática Uma solução analítica completa da Equação (24) para v(t), durante a queima do propelente, foi obtida por Rodrigues et al. (2014) [5]. O autor demonstrou, utilizando técnicas de cálculo diferencial, que a Equação (3) pode ser reescrita em termos de uma equação diferencial cujas soluções são funções de Bessel. Veja inicialmente que a Equação (24), agora com m = mo−Rt, pode ser escrita como: (31) d ⁢ v d ⁢ t + b m ⁢ v 2 - R ⁢ v e m + g = 0 , Durante as etapas para a solução da Equação (31), é possível realizar uma mudança da variável v, para uma nova variável χ [5]. Neste caso, esta equação se transforma em: (32) ξ 2 ⁢ d 2 ⁢ χ d ⁢ t 2 + ξ ⁢ d ⁢ χ d ⁢ t + ( ξ 2 - ν 2 ) ⁢ χ = 0 , que é uma Equação de Bessel de ordem ν [5, 19]. As constantes ξ e ν são definidas como: (33) ξ = ξ ⁢ ( m ) = 2 R ⁢ b ⁢ m ⁢ g , (34) ν = 2 ⁢ b ⁢ v e R . A solução geral da Equação (32) é uma combinação linear do tipo: (35) χ ⁢ ( ξ ) = C 1 ⁢ J ν ⁢ ( ξ ) + C 2 ⁢ Y ν ⁢ ( ξ ) , no qual Jν(ξ) e Yν(ξ) são funções de Bessel de primeira e segunda ordens respectivamente [5, 19, 20], C1 e C2 são constantes reais. Após algumas álgebras utilizando propriedades convenientes das funções de Bessel [5], é possível escrever a solução geral para a velocidade como: (36) v I ⁢ V ⁢ A ⁢ ( ξ ) = R 2 ⁢ b ⁢ [ ξ ⁢ C ⁢ Y ν + 1 ⁢ ( ξ ) + J ν + 1 ⁢ ( ξ ) C ⁢ Y ν ⁢ ( ξ ) + J ν ⁢ ( ξ ) - v o ] , (37) c ⁢ o ⁢ m ⁢ C = ( 2 ⁢ b ⁢ v o + R ⁢ ν ) ⁢ J ν ⁢ ( ξ o ) - R ⁢ ξ o ⁢ J ν + 1 ⁢ ( ξ o ) R ⁢ ξ o ⁢ Y ν + 1 ⁢ ( ξ o ) - ( 2 ⁢ b ⁢ v o + R ⁢ ν ) ⁢ Y ν ⁢ ( ξ ) na qual vo é a velocidade inicial do foguete, que pode ser definida como zero e ξo=2R⁢b⁢mo⁢g. Tendo em vista que ξ = ξ(m), Equação (36), e que m = m(t), Equação (5), é possível por meio da Equação (36) construir um gráfico da velocidade em função do tempo para (0≤t≤tq), etapa A, conhecendo-se os parâmetros mo,g,R,b e ve, previamente obtidos de cálculos teóricos/experimentais. A construção do gráfico y(t) pode ser obtida a partir de v(t) utilizando programas de computador adequados. Para este modelo, foi utilizado o programa GNU Octave [21] para simular as funções v(t) e y(t); a partir da velocidade foi obtida a altura através de um processo de integração numérica denominado quadratura de Gauss. No Apêndice I estão descritas as funções utilizadas no GNU Octave, bem como os scriptsimplementados para simulação do modelo apresentado nesta seção. No Apêndice II foi deixada uma explicação sucinta sobre a quadratura de Gauss. Para a etapa B neste modelo as Equações (28) e (29) podem ser aplicadas. Nesta condição deve-se considerar que m = m(tq). , é o modelo mais completo, dentre as limitações dos modelos aqui estudados, e assim, para a simulação do voo do Epsilon-8, é o mais adequado dentre os modelos.

Uma grandeza física de grande interesse no movimento de foguetes é o coeficiente de arrasto, CD. Ele é importante na determinação da constante bQ em movimentos cuja força de arrasto é proporcional ao quadrado da velocidade. As grandezas b e CD estão relacionadas de acordo com a Equação [1717. J. Waters, Undergraduate Journal of Mathematical: Modeling onde + Two 6, 1 (2014).]:

(44) b Q = ρ C D A 2 ,

na qual ρ é a massa específica e A é a área de seção transversal do foguete, normal à direção do movimento. Para um modelo aproximado, no qual o ar é considerado um gás ideal, sua densidade pode ser calculada como:

(45) ρ = P R e s p T ,

em que P e T representam os valores absolutos de pressão e temperatura local e Resp o valor específico da constante dos gases ideais. No momento do lançamento do Epsilon-8, P = 90.760 Pa, T = 296,05 K. Com Resp = 286,9 J/kg.K obtém-se ρ≈1,069 kg/m3. Com o valor de bQ = 0,000306 kg/m calculado do modelo IV e considerando que o Epsilon-8 possui A = 0,000324 m2, o valor calculado para o coeficiente de arrasto, utilizando a Equação (44) é CD = 1,77. Este valor é mais realista que o valor genérico comumente utilizado ≈0,75 [1717. J. Waters, Undergraduate Journal of Mathematical: Modeling onde + Two 6, 1 (2014).] pois o Epsilon-8 apresentava três empenas grandes que influenciam no valor deste coeficiente.

5. Considerações Finais

Neste trabalho foram apresentados 4 modelos para o voo de foguetes na atmosfera. Nas simulações de voo do minifoguete Epsilon-8 foram mantidas todas as características experimentais deste minifoguete, possibilitando descrever funções teóricas para v(t) e y(t) para os modelos de II a IV. O modelo I, sem qualquer força de arrasto, foi útil para descrever um movimento onde não existe atmosfera. Ao ser comparado com os demais modelos, que consideram uma força de arrasto, foi observado que os valores de velocidade e altitudes estão bem acima dos modelos que consideram o arrasto, mesmo para baixos valores da constante b. No modelo II, a função v(t) teórica para a constante bL = 0,0358 kg/s escolhida não se ajustou bem aos valores experimentais. Isto prevê altitudes menores quando comparado à curva experimental para o minifoguete Epsilon-8. A diferença entre a curva teórica e a experimental se agrava ainda mais para instantes no tempo próximo ao apogeu, indicando que este modelo não é adequado para descrever o movimento vertical ascendente de um minifoguete.

Os modelos III e IV, mais realísticos por representarem um arrasto quadrático com a velocidade, resultam em erros menores no apogeu do minifoguete. A curva v(t) para o modelo III passa ligeiramente acima dos pontos experimentais, prevendo apogeu com erro de 8,6% para a escolha da constante bQ em 0,000306 kg/m. Já a curva de v(t) para o modelo IV, mais completo para as equações de movimento, passa mais próxima aos pontos experimentais quando comparada ao modelo III, prevendo um erro de apenas 1,5% no apogeu tornando este modelo mais adequado para descrever o movimento ascendente do minifoguete na atmosfera. Mesmo assim, o modelo III por considerar cálculos menos complexos, torna-se útil no ensino das equações teóricas do movimento de um minifoguete.

6. Trabalhos futuros

Pretende-se incluir nas simulações das equações de v(t) e y(t) o método numérico que utiliza a curva de empuxo em função do tempo (T(t)). Neste método estima-se o empuxo do minifoguete em cada ponto da curva T(t), e a partir da equação diferencial de movimento do foguete, é possível estimar a velocidade e a altura em cada instante de voo. Neste novo trabalho serão considerados os modelos nos quais a massa do foguete é constante (como o modelo III aqui descrito), outro no qual a massa varia linearmente com o tempo em todos os instantes, e um modelo mais realista no qual a massa é variável e extraída da curva de empuxo do motor.

Agradecimentos

A equipe de trabalho agradece à PROGRAD/UFES e a FAPES pelo financiamento desta pesquisa. O professor Carlos Henrique Marchi agradece em especial ao CNPq pela concessão de sua bolsa de pesquisador.

References

Datas de Publicação

  • Publicação nesta coleção
    15 Fev 2021
  • Data do Fascículo
    2021

Histórico

  • Recebido
    20 Nov 2020
  • Revisado
    06 Jan 2021
  • Aceito
    09 Jan 2021
Sociedade Brasileira de Física Caixa Postal 66328, 05389-970 São Paulo SP - Brazil - São Paulo - SP - Brazil
E-mail: marcio@sbfisica.org.br