Modelado y simulación simplificada de reentradas atmosféricas

UNIVERSIDAD DE VALLADOLID ESCUELA DE INGENIERIAS INDUSTRIALES Grado en Ingeniería en Tecnologías Industriales Modelado y Simulación Simplificada de Reentradas Atmosféricas Autor: Galván Lobo, Jorge Tutor: Portillo de la Fuente, Ana Valladolid, Julio 2021. AGRADECIMIENTOS Me gustaría comenzar agradeciendo a mi tutora, la Profª. Dª. Ana Portillo de la Fuente, por la confianza depositada para realizar este trabajo. Quiero agradecer que me permitiera investigar y trabajar sobre una materia de un gran interés para mí a nivel personal. Quiero agradecer también su disposición completa para ayudarme a progresar, mejorar y aprender a lo largo del desarrollo de este trabajo. Si este proyecto ha llegado a buen puerto es, en gran parte, por las sugerencias, correcciones e ideas que a lo largo del desarrollo de este trabajo ha realizado. Quiero agradecer también la ayuda prestada por el Prof. D. Francisco Castro Ruiz en materia de aerodinámica y de la Profª Dª. Smaranda Podar Cristea por sus sugerencias para enfocar el diseño del controlador. Finalmente, quisiera agradecer el apoyo constante, no solo durante el desarrollo de este trabajo, sino durante todas las etapas que han venido antes, a mis padres y hermanos, así como mi amigo y compañero de mil batallas durante estos cuatro años, Álvaro Alonso. No podría haber completado este trabajo y esta etapa sin vosotros. Resumen Desarrollar las herramientas de diseño y simulación de entradas atmosféricas es vital para el planteamiento de misiones tripuladas a Marte y otros planetas. Aunque ya existen controladores complejos y robustos para estas situaciones, este trabajo trata de abordar el problema buscando una solución simplificada que permita dar buenas aproximaciones que puedan servir de primeras estimaciones de métodos y simuladores más potentes, reduciendo el coste operacional de estos. Con este objetivo, se analizará la teoría de control en entradas atmosféricas y datos empíricos de misiones anteriores para desarrollar un programa capaz de generar un controlador sencillo pero eficaz. También se simularán y validarán este programa y sus productos. Palabras claves Entrada atmosférica, Sistema de control, Simulación, Método de Dormand-Prince, Coordenadas preliminares de reentrada Abstract Developing the design and simulation tools for atmospheric entry is vital for the approach of manned missions to Mars and other planets. Even though more complex and robust controllers already exist, this project attempts to approach the problem searching for a simplified solution that allows to obtain good approximations that could serve as initial estimations for the more powerful methods and simulators, reducing the operational costs of these. To this avail, the theory for atmospheric entry and empiric data from past missions will be analyzed to develop a program capable of generating a simple but efficient controller. This program and its products will be simulated and validated as well. Keywords Atmospheric entry, Control system, Simulation, Dormand-Prince Method, Preliminary reentry coordinates INDICE Contenido 1. Introducción.................................................................................................................................... 1 1.1. Prefacio ................................................................................................................................... 1 1.2. Objetivos ................................................................................................................................. 2 1.3. Planteamiento del problema y Metodología ........................................................................ 2 2. Modelado del Sistema ................................................................................................................... 5 2.1. Modelado atmosférico ........................................................................................................... 5 2.2. Modelo gravitatorio .............................................................................................................. 12 2.3. Modelo Termodinámico ....................................................................................................... 13 2.3.1. Modelo Dinámico de 4 variables ................................................................................ 13 2.3.2. Modelo Dinámico de 6 variables ................................................................................ 15 2.3.3. Modelo Térmico ............................................................................................................ 21 2.4. Modelado de Vehículo ......................................................................................................... 22 2.5. Implementación del modelo ................................................................................................ 23 3. Control .......................................................................................................................................... 25 3.1. Método de entrada .............................................................................................................. 25 3.2. Etapas de entrada ............................................................................................................... 25 3.3. Restricciones y Limitaciones ............................................................................................... 31 3.4. Variables de Control ............................................................................................................. 32 3.5. Arquitectura de programa ................................................................................................... 33 3.6. Validación de control ........................................................................................................... 37 4. Configuración y ajuste del Simulador ......................................................................................... 40 4.1. Simulador ............................................................................................................................. 40 4.2. Modulación. Evaluación de parámetros ............................................................................. 42 5. Resultados de Simulación ........................................................................................................... 58 5.1. Caso 1 ................................................................................................................................... 58 5.2. Caso 2 ................................................................................................................................... 68 5.3. Caso 3 ................................................................................................................................... 79 6. Conclusiones, Repercusiones y líneas futuras .......................................................................... 92 6.1. Conclusiones ........................................................................................................................ 92 6.2. Repercusiones ..................................................................................................................... 92 6.3. Líneas Futuras ...................................................................................................................... 93 7. Bibliografía ................................................................................................................................... 94 7.1. Lista de referencias ............................................................................................................. 94 7.2. Comentario Bibliográfico ..................................................................................................... 95 ANEXO A: Código de implementación de modelo .............................................................................. 98 ANEXO B: Código de algoritmo de diseño de control ....................................................................... 101 ANEXO C: Código de simulación ........................................................................................................ 106 ANEXO D: Código de implementación de métodos numéricos ....................................................... 108 Indice de Ilustraciones Ilustración 2. 1 Comparativa de Densidad Mars-GRAM NASA ......................................................... 7 Ilustración 2. 2 Comparativa de Temperatura Mars-GRAM NASA ................................................... 7 Ilustración 2. 3 Comparativa de Densidad Mars-GRAM MCD ........................................................ 11 Ilustración 2. 4 Comparativa de Temperatura Mars-GRAM MCD .................................................. 11 Ilustración 2. 5 Sistema de Coordenadas del modelo de 4 variables ........................................... 14 Ilustración 2. 6 Sistema de coordenadas fijo al centro planetario ................................................ 16 Ilustración 2. 7 Sistema de coordenadas solidario a la rotación planetaria ................................. 17 Ilustración 2. 8 Sistema de coordenadas solidario al vehículo ..................................................... 18 Ilustración 2. 9 Sistema de coordenadas referenciado a la velocidad del vehículo ..................... 19 Ilustración 2. 10 Declaración de la función System ...................................................................... 23 Ilustración 2. 11 Codigo para la entrada en la tabla de referencia de variables de estado atmosférica .................................................................................................................................... 23 Ilustración 2. 12 Codigo de implementación gravedad, arrastre y sustentación. ......................... 24 Ilustración 2. 13 Codigo de implementación del sistema de ecuaciones diferenciales ................ 24 Ilustración 3. 1 Etapas presentes en las trayectorias "Skip-Out" ................................................... 26 Ilustración 3. 2 Historico de variación del ángulo de inclinación lateral en la reentrada de la Apollo10 ......................................................................................................................................... 27 Ilustración 3. 3 Histórico de variación del ángulo de inclinación lateral en la reentrada de la Apollo10 ......................................................................................................................................... 28 Ilustración 3. 4 Perfil de velocidad de la cápsula Apollo10 durante la reentrada ......................... 30 Ilustración 3. 5 Declaración de la función Control ......................................................................... 33 Ilustración 3. 6 Código de generación del vector de valores iniciales ........................................... 34 Ilustración 3. 7 Código de declaración de las variables de control ............................................... 34 Ilustración 3. 8 Código de implementación del controlador .......................................................... 37 Ilustración 3. 9 Comparativa de perfiles de altitud Simulador vs Apollo10 ................................... 38 Ilustración 3. 10 Comparativa de perfiles de velocidad Simulador vs Apollo10 ........................... 39 Ilustración 4. 1 Tablero de Butcher Método Dormand-Prince ....................................................... 40 Ilustración 4. 2 Comparativa de perfiles de altitud creados por diferentes solvers numéricos .... 42 Ilustración 4. 3 Influencia de la primera actuación de control en la trayectoria ........................... 43 Ilustración 4. 4 Influencia de la segunda actuación de control en la trayectoria .......................... 44 Ilustración 4. 5 Influencia de la tercera actuación de control en la trayectoria ............................ 45 Ilustración 4. 6 Influencia de la cuarta actuación de control y porsteriores en la trayectoria ....... 46 Ilustración 4. 7 Influencia de la segunda actuación de control en la longitud .............................. 47 Ilustración 4. 8 Influencia de la segunda actuación de control en la latitud ................................. 47 Ilustración 4. 9 Influencia de la tercera actuación de control en la longitud................................. 47 Ilustración 4. 10 Influencia de la tercera actuación de control en la latitud ................................. 48 Ilustración 4. 11 Influencia de la cuarta actuación de control y posteriores en la longitud .......... 48 Ilustración 4. 12 Influencia de la cuarta actuación de control y posteriores en la latitud ............ 48 Ilustración 4. 13 Influencia en la trayectoria del parametro de modulación Tp ............................ 49 Ilustración 4. 14 Influencia en la longitud del parametro de modulación Tp ................................ 50 Ilustración 4. 15 Influencia en la latitud del parametro de modulación Tp ................................... 51 Ilustración 4. 16 Influencia en la velocidad del parametro de modulación Tp .............................. 52 Ilustración 4. 17 Influencia de la altura de interfase en la trayectoria .......................................... 53 Ilustración 4. 18 Influencia de la velocidad inicial en la trayectoria .............................................. 54 Ilustración 5. 1 Influencia del ángulo de rumbo inicial en el perfil de longitud en el Caso 1 ........ 59 Ilustración 5. 2 Influencia del ángulo de rumbo inicial en el perfil de latitud en el Caso 1 ........... 59 Ilustración 5. 3 Perfil de variación del ángulo de inclinación lateral en el Caso 1 ........................ 60 Ilustración 5. 4 Trayectoria de entrada 3D en el Caso 1 ............................................................... 61 Ilustración 5. 5 Perfil de altura en el Caso 1 .................................................................................. 62 Ilustración 5. 6 Perfil de velocidad en el Caso 1 ............................................................................ 63 Ilustración 5. 7 Perfil de longitud en el Caso 1 .............................................................................. 64 Ilustración 5. 8 Perfil de latitud en el Caso 1 ................................................................................. 65 Ilustración 5. 9 Perfil del ángulo de rumbo en el Caso 1 ............................................................... 66 Ilustración 5. 10 Perfil del ángulo de trayectoria en el Caso 1 ...................................................... 67 Ilustración 5. 11 Flujo térmico al vehículo en el Caso 1 ................................................................ 68 Ilustración 5. 12 Influencia del ángulo de rumbo inicial en el perfil de longitud en el Caso 2 ...... 69 Ilustración 5. 13 Influencia del ángulo de rumbo inicial en el perfil de latitud en el Caso 2 ........ 70 Ilustración 5. 14 Perfil de variación del ángulo de inclinación lateral en el Caso 2 ...................... 71 Ilustración 5. 15 Trayectoria de entrada 3D en el Caso 2 ............................................................. 72 Ilustración 5. 16 Perfil de altura en el Caso 2 ............................................................................... 73 Ilustración 5. 17 Perfil de velocidad en el Caso 2 ......................................................................... 74 Ilustración 5. 18 Perfil de longitud en el Caso 2 ............................................................................ 75 Ilustración 5. 19 Perfil de latitud en el Caso 2 ............................................................................... 76 Ilustración 5. 20 Perfil del ángulo de rumbo en el Caso 2 ............................................................. 77 Ilustración 5. 21 Perfil del ángulo de trayectoria en el Caso 2 ...................................................... 78 Ilustración 5. 22 Flujo térmico al vehículo en el Caso 2 ................................................................ 79 Ilustración 5. 23 Influencia del ángulo de rumbo inicial en el perfil de longitud en el Caso 3 ...... 80 Ilustración 5. 24 Influencia del ángulo de rumbo inicial en el perfil de latitud en el Caso 3 ........ 81 Ilustración 5. 25 Perfil de variación del ángulo de inclinación lateral en el Caso 3 ...................... 82 Ilustración 5. 26 Trayectoria de entrada 3D en el Caso 3 ............................................................. 83 Ilustración 5. 27 Perfil de altura en el Caso 3 ............................................................................... 84 Ilustración 5. 28 Perfil de velocidad en el Caso 3 ......................................................................... 85 Ilustración 5. 29 Perfil de longitud en el Caso 3 ............................................................................ 86 Ilustración 5. 30 Perfil de latitud en el Caso 3 ............................................................................... 87 Ilustración 5. 31 Perfil del ángulo de rumbo en el Caso 3 ............................................................. 88 Ilustración 5. 32 Perfil del ángulo de trayectoria en el Caso 3 ...................................................... 89 Ilustración 5. 33 Flujo térmico al vehículo en el Caso 3 ................................................................ 90 Indice de Tablas Tabla 1. 1 Datos Atmosféricos MARS-GRAM.................................................................................... 6 Tabla 1. 2 Densidad-Altura MCD ...................................................................................................... 9 Tabla 1. 3 Temperatura-Altura MCD .............................................................................................. 10 Tabla 2. 1 Datos técnicos Capsula ORION ..................................................................................... 22 Tabla 4. 1 Comparativa de tiempos de ejecucción de solver ........................................................ 41 1 1. Introducción 1.1. Prefacio Desde la antigüedad, la humanidad ha mirado al cielo y ha tratado de explicar lo que veía. En las constelaciones se leían historias, en las estrellas se veía la inmutabilidad de una bóveda fija y en los cometas presagios del futuro. Con el tiempo, hemos aprendido que las estrellas son grandes esferas gaseosas encendidas por la fusión de su núcleo; que el espacio es un lugar de enormes dimensiones, en constante cambio y expansión del que conformamos la infinitésima parte. Nuestro concepto del espacio que rodea la esfera rocosa que llamamos hogar ha cambiado con los siglos, pero nuestra fascinación con él no ha hecho sino aumentar. Afortunadamente, desde 1957, cuando el satélite Sputnik 1 fue puesto en órbita, la humanidad ha comenzado a explorar el espacio además de observarlo. Este afán exploratorio ha acelerado nuestros descubrimientos y alimentado el desarrollo de una de las industrias que actualmente ostenta proyecciones de demanda y crecimiento al alcance de muy pocos sectores [1]. Desde mediados del siglo XX, el sector espacial ha ido progresando, abriéndose hueco a nivel industrial y a nivel social. Hoy en día, la sociedad ve la exploración espacial como un tema de interés y no una curiosidad. Y las empresas y organizaciones del sector como la ESA o la NASA, antes controladas casi por completo por gobiernos, compiten ahora contra y con la iniciativa privada. Empresas como Blue Horizon, SpaceX o Thales Alenia Space, compiten ahora por ser las siguientes en hacer avanzar a la humanidad en su odisea por el espacio. Sin embargo, la incorporación de nuevos agentes a la actualidad del sector espacial no ha cambiado el planteamiento fundamental de sus acciones. La industria espacial tiene por objetivo idear, perfeccionar y fabricar tecnología, con la intención de poner parte de ésta en órbita o, en su defecto, en algún otro lugar del espacio; y en ciertos casos, hacer que esta tecnología retorne a la Tierra [2]. En particular, esta tecnología se compone mayoritariamente de satélites, telescopios y otros instrumentales diseñados para funcionar en órbita. Pero, en los últimos años, un tipo de misión, que otrora fue protagonista con misiones como Vostok 1 (1961) o Apollo 11 (1969), ha vuelto a tomar relevancia. Las misiones tripuladas, que durante años han consistido en viajes y misiones a la Estación Espacial Internacional (ISS), han tomado otra perspectiva con misiones como ExoMars RSP (2022) [3] o Artemis (2024) [4]. La misión ExoMars es una misión desarrollada por la ESA que pretende desplegar en Marte un Rover y una plataforma de superficie para explorar la superficie marciana y buscar señales de vida. La misión Artemis es un conjunto de misiones y operaciones desarrolladas por NASA que tienen por objetivo llevar astronautas de vuelta a la Luna y preparar la infraestructura necesaria para llevar humanos a Marte en el 2030. Es observando estas misiones, así como los viajes a la ISS y el despliegue de Rovers como Soujourner, Pathfinder (1997), Curiosity (2011) o Perseverance (2021), que se puede encontrar un punto en común de todas ellas. Una conexión que no comparten las misiones basadas en satélites. Todas estas misiones requieren en uno u otro momento de una entrada atmosférica para tener éxito. Las entradas atmosféricas son operaciones de una elevada complejidad y crucial importancia. El correcto desarrollo de una entrada atmosférica es crítico para el éxito de una misión. Es por 8 mayores, no modela correctamente la temperatura y tampoco tiene en cuenta las variaciones debido al día y la noche. Es, sin embargo, a partir de esta comprobación que se procede a buscar el siguiente modelo. Nos fijamos ahora en el propio modelo de comprobación. El MARS-GRAM. Tras una búsqueda exhaustiva, usando como palabras clave: MARS-GRAM, Modelo atmósfera, Marte. Los resultados resultan bastante decepcionantes. Tras una revisión de más de una decena de artículos encontrando poco más que las tablas para condiciones estáticas y particulares de las que ya se disponía, se accede a una página web perteneciente a la NASA. En ella se explica la imposibilidad de acceder al modelo sin el permiso de esta organización o de ciertos departamentos del gobierno estadounidense. En esta página se pone de manifiesto la poca probabilidad de obtener acceso al modelo sin pertenecer a la organización. Es por ello por lo que se desiste de obtener acceso al modelo. Sin embargo, tras descartar el MARS-GRAM, se pasa a estudiar la posibilidad de emplear el MCD (Mars Climate Database) como el modelo a emplear. El MCD es un modelo atmosférico de elevada complejidad a la altura en rigor y completitud del MARS-GRAM [8, 9]. El MCD se encuentra disponible gratuitamente y con libre acceso en [10]. En esta página web se encuentran tanto los diferentes artículos científicos desarrollados paralelamente al modelo como el modelo en sí. Este modelo permite obtener datos teniendo en cuenta Latitud, Longitud, Día, Hora, Altura, … Se procede a estudiar tanto los datos disponibles directamente como la posibilidad de obtener acceso al modelo en sí. Las tabulaciones de datos entregadas por la interfaz web cuentan con el siguiente formato: 9 Tabla 1. 2 Densidad-Altura MCD 10 Tabla 1. 3 Temperatura-Altura MCD En este caso se ha empleado la secuencia de datos predichos por el MCD para la misión Pathfinder. Haciendo la misma comprobación empleada para el modelo propuesto por NASA obtenemos las siguientes relaciones 11 Ilustración 2. 3 Comparativa de Densidad Mars-GRAM MCD Ilustración 2. 4 Comparativa de Temperatura Mars-GRAM MCD Como se puede observar claramente, la correlación de densidades es algo mejor en este caso y la de temperaturas es mucho mejor. A mayores, este modelo cumple la mayoría de los objetivos secundarios planteados, dando flexibilidad al simulador sin perder precisión. 12 Solo estos resultados resultarían justificación suficiente para emplear el MCD sobre el resto de las opciones. Sin embargo, hay una razón adicional que insta a elegir este modelo, se obtuvo permiso para trabajar con el modelo directamente, sin requerir la interfaz web. Con la intención de obtener acceso se contactó la organización al cargo del modelo. Se expuso la relevancia científica del proyecto, se propuso el empleo de este modelo dada la intención de simular las entradas en Marte. Afortunadamente, esta organización encabezada por el Laboratorio de Meteorología Dinámica de París accedió a compartir el código del modelo. Gracias a ello, se pudo mejorar el modelo. Con el código, se obtuvo independencia de la versión web y se pudo desarrollar un controlador autocontenido. Elegido el modelo, analicemos la extracción de información y su tratamiento para su empleo en el código de simulación. El procedimiento a seguir será, primeramente, el traspaso de los datos tabulados por el MCD a sendas matrices conteniendo cada una las propiedades de densidad, presión, conductividad térmica y temperatura indexadas con sus respectivas alturas. A continuación, se generarán cuatro matrices, una por propiedad, de dimensiones (5n) donde n es el número de tramos de altura otorgados por el MCD. La estructura y objeto de cada columna será la siguiente: • Columna 1: Valor de la propiedad tabulada al principio del tramo n. • Columna 2: Valor de la propiedad tabulada al principio del tramo n+1. • Columna 3: Pendiente de la recta generada por los puntos determinados por los pares de propiedad y altura al principio y final del tramo. • Columna 4: Valor de la altura al principio del tramo n. • Columna 5: Valor de la altura al principio del tramo n+1. Con esto proceso obtenemos 4 matrices que contienen toda la información necesaria para modelar las propiedades de estado de la atmosfera para las coordenadas, fecha y hora que se deseen simular. (Nota: El código de implementación se encuentra disponible en los anexos). 2.2. Modelo gravitatorio El modelado de la gravedad en el planeta se realiza aplicando directamente la Ley de Gravitación Universal. Para ello, se hace una serie de suposiciones importantes: • El planeta en cuestión, Marte, es una esfera perfecta. • El planeta es una masa de densidad constante y homogénea. Tenido esto en cuenta, se calcula la gravedad superficial del planeta, lo que llamaremos g0. 𝑔0=𝐺 𝑀 𝑅02, (2.6) donde G=6.67·10-11 N2m2/kg, R0=3389 km, M=6.39·1023 kg [11, 12] 𝑔0=3.71 (𝑚/𝑠2). (2.7) 13 Hecho el cálculo de la gravedad en la superficie solo resta realizar la relación para cualquier altura con respecto a la superficie, teniendo en cuenta que la distancia influye cuadráticamente en la gravedad 𝑔(ℎ)=𝑔0(𝑅0 𝑅0+ℎ)2. (2.8) 2.3. Modelo Termodinámico 2.3.1. Modelo Dinámico de 4 variables Realizado el modelo del entorno, es preciso estudiar el modelado de la termodinámica de la primera etapa de reentrada. Dado el problema a estudiar, se opta por un modelo que desacople e independice el cálculo del flujo térmico del cálculo de las variables dinámicas. Para ello se valora la implantación de dos modelos diferentes. Un modelo de 4 variables dinámicas y un modelo de 6 variables dinámicas El modelo de 4 variables dinámicas se analiza primero con el objeto de estudiar la viabilidad de ampliar la complejidad del simulador con un sistema con más variables. Dada la menor complejidad, es interesante valorar y estudiar el desempeño del modelo más sencillo para poder tomar una decisión rigurosa sobre la adecuidad de emplear un modelo más complejo en el simulador. El modelo de 4 variables se obtiene directamente de la literatura, en particular de [12]. En esta publicación se discute y justifica en detalle el modelo, las suposiciones que emplea y el sistema de coordenadas empleado. Dado que este modelo no será el modelo finalmente empleado en el simulador sino la justificación para la implementación del modelo de 6 variables, no se procederá a explicar en detalle el modelo de 4 variables, reduciéndose la explicación a los aspectos más importantes. El modelo de 4 variables emplea el siguiente sistema de coordenadas: 14 Ilustración 2. 5 Sistema de Coordenadas del modelo de 4 variables De esta manera, podemos descomponer las fuerzas, velocidades, aceleraciones y posiciones en los ejes solidarios al cuerpo, lo que va a simplificar enormemente las ecuaciones. Se emplean las siguientes suposiciones: • Planeta esférico • Planeta con rotación nula • Aceleración inicial de vehículo = 0. Finalmente, se plantea el siguiente sistema de ecuaciones diferenciales 𝑓(1)=𝑥󰇗=𝑉𝑐𝑜𝑠𝛾, (2.9) siendo x la distancia horizontal considerando la superficie del planeta como plana en la superficie de reentrada. Esta ecuación representa la variación temporal de esta variable, por tanto, la velocidad horizontal. 15 𝑓(2)=ℎ󰇗=𝑉𝑠𝑒𝑛𝛾, (2.10) siendo h la altura con respecto a la superficie. Esta ecuación representa la variación temporal de esta variable, por tanto, la velocidad vertical perpendicular a la superficie. 𝑓(3)=𝑉󰇗=−𝜌(ℎ)𝑉2 2( 𝑚 𝐶𝐷𝑆)+𝑔(ℎ)𝑠𝑖𝑛𝛾, (2.11) siendo V la velocidad como variable vectorial cuyas componentes son 𝑥󰇗 y ℎ󰇗. Esta ecuación representa la variación temporal de V, por tanto, la aceleración de la cápsula. 𝑓(4)=𝛾󰇗=−𝜌(ℎ)𝑉(𝐿 𝐷) 2( 𝑚 𝐶𝐷𝑆)+𝑔(ℎ)𝑐𝑜𝑠𝛾 𝑉−𝑉𝑐𝑜𝑠𝛾 𝑅0, (2.12) siendo 𝛾 el ángulo de incidencia con la atmósfera, representativo de cómo de agresiva es la reentrada. Esta ecuación representa la variación temporal de esta variable. Este modelo emplea a mayores los coeficientes de arrastre (𝐶𝐷), la proporción de empuje y arrastre (L/D), la masa del vehículo (m), el radio volumétrico medio (𝑅0) y área de referencia para el cálculo aerodinámico (S). Este modelado dinámico emplea la relación L/D como parámetro de control. Manipulando esta relación, se altera el valor de 𝛾 , lo que a su vez altera la aceleración y finalmente altera las velocidades y posiciones. Este sistema es relativamente fácil de implantar en nuestro entorno de simulación y resolver empleando el modelo atmosférico ya citado. Sin embargo, como ya hemos comentado anteriormente, el objeto de este modelo es justificar el empleo del otro modelo. El sistema de 6 variables nos aporta una dimensión adicional y permite trabajar el problema desde un enfoque tridimensional. A nivel de diseño, esto dota al simulador de una gran potencia. De esta manera, mientras la implementación no sea mucho más compleja que la del modelo de 4 variables, el modelo de 6 variables será una mejor opción. Comprobemos la complejidad del modelo de 6 variables. 2.3.2. Modelo Dinámico de 6 variables El modelo de 6 variables se desarrolla siguiendo el desarrollo teórico expuesto en [13]. En particular se emplea el modelo basado en las siguientes suposiciones: • Planeta esférico • Planeta no rotatorio • Masa puntual • Masa constante • Arrastre siempre actúa en sentido contrario a velocidad • Empuje actúa siempre perpendicular a velocidad • La gravedad actúa siempre como un vector entre la masa puntual y el centro del planeta • Sistema de coordenadas solidario a planeta y solidario a vehículo comienzan alineados a t=0(s) • Los coeficientes de arrastre y empuje son constantes 16 • El vehículo se mantiene en Hipersonicidad (M>3). Vistas las suposiciones, se plantea el sistema de coordenadas. Nótese que se empleará una notación similar al modelo en el que se basa nuestro propio modelo. El primer sistema de coordenadas que se empleará será un sistema centrado en el planeta y no giratorio. Este sistema se conocerá como OXYZ. Ilustración 2. 6 Sistema de coordenadas fijo al centro planetario El segundo sistema a considerar es un sistema solidario al planeta. Esto implica que el sistema rota con el planeta y, por tanto, existe una velocidad angular relativa entre OXYZ y este nuevo sistema que se llamará O1X1Y1Z1. Será importante por tanto plantear el cambio de base para pasar de un sistema a otro. [𝑒1]=(cos (𝜔∆𝑡) sen (𝜔∆𝑡) 0 −sen (𝜔∆𝑡) cos (𝜔∆𝑡) 0 0 0 1)[𝑒]. (2.13) En particular, dada nuestra suposición de rotación nula, ω=0: 17 [𝑒1]=(100 010 001)[𝑒]. (2.14) Ilustración 2. 7 Sistema de coordenadas solidario a la rotación planetaria A continuación, se procede a plantear el sistema solidario al vehículo. Para ello, debemos ahora incorporar los ángulos de latitud y longitud que posicionan el vehículo en el sistema de coordenadas planetarias. De manera similar a como se hizo anteriormente, se plantea la longitud y la latitud como un desplazamiento de los ejes en torno a uno fijo. Esto queda patente en el diagrama siguiente: 24 Ilustración 2. 12 Código de implementación gravedad, arrastre y sustentación. Se plantean los bloques de control, en los que se contienen las actuaciones del controlador. La arquitectura de estos bloques será comentada más adelante Finalmente, se plantea el sistema de ecuaciones, empleando el vector “f” para representar la derivada del vector de estados e “y” para representar el vector de estados. El resultado es el siguiente Ilustración 2. 13 Código de implementación del sistema de ecuaciones diferenciales 25 3. Control 3.1. Método de entrada Existen diversos métodos de (re)entrada atmosférica. Cada uno cuenta con extensos estudios y literatura dedicados a ellos. De entre las múltiples opciones, tres conjuntos destacan particularmente: El planeo, el “Skip-Out” y el balístico. Es importante analizar las implicaciones que comportan cada una de estas opciones, estudiar cual de entre ellas es más adecuada para el problema planteado y finalmente elegir el método a emplear. Las trayectorias de planeo, clasificadas según su agresividad como superficial, media o pronunciada, son trayectorias de entrada que emplean planeos para atravesar la atmósfera y controlar su descenso. Generalmente implica realizar un descenso con ángulo de inclinación lateral constante. Este tipo de trayectorias suelen encontrarse limitadas en su alcance, aunque son más sencillas de realizar a nivel práctico. Las trayectorias balísticas son entradas muy agresivas, que evitan los vuelos prolongados y que buscan llevar la cápsula a tierra lo antes posible. Esto conlleva que la cápsula sufra una aceleración muy elevada. Es un método con alcance limitado y empleado generalmente solo como medida de emergencia. Este es el caso del sistema de reentrada de emergencia de las cápsulas Soyuz rusas. En particular, este método fue empleado con éxito en la reentrada de la misión Soyuz TMA 10. También son el método más empleados en cargas no vivas, dado que la restricción por aceleración es mucho más permisiva Las trayectorias “Skip-Out” son, posiblemente, las trayectorias más extendidas y conocidas. Implican una primera fase en la que la atmósfera actúa de resorte haciendo rebotar al vehículo que, con el control adecuado, consigue recuperar el rumbo y entrar de nuevo a la atmósfera, habiendo perdido energía cinética tras el rebote. La nave rebota un cierto número de veces hasta que las condiciones operacionales son tales que el vehículo se comporta como un proyectil balístico normal. Estos métodos permiten mayor control, mayor flexibilidad y mayor alcance. Comportan aceleraciones y calores controlables y razonables. A cambio, estas trayectorias requieren de controladores más complejos. En vista de los beneficios tanto operacionales como de seguridad se considera justificado la elección de las trayectorias “Skip-Out”, a pesar de la mayor dificultad de control. 3.2. Etapas de entrada Presentados y discutidos los diversos métodos de reentrada atmosférica y justificada la elección de las trayectorias tipo “Skip-Out” es imperativo estudiar ésta en detalle para poder extraer de ella la información necesaria para construir el algoritmo de control que permita 26 reproducir de manera consistente y adecuada esta trayectoria para el mayor espectro de casos posibles. En particular se estudiará el caso “Skip-Out” de un salto. Es decir, un único rebote. Se realiza este análisis por dos motivos: • Se reduce enormemente la complejidad de programación teniendo en cuenta en el control que se va a plantear • El primer rebote es el más crítico y por tanto el que es más relevante estimar. Toda trayectoria Skip-Out consta de una serie de etapas que la definen. El estudio de estas etapas, así como los preceptos de control a emplear en cada una es una materia ampliamente estudiada tanto a nivel teórico como a partir de datos de carácter empírico [17, 18]. La trayectoria Skip-Out es la trayectoria empleada en misiones como las Apollo (1960-75), Pathfinder (1976), Curiosity (2011) o Perseverance (2020) entre otras. En todas ellas, las siguientes etapas no solo se encuentran presentes, sino que son de vital importancia para el correcto funcionamiento de los sistemas de control. Se presenta a continuación un diagrama de las etapas de una trayectoria Skip-Out: Ilustración 3. 1 Etapas presentes en las trayectorias "Skip-Out" 27 Se presentan también la modulación del ángulo de inclinación lateral del Apollo 4 y Apollo 10, que se emplearán en el diseño del controlador [19, 20]: Ilustración 3. 2 Histórico de variación del ángulo de inclinación lateral en la reentrada de la Apollo10 28 Ilustración 3. 3 Histórico de variación del ángulo de inclinación lateral en la reentrada de la Apollo10 Etapa 0: Descenso inicial Esta etapa es la que comienza una vez el vehículo empieza el descenso, pero todavía no ha acelerado lo suficiente como para poder actuar de manera consecuente sobre el vehículo. Toda la literatura consultada y los datos empíricos observados indican que en esta etapa la acción a tomar es mantener el ángulo de inclinación lateral constante. Si se analizan los valores de estas acciones se puede inferir que se suele optar por un ángulo que o bien es cercano a cero, o es directamente 0. Por tanto, mientras dure esta etapa, se solicitará un ángulo de inclinación lateral de 0º. Etapa 1: Arrastre Constante La siguiente etapa comienza una vez el vehículo toma la suficiente aceleración. Una característica importante de esta etapa es la poca variación del arrastre experimentado por el vehículo, así como contener el primer mínimo de la trayectoria. Éste se da por la cancelación momentaria de la sustentación y el efecto gravitatorio. Esto ocasiona una situación en la que, si el ángulo de incidencia no es muy agresivo o se controla adecuadamente, la sustentación 29 se impondrá a la gravedad. Esta etapa engloba por tanto el mínimo y parte de la trazada anterior. Analizando los datos empíricos, se observa que la actuación sobre el ángulo de inclinación lateral comienza un cierto tiempo tras este mínimo. Por tanto, solo es necesario obtener el mínimo. Durante esta etapa mantendremos el ángulo de inclinación de la etapa anterior. Por tanto 𝜎=0º. Etapa 2: Control de Ascenso La etapa de control de ascenso comienza una vez se pasa el primer mínimo de la trayectoria y representa por una parte el momento más característico de esta trayectoria como también el más crucial. Esta etapa concluye una vez se alcanza el máximo que indica que el “skip” ha concluido y, por tanto, que la gravedad vuelve a dominar a la sustentación. En esta etapa, la imposición por defecto de la sustentación sobre la atracción gravitatoria es de sumo interés. Es importante analizar el tipo de actuaciones que permiten convertir lo que por defecto sería una trayectoria supercircular, es decir, una trayectoria que toma una velocidad y dirección tal, que vuelve a entrar en órbita en una trayectoria de descenso. Observando tanto los datos empíricos como las conclusiones alcanzadas en [18] se pueden extraer una serie de preceptos que se pueden seguir en el diseño del control de esta etapa. • Durante un cierto tiempo al principio de la etapa se debe aplicar un ángulo de inclinación lateral superior a los 100º, y si es posible, de 180º. • Pasado el primer cambio de ángulo hacia valores superiores a 100º se debe aplicar, en un periodo similar un cambio a en torno a -60º. Para cumplir con estos preceptos se opta por plantear un parámetro de diseño conocido en la sintáctica del código como Tp (Tiempo predicción). Este parámetro cuyo ajuste se tratará más adelante, indicará cuánto tiempo se permanecerá en cada una de las dos acciones planteadas. Tp representa un horizonte de predicción empleado para precalcular en los bucles de definición del controlador. Se decide actuar con 𝜎=180º durante dos veces el Tp, desde el comienzo de la etapa. Acabado este periodo se procede a actuar con 𝜎=-60º hasta el fin de la etapa. El fin de la etapa se alcanza una vez se alcanza el siguiente máximo de la trayectoria. Etapa 3: Balística La siguiente etapa a tratar es la etapa balística o tipo Kepler. Durante esta fase, el control suele ser bastante poco eficaz, dado que el vehículo toma altura y la atmósfera ve reducida su influencia. Sin embargo, en este control se va a buscar un skip más suave que el presentado en misiones y trayectorias anteriores al diseño de este controlador. Por ello, la variación de las variables de estado atmosféricas no será tan grande y retendrá influencia. Observando los históricos, se opta por plantear un ángulo constante durante la etapa balística de 𝜎=-70º. La etapa comienza en el máximo de trayectoria y termina una vez transcurre el mismo tiempo que se actuó con -60º en la etapa anterior. 30 Etapa 4: Seguimiento de referencia Se alcanza finalmente la última etapa de la trayectoria “Skip Out”. Esta etapa también se conoce como planeo y, como su nombre indica, implica el planeado del vehículo desde el final del “Skip” hasta alcanzar la altura de despliegue del paracaídas. El objetivo de esta fase es permitir al vehículo alcanzar la posición geográfica óptima para que pueda alcanzar su objetivo de aterrizaje con la mayor exactitud posible. Durante esta etapa se realizan las decisiones de control en función de los históricos y el perfil de velocidades. Se realiza esto por dos motivos: • Esta etapa, en caso de emplearse métodos y programas más robustos y complejos, emplea una referencia y optimizaciones para alcanzar los objetivos de trayectoria punto a punto. Este control no emplea dicha trayectoria dado que, entre otras cosas, pretende otorgar información que facilité la estimación y construcción de dicha referencia • Se observa una posible relación entre la variación temporal de la velocidad y la actuación en esta etapa en los casos empíricos. Por ende, se considera oportuno observar y determinar la existencia o no de esta posible causalidad y, en caso positivo, aplicarla. Se presenta a continuación el perfil de velocidades del Apollo 10 Ilustración 3. 4 Perfil de velocidad de la cápsula Apollo10 durante la reentrada 31 Examinando en conjunto las actuaciones de dicha reentrada se observa un cierto patrón de comportamiento. Cuando la pendiente de la velocidad es elevada, se actúa con ángulos positivos, cuando la pendiente es pequeña, se actúa con ángulos negativos. Esto implica que, pasada la etapa balística, en torno a los 300 segundos, cuando la velocidad parece trazar curvas de 120º pasando de casi horizontal (constante) a rápidamente variante, se actúa con 𝜎=-50º. Cuando la velocidad varía rápidamente, se actúa con 𝜎=90º. Una vez la velocidad vuelve a tornarse constante, en este caso, un poco después del segundo 500, se realiza una última actuación de 0º, con la intención de esperar a alcanzar la altitud de apertura del paracaídas. 3.3. Restricciones y Limitaciones El control debe respetar también una serie de restricciones operacionales para garantizar un funcionamiento adecuado. Éstas deben ser planteadas en el propio controlador a fin de poder obtener vectores de actuación que garanticen que la trayectoria no generará fallos. Es importante entender, sin embargo, que las diversas restricciones afectan al control de diversa manera y como tal deberán ser tratadas de manera conveniente cada una. Estas son las restricciones más importantes a tener en cuenta: • La trayectoria ha de ser “Skip-Out”, por lo que al final de la simulación, el vehículo no puede encontrarse en una trayectoria supercircular • El flujo térmico que recibe el vehículo no puede sobrepasar en ningún momento el límite térmico establecido. • La simulación debe parar una vez la capsula alcanza los 10km de altura • El número de Mach ha de ser observado, pues no debe caer por debajo de 3 para que el modelado mantenga su validez teórica. Para tratar la primera restricción se realiza un primer control de 1 actuación alterando el ángulo de incidencia hasta obtener una trayectoria que garantiza la entrada en la atmosfera Para atender la segunda restricción, primero se resuelve el problema de diseño trabajando con una trayectoria agresiva de convergencia a superficie espontánea. Una vez resuelto se comprueba si se excede el límite térmico en algún momento. De excederse, se vuelve a resolver el problema de diseño, pero empleando ahora un ángulo menos agresivo que el del control ya diseñado como punto de partida. Se repite este proceso iterativo hasta obtener un conjunto de ángulos y actuaciones que no excedan el límite térmico. Para tratar la tercera, se genera código capaz de recortar el último tramo de la trayectoria hasta obtener el primer punto donde se alcanza la altura de apertura del paracaídas. Obtenido este punto, se acorta de igual manera el tiempo asignado a la simulación. Para la última restricción se obtiene la conductividad térmica de la atmosfera dependiente de la altura para Marte, se calcula el Mach en cada punto de la simulación y se obtiene el tiempo donde el Mach cae por debajo de 3. Esta restricción será de las más relevantes durante el diseño. Esto se debe a que se trabajará con la hipótesis de que, para el caso particular de Marte, las ecuaciones dinámicas del modelo mantienen su validez por debajo de la hipersonicidad suficientemente. Sin embargo, esto no es aplicable al modelo térmico, que trabaja con ciertas hipótesis que requieren hipersonicidad de una manera mucho más fuerte. 32 Esto implica que la evaluación térmica solo abarcará el tiempo en el que M>3. Afortunadamente, esto suele incluir con bastante margen el tramo de máximo flujo térmico de la trayectoria, Por tanto, uno de los aspectos que se evaluaran en el modelo es la validez o no de simular con el modelo dinámico por debajo de M=3. 3.4. Variables de Control El control del vehículo se realiza a partir del manejo de las variables manipuladas. El diseño de control debe generar las actuaciones sobre estas variables que posibilitan obtener la trayectoria que mejor cumple los objetivos de la maniobra. Dado el modelo elegido, se cuenta con tres variables manipuladas, de las cuales dos serán manipuladas automáticamente por el controlador. El ángulo de inclinación lateral y el ángulo de trayectoria serán manipulados automáticamente, mientras que el ángulo de rumbo inicial será manipulado manualmente. Para manipular estas variables, el vehículo elegido (Orión Command Capsule) emplea un sistema de 12 reactores RCS (Reaction Control System). Estos reactores permiten generar pulsos intermitentes de empuje que permiten alterar la orientación del vehículo. Este método ha sido empleado en múltiples misiones en el pasado y su validez y viabilidad está demostrada. Dado que la manipulación de las variables se realiza con pulsos de reacción habrá que tener en cuenta dos aspectos: • El coste de maniobrar del vehículo está fuertemente ligado al consumo de combustible que realiza el RCS. Cuantas más manipulaciones se realicen, más combustible será necesario, lo que incrementará el peso del módulo y el coste total. • Dada la naturaleza de la acción de control, el gasto de combustible será mayor cuanto más brusca sea la manipulación. Por ende, se debe prestar atención no solo al número de actuaciones sino también a la magnitud del incremento. Para la manipulación manual, el coste operativo es muy pequeño dado que es una elección que se ha de tomar en el diseño de misión, no un aspecto que el vehículo deba manipular solo. Es por esta razón que no se dota al controlador de la capacidad de manipularlo. Para elegir el ángulo sencillamente se realiza la simulación variando de manera secuencial el ángulo de rumbo inicial, se analiza el rumbo resultante y se obtiene el intervalo de ángulos válidos. De entre estos, se extrae el que se empleará. Se plantea la capacidad de manipulación de las variables. El ángulo de trayectoria es una variable de estado en el modelo termodinámico elegido que representa cómo de agresiva es o está siendo una cierta trayectoria. La actuación sobre este ángulo es, en realidad, el objetivo último de toda la acción de control. Sin embargo, la manipulación directa de este ángulo es muy compleja. Esto se debe a que su variación depende de otros factores a parte del ángulo de inclinación lateral 𝑑𝛾 𝑑𝑡 =𝛾󰇗=𝜌𝐶𝐿𝑆 𝑉 2𝑚 𝑐𝑜𝑠𝜎−𝑔𝑐𝑜𝑠𝛾 𝑉+𝑉𝑐𝑜𝑠𝛾 𝑟. (3.1) Es por ello por lo que solo se podrá actuar directamente sobre el valor inicial del ángulo de trayectoria. 33 Por otra parte, como ya se planteó anteriormente, la manipulación del ángulo de trayectoria se puede realizar de manera indirecta manipulando el ángulo de inclinación lateral. En particular el sistema RCS permite manipular esta variable mid-flight, lo que da pie a que se pueda diseñar el control. Existen varias restricciones de actuación sobre el ángulo de inclinación lateral • El ángulo de inclinación lateral cuenta con un límite inferior y superior. Su rango es: [-180º,180º] • La actuación cuenta con un cierto retardo. Por tanto, es aconsejable que exista un cierto tiempo entre actuaciones dado que en el caso real la actuación no es instantánea. El control debería ser capaz de anticipar la actuación prediciendo el resultado de ésta y aplicando la manipulación con tiempo. 3.5. Arquitectura de programa Elegido el tipo de trayectoria, analizadas las etapas de ésta y planteadas las restricciones operacionales, se debe pasar a programar el controlador. Se expone a continuación la implementación tanto de las etapas como de las restricciones. La programación se realiza en MATLAB 2020b. Se genera una función de MATLAB que toma como entradas la latitud y longitud objetivo, el ángulo de trayectoria de vuelo inicial y el ángulo de rumbo elegido por el diseñador. Como salidas se obtendrá el vector de valores iniciales calculados por el programa, el cálculo del calor experimentado por el vehículo y el vector de tiempos que permite graficar este calor. Ilustración 3. 5 Declaración de la función Control El controlador comienza con una fase de inicialización de variables, Se realiza una carga de las constantes vehiculares y planetarias. El manejo de estas constantes de esta manera permitirá manipular de manera completamente independiente al controlador estas variables en el futuro. Esto es especialmente interesante dada nuestra intención de poder aplicar este programa a otros planetas y vehículos. A continuación, se cargan las matrices de modelado atmosférico, con el diseño ya expuesto en el apartado de modelado, correspondientes a la latitud y longitud objetivo. Será necesario cargar la matriz de densidad, de temperatura y de conductividad térmica. Se procede a generar la primera versión de x0, el vector de valores iniciales que empleará el simulador. Para ello, se toman los valores angulares dados a la función como entradas. Estos ángulos son convertidos a radianes. Nótese que se toma Latitud y Longitud objetivo como valor inicial. Esto se debe a que se trata de una primera aproximación. Una vez el controlador este diseñado, se cambiarán estos dos ángulos adecuadamente. Nótese también que el valor inicial de la altura y de la velocidad se encuentran parametrizados para poder ser manipulados cómodamente desde la función de carga de constantes. 40 4. Configuración y ajuste del Simulador Elegido el modelo y el control, se puede proceder a comprobar las diferentes hipótesis planteadas, observar la influencia de los diferentes parámetros de modulación y, finalmente, analizar los resultados otorgados por el simulador y las implicaciones de estos. 4.1. Simulador Primeramente, se debe decidir qué solver emplear para realizar los cálculos. Se valora emplear diferentes solvers de sistemas de ecuaciones diferenciales ordinarias (ODEs) [22, 23]. El objetivo es evaluar los tiempos de ejecución del algoritmo de diseño según el solver empleado y escoger la combinación más adecuada. Se evalúa emplear la función predefinida de MATLAB ode45, así como el método Runge-Kutta de orden 4 clásico (rk4) y un método predictor corrector de orden 4 (PECE4). Se emplean las instrucciones TIC-TOC de MATLAB para medir el tiempo de procesamiento tanto del bucle diseño como del bucle de diseño más la ejecución de la simulación final. El código de implementación de estos métodos se puede encontrar en el anexo D. El método ode45 de MATLAB es el método conocido como método de Dormand-Prince. Es un método explícito de paso variable, formado por dos métodos de Runge-Kutta encajados, de órdenes 4 y 5. Tiene 7 etapas, porque los coeficientes del método de orden 5 se eligen para minimizar su parte principal del error, pero que en último término sólo requiere 6 evaluaciones funcionales por paso. Los dos métodos encajados se utilizan para calcular una estimación del error que permite elegir el tamaño de paso en función de la tolerancia solicitada. Finalmente se avanza con el método de orden 5. Estos coeficientes se presentan en el siguiente tablero de Butcher. Ilustración 4. 1 Tablero de Butcher Método Dormand-Prince Este método resulta indicado para problemas que no presentan una gran rigidez. 41 Se realizan 10 mediciones para cada método. Se obtienen los siguientes resultados: ode45 rk4 PECE4 1 3.332843 18.588857 10.544599 2 3.386435 19.757945 10.294613 3 3.283669 18.998639 10.329623 4 3.384126 19.202077 10.360399 5 3.321538 20.421473 10.385827 6 3.378544 18.887958 10.597939 7 3.312898 19.019099 10.824934 8 3.344411 18.543586 10.601495 9 3.332625 18.625486 10.627487 10 3.293297 18.945975 10.540889 Tabla 4. 1 Comparativa de tiempos de ejecución de solver Existe una clara diferencia en los tiempos de ejecución del programa con uno u otro método. Sin embargo, antes de tomar una decisión, es importante evaluar la respuesta que cada método proporciona. Es importante que la ejecución sea ágil, pero eso no debe poner en peligro la validez de los resultados. Por ello, se simula una trayectoria empleando cada uno de los métodos y se procede a comparar los resultados: 42 Ilustración 4. 2 Comparativa de perfiles de altitud creados por diferentes solvers numéricos Se puede observar una gran similaridad entre las soluciones obtenidas con cada método. En particular los primeros 100 segundos, la zona más crítica de la simulación, otorga resultados casi idénticos para los tres métodos. Por ello, y teniendo en cuenta la mayor rapidez y similar resultado, el método que se empleará en las simulaciones será el método ode45 preprogramado en MATLAB. A mayores, se conoce que el sistema de ecuaciones diferenciales no presenta problemas de rigidez, lo que elimina la única restricción del método. Finalmente, el entorno de simulación empleado será MATLAB R2020b. 4.2. Modulación. Evaluación de parámetros Elegido el método, el simulador está preparado para ser utilizado. Sin embargo, a pesar de estar preparado, no está ajustado y modulado. Como se ha destacado repetidas veces durante este trabajo, varios parámetros cuentan por defecto con ciertos valores que deben ser modulados y cuya relevancia e influencia debe ser estudiada para determinar si las hipótesis de diseño que les conciernen son válidas, así como para encontrar los valores más adecuados para estos parámetros. 43 Se empieza analizando la influencia de los ángulos empleados en el algoritmo de control para determinar la criticalidad de cada actuación. Ilustración 4. 3 Influencia de la primera actuación de control en la trayectoria En esta gráfica se puede observar que la primera etapa, como se predijo, no es controlable. Los primeros 100 segundos de la simulación, que corresponden a esta etapa, no se ven alterados por el valor del ángulo de inclinación lateral. Esto se debe a la poca importancia de la atmósfera a gran altura, dado que la influencia del ángulo de inclinación lateral está ligada a la densidad atmosférica. Sin embargo, hacia el final de la etapa, la densidad empieza a aumentar considerablemente. Por ello, la elección del ángulo constante de la etapa 0 influye en la trayectoria final, dado que el ángulo altera el punto en el que comienza el rebote. Se decide mantener el ángulo de esta etapa en 0º, dado que resulta en un rebote más suave. Sin embargo, si se desea un rebote más alto y mayor alcance se considera emplear 135º. Se elegiría 135º sobre 180º por la mayor facilidad de controlar el primero al ser un valor menos extremo. 44 Ilustración 4. 4 Influencia de la segunda actuación de control en la trayectoria En esta gráfica se puede apreciar el efecto que tiene sobre la simulación la elección de la segunda actuación de control. Como se puede apreciar, esta elección es crítica. El ángulo correcto otorga resultados válidos, mientras que el ángulo incorrecto otorga trayectorias supercircular. Por tanto, la elección de cualquier ángulo en el intervalo (-90º, 90º) es de entrada, no válido. A mayores, la elección de cualquier ángulo entre 90º y 135º comporta bastante riesgo de, en ciertas situaciones, causar supercircularidad. Por tanto, se opta por mantener la actuación en 180º. 45 Ilustración 4. 5 Influencia de la tercera actuación de control en la trayectoria Se analiza ahora la influencia de la tercera actuación en el ángulo. Como se puede observar en la gráfica, este ángulo tiene una influencia pronunciada en el resultado final. La trazada realizada por la trayectoria con ángulo 0º se parece mucho a la referencia buscada, mientras que el valor de diseño que se había planteado resulta un poco insuficiente. Es por ello por lo que se opta por modificar de esta actuación y pasar de -60º a 0º. 46 Ilustración 4. 6 Influencia de la cuarta actuación de control y posteriores en la trayectoria Para el resto de los ángulos no se aprecia diferencia entre unos y otros ángulos. Esto se debe a que el ángulo de inclinación lateral afecta al arrastre y la sustentación. Estos a su vez dependen directamente de la velocidad. Pasado un cierto instante, la velocidad es suficientemente pequeña como para que no haya cambios sustanciales en la trayectoria. Se analiza también la influencia de los ángulos en el rumbo. 47 Ilustración 4. 7 Influencia de la segunda actuación de control en la longitud Ilustración 4. 8 Influencia de la segunda actuación de control en la latitud Ilustración 4. 9 Influencia de la tercera actuación de control en la longitud 48 Ilustración 4. 10 Influencia de la tercera actuación de control en la latitud Ilustración 4. 11 Influencia de la cuarta actuación de control y posteriores en la longitud Ilustración 4. 12 Influencia de la cuarta actuación de control y posteriores en la latitud 49 Observando la influencia de los angulos de la segunda y tercera actuación, así como la invariabilidad en el resto de las actuaciones, se extrae que la elección de ángulos propuesta anteriormente no causa problemas de posicionamiento que invalidarían la decisión. Modulada la acción de control se debe modular el parámetro Tp, que afecta a la generación de los intervalos de actuación. Por ello, se simula una trayectoria teniendo en cuenta diferentes tiempos de predicción. Ilustración 4. 13 Influencia en la trayectoria del parámetro de modulación Tp 56 Ilustración 4. 20 Influencia de la velocidad inicial en el perfil de velocidad Observando estas gráficas, se comprueba que la elección de los valores de la Perseverance da buenos resultados. Por tanto, será la empleada en los diferentes casos. De esta manera, quedan elegidos los parámetros de diseño y modulación del programa. • Tp=20 • σ= [0,180,0,-70,0,90,0] 57 58 5. Resultados de Simulación 5.1. Caso 1 Habiendo elegido los parámetros y ultimado el simulador, es necesario ahora emplear el simulador para extraer resultados y evaluar su utilidad. Se estudia primeramente el caso de prueba que se ha empleado a lo largo de todo el proceso de diseño. Los datos operacionales son: Longitud Objetivo: 174.24384º Este Latitud Objetivo: 23. 51457º Sur En este caso, las coordenadas corresponden a las coordenadas de entrada de la misión Apollo 10 en la Tierra. En otras situaciones se emplearán coordenadas y ángulos descritas con números enteros. Se ejecuta el diseñador de control empleando como altura de interfase 130000 metros y velocidad inicial 6000 m/s Se pretende emplear el ángulo que se empleó en la misión Apollo 10. Sin embargo, es necesario comprobar que este ángulo no da problemas al algoritmo de selección de coordenadas. Ángulo de rumbo inicial: 18.0683º. 59 Ilustración 5. 1 Influencia del ángulo de rumbo inicial en el perfil de longitud en el Caso 1 Ilustración 5. 2 Influencia del ángulo de rumbo inicial en el perfil de latitud en el Caso 1 60 Observando estas gráficas, las coordenadas objetivo no dan problemas al algoritmo, por lo que no hay ningún problema por emplear el ángulo de rumbo inicial propuesto. Ejecutando el programa de diseño del controlador se obtiene la siguiente política de actuación: Ilustración 5. 3 Perfil de variación del ángulo de inclinación lateral en el Caso 1 Y el siguiente vector de valores iniciales: [3519000, 2.658711, 0.556245976, 6000, -0.200712898, 0.315351325]. Se pasa a simular la entrada atmosférica empleando este controlador y vector de valores iniciales. Se obtienen las siguientes gráficas: 61 Ilustración 5. 4 Trayectoria de entrada 3D en el Caso 1 62 Ilustración 5. 5 Perfil de altura en el Caso 1 63 Ilustración 5. 6 Perfil de velocidad en el Caso 1 64 Ilustración 5. 7 Perfil de longitud en el Caso 1 65 Ilustración 5. 8 Perfil de latitud en el Caso 1 72 Ilustración 5. 15 Trayectoria de entrada 3D en el Caso 2 73 Ilustración 5. 16 Perfil de altura en el Caso 2 74 Ilustración 5. 17 Perfil de velocidad en el Caso 2 75 Ilustración 5. 18 Perfil de longitud en el Caso 2 76 Ilustración 5. 19 Perfil de latitud en el Caso 2 77 Ilustración 5. 20 Perfil del ángulo de rumbo en el Caso 2 78 Ilustración 5. 21 Perfil del ángulo de trayectoria en el Caso 2 79 Ilustración 5. 22 Flujo térmico al vehículo en el Caso 2 La posición final según la simulación, a 10 km de la superficie es 87.2844ºNorte y 24.9982ºEste. El calor máximo soportado 234.14 KW/m2. De nuevo, esto otorga un enorme margen a nivel térmico, confirmando la necesidad de reducir el escudo térmico del vehículo. Por otra parte, existe un desvío con respecto al objetivo de 0.28º en dirección Norte y 0.001º en dirección Este. Esto comporta una desviación aproximada de 16.6 km. Este cálculo hace las mismas suposiciones que el caso anterior. 5.3. Caso 3 Se evalúa una tercera posición para observar cómo se desenvuelve el programa con diferentes coordenadas y confirmar que funciona correctamente. En este caso, los datos operacionales son: Longitud Objetivo: 215º Este Latitud Objetivo: 35º Norte 80 Se emplea 130km como altura de interfase atmosférica y 6000 m/s como la velocidad inicial. La latitud objetivo localiza la maniobra cerca del ecuador, esto implica que algoritmo no debería tener problemas para encontrar las coordenadas validas, para cualquier ángulo de rumbo inicial. Se comprueba la influencia de este ángulo en la longitud y latitud para confirmarlo. Ilustración 5. 23 Influencia del ángulo de rumbo inicial en el perfil de longitud en el Caso 3 81 Ilustración 5. 24 Influencia del ángulo de rumbo inicial en el perfil de latitud en el Caso 3 Analizando estas gráficas, se observa que como se suponía, la elección del ángulo de rumbo inicial no da grandes problemas al algoritmo de selección de coordenadas. Por tanto, se elige un ángulo que en vista de la influencia que se acaba de observar de buenos resultados. Ángulo de rumbo inicial: 40º 88 Ilustración 5. 31 Perfil del ángulo de rumbo en el Caso 3 89 Ilustración 5. 32 Perfil del ángulo de trayectoria en el Caso 3 90 Ilustración 5. 33 Flujo térmico al vehículo en el Caso 3 La posición final según la simulación, a 10 km de la superficie es 35.3801º Norte y 213.9863º Este. El calor máximo soportado 280.57 KW/m2. Por otra parte, existe un desvío con respecto al objetivo de 0.35º en dirección Norte y 1.1º en dirección Oeste. Esto comporta una desviación aproximada de 69 km. Este cálculo asume la aproximación de la distancia al objetivo como el arco de circunferencia para el error angular, con radio igual al radio planetario más 10 km. 91 92 6. Conclusiones, Repercusiones y líneas futuras 6.1. Conclusiones Extraídos los resultados y construido el programa, se debe evaluar el cumplimiento de los objetivos del trabajo. • La elección del MCD como modelo atmosférico, el sistema de 6 variables de Allen-Egger el modelo térmico de los mismos autores para el modelo termodinámico y los parámetros técnicos ofrecidos por NASA son capaces de representar de manera fidedigna las condiciones de entrada atmosférica y cuentan con el rigor, validez y grado de actualización como para ser relevantes en la actualidad del sector. • El sistema de control diseñado es capaz de generar actuaciones, sin requerir optimizadores ni herramientas matemáticas complejas, que permiten realizar entradas exitosas que cumplen razonablemente con los objetivos y restricciones operacionales. • Se ha comprobado que el simulador es capaz de implementar correctamente la actuación diseñada y que esta causa la respuesta deseada y esperada. • El programa es capaz de generar coordenadas de entrada atmosférica que permiten alcanzar el objetivo, a falta de 10km, con un margen de error menor de 100km, lo que se considera suficiente como primera aproximación. • El programa genera un vector de actuación válido que el vehículo puede emplear para realizar la maniobra y entrega otros datos de interés como el flujo térmico máximo soportado por el vehículo. • El programa ha sido desarrollado tomando el planeta objetivo en consideración, empleando un modelo atmosférico de vanguardia y un vehículo de reentrada que se pretende emplear en la próxima década en el planeta Marte. 6.2. Repercusiones El programa creado en este trabajo tiene su razón de ser en el actual interés cultural, científico, político y económico de la exploración espacial. El crecimiento acelerado del sector ha generado una importante demanda de herramientas científicas y técnicas para democratizar el sector. Este trabajo representa en cierto modo esta evolución de la situación. Empleado como herramienta de diseño, permite generar estimaciones dimensionales y operacionales con muy poca rigidez. Esto implica, que se pueden hacer estimaciones groseras y obtener una respuesta razonable que permita descartar situaciones y configuraciones no válidas sin requerir de simuladores más potentes. A mayores, el modelado y el controlador sencillo permiten que se pueda abordar el problema de la entrada con los conocimientos mínimos sobre la materia y encontrar soluciones razonables. Esto permite que puedan introducirse en el sector con mayor facilidad investigadores, estudiantes, etc. sin requerir primero una formación extensa en mecánica de fluidos, diseño de control, programación, etc. Esta herramienta hace el problema más accesible y abordable. Adicionalmente, al estar desarrollada para Marte en particular, esta herramienta puede ayudar a abordar el diseño de misiones durante las próximas décadas e incrementar el interés general, sobre todo entre estudiantes, en el sector y, en particular, en la exploración de Marte. 93 6.3. Líneas Futuras Si bien se ha buscado generar un trabajo y un programa lo más coherente y completo posible, ningún trabajo está nunca realmente completo. Siempre tiene cabida la mejora y la expansión. Estos son algunos de los puntos de los puntos sobre los que se puede mejorar o profundizar. La programación ha sido desarrollada con el objetivo de que diera resultados. El código generado no está optimizado y la eficiencia del simulador en su conjunto podría beneficiarse enormemente de la visión de la ingeniería informática. Este programa ha sido desarrollado con el planeta Marte en mente, sin embargo, las misiones de las próximas décadas tienen otros destinos aparte de Marte. Europa, la Luna o Titán son algunos ejemplos. Este programa sería fácilmente adaptable a estos escenarios, requiriendo sobre todo modificar los modelos atmosféricos de cada escenario y el vehículo empleado. Este trabajo se ha centrado en el desarrollo de un programa que controlara la maniobra de entrada sin utilizar el dominio S o el dominio Z empleados en la Ingeniería de Control Automático y sin usar controladores basados en respuesta óptima como el MPC y sus derivados. Sin embargo, el desarrollo de programas que los empleen es también de un elevado interés en el sector. La integración de este programa de estimación con controles, lineales o no lineales, es posiblemente el paso lógico a dar si se desea ampliar las capacidades de este programa. 94 7. Bibliografía 7.1. Lista de referencias [1] Morgan Stanley, «Morgan Stanley» [En línea]. Available: https://www.morganstanley.com/ideas/investing-in-space. [Último acceso: 11 Julio 2021]. [2] JPL,NASA, «Jet Propulsion Laboratory» [En línea]. Available: https://www.jpl.nasa.gov/missions. [Último acceso: 8 Julio 2021]. [3] ESA, «Exploration ESA» 12 Marzo 2020. [En línea]. Available: https://exploration.esa.int/web/mars/-/48088-mission-overview. [Último acceso: 11 Julio 2021]. [4] NASA, «NASA’s Lunar Exploration Program Overview» September 2020. [En línea]. Available: https://www.nasa.gov/specials/artemis/. [5] Wolf, Aron & Graves, Claude & Powell, Richard & Johnson, Wyatt., «Systems for pinpoint landing at Mars.» Advances in the Astronautical Sciences., vol. 119, pp. 2-3, 2005. [6] NASA, «NASA» 13 Mayo 2021. [En línea]. Available: https://www.grc.nasa.gov/www/k- 12/airplane/atmosmrm.html. [Último acceso: 10 Julio 2021]. [7] Jules et al., «Atmospheric Environments for Entry, Descent and Landing (EDL)» p. 21, 2007. [8] Forget et al., «Improved general circulation models of the Martian atmosphere from the surface to above 80 km,» Journal of Geophysical Research, vol. 104, nº 10, pp. 24155-24176, 1999. [9] Millour e. Al., «The Mars Climate Database 5.2» 2018. [10] Mars Climate Database, «Mars Climate Databse» [En línea]. Available: http://wwwmars.lmd.jussieu.fr/mars/access.html. [Último acceso: 11 Julio 2021]. [11] NASA, «NASA» [En línea]. Available: https://nssdc.gsfc.nasa.gov/planetary/factsheet/marsfact.html. [Último acceso: 11 Julio 2021]. [12] Derrick G. Tetzman, «Simulation and Optimization of Spacecraft Re-entry Trajectories» 2010. [13] Hicks, Lt. Col. Kerry D., Introduction To Astrodynamic Reentry, 2009. [14] Bibb, Karen & Alter, Stephen & Mcdaniel, Ryan, « Aerodynamic Analysis of Simulated Heat Shield Recession for the Orion Command Module. 10.2514/6.2008-356.» 46th AIAA Aerosciences Meeting and Exhibit, 2008. 95 [15] Stanley ,Doug & As., «NASA's Exploration Architecture Study» 2005. [16] Murphy, Kelly & Brauckmann, Greg & Bell, James & Rhode, Matthew & Owens, D. & Wilson, Thomas & Bibb, Karen & Chan, David & Walker, Eric., «Orion Crew Module Aerodynamic Testing. 10.2514/6.2011-3502.» de 29th AIAA Applied Aerodynamics Conference, Honolulu, 2011. [17] Martin, Dillon, «Atmospheric Entry», University of central Florida, 2017. [18] Crawford, D.J. Lickly H.R. Morth B.S., Apollo Reentry Guidance, 1963. [19] D. Szelc, APOLLO 10 Mission Report. Entry PostFlight Analysis, Houston, Texas, 1969. [20] Hillje, Ernest R., ENTRY AERODYNAMICS AT LUNAR RETURN CONDITIONS OBTAINED FROM THE FLIGHT OF APOLLO 4, Houston, Texas, 1969. [21] NASA, «NASA» 13 Mayo 2021. [En línea]. Available: https://www.grc.nasa.gov/www/k- 12/airplane/atmosmet.html. [Último acceso: 11 Julio 2021]. [22] J. Lambert, Numerical Methods for Ordinary Differential Systems, 1997. [23] E. Hairer, S. Norset y G. Wanner, Solving Ordinary Differential Equations I and II, Springer- Verlag, 1991. [24] NASA, «NASA MArs 2020» Febrero 2020. [En línea]. Available: https://eyes.nasa.gov/apps/mars2020/#/home. [Último acceso: 2021 Julio 11]. [25] E. Mooij, «Characteristic Motion of Re-entry Vehicles» de AIAA Atmospheric Flight Mechanics (AFM) Conference, 2013. [26] Guerreiro, Luis André Carraça, Development of a Guidance and Control Design Tool for Entry Space Vehicles with Different Lift-over-Draft Ratios, 2011. 7.2. Comentario Bibliográfico [4] Es un informe sobre el planteamiento del conjunto de misiones Artemis. La lectura de este informe puede aportar una perspectiva general del planteamiento a futuro de la NASA para misiones tripuladas y puede ser un buen punto de partida en trabajos enfocados en ese tema. [5] Es un artículo que explora los diferentes sistemas que permitirían realizar aterrizajes exactos y precisos en Marte. Se recomienda si se desea profundizar en la mejora de sistemas de entrada para misiones tripuladas. 96 [8, 9, 10]Corresponden al MCD. Este modelo es muy extenso y completo. Es perfecto para el estudio de problemas de entrada en Marte. Sobre todo, si no se tiene acceso al MARS-GRAM de NASA. A mayores, el MCD cuenta con muchas herramientas y datos, de interés en la resolución de otros problemas de ingeniería espacial y física. Se recomienda esta herramienta para cualquier trabajo que requiera de un modelo de las condiciones atmosféricas de Marte. Finalmente, la organización al cargo del MCD también presenta conferencias de renombre internacional sobre la materia de la atmósfera marciana, y su comunidad científica puede ser de gran apoyo en trabajos relacionados con la materia. [12] Esta tesis doctoral es la inspiración de este trabajo, y su lectura se recomienda encarecidamente si desea proseguir profundizando y trabajando sobre el tema. El texto resulta muy claro y accesible y puede actuar tanto de guía de trabajo como de entrada al campo de la entrada atmosférica. [13] Este libro de las fuerzas aéreas americanas ha servido de principal referencia en el desarrollo de este trabajo. Su lectura está muy recomendada si pretende abordar cualquier trabajo relacionado con las entradas atmosféricas. Es un texto muy completo, que presenta modelos y soluciones particulares con distintos grados de complejidad, facilitando la elección del modelo y permitiendo implementar mejoras con facilidad. [14, 15, 16] Estos tres documentos fueron empelados para extraer la información técnica del vehículo. En particular [15] es un informe completo sobre la arquitectura operacional de un cierto sector de NASA. Por ello se recomienda solo atender al apartado 5, en particular lo concerniente a la aplicación de la ORION a entradas extraterrestres. Sin embargo, también se recomienda atender a las conclusiones que se plantean en este apartado sobre política de control en las entradas atmosféricas. [17] Esta tesis fue empleada para obtener nociones sobre los distintos planteamientos disponibles para el control en la entrada atmosférica. En particular presenta un enfoque basado en redes neuronales. El planteamiento puede resultar muy interesante si se consigue integrar exitosamente dichas redes en un planteamiento sencillo como el buscado en este trabajo. [18] Este informe técnico fue empleado para observar los planteamientos de control realizados en las misiones Apollo. Este planteamiento es interesante si se desea generar controladores que no estén basados fuertemente en controlo predictivo. [19, 20] Estos dos informes fueron empleados para extraer conclusiones sobre sus datos empíricos y corroborar conclusiones teóricas. Por desgracia no hay datos accesibles para realizar estas deducciones con misiones a Marte. Sin embargo, se recomienda encarecidamente que en caso de comenzar un trabajo en este campo se trate de obtener dicha información, dado que existen ciertas diferencias en el control particular de la entrada a cada cuerpo celeste. [22, 23] Estos dos libros han sido la referencia principal en materia de resolución numérica del modelo y de su aplicación en el simulador. [25, 26] Estos dos trabajos presentan diferentes enfoques para el abordaje de la entrada atmosférica desde un enfoque más centrado en el control automático. Si se desea profundizar por este camino se recomienda su lectura como punto de partida. 97 104 h=yy(end,1)-r0; if(h<0) h=0; end h_index=round(h/paso_h2)+1; if (h_index >34) h_index=34; end h_tramo=h-d_par2(h_index,4); T=t_par2(h_index,1)+h_tramo*t_par2(h_index,3); heat=h_par2(h_index,1)+h_tramo*h_par2(h_index,3); a=sqrt(heat*T*8.314); M=yy(end,4)/a; if(M<=3&&flag7==0) tq=i+Tp; flag7=1; end end % bucle de calor [ttq,yyq]=ode45(@System,[1 tq],x00); for j=1:1:size(yyq,1) h=yyq(j,1)-r0; if(h<0) h=0; end h_index=round(h/paso_h2)+1; if (h_index >34) h_index=34; end h_tramo=h-d_par2(h_index,4); d=d_par2(h_index,1)+h_tramo*d_par2(h_index,3); q(j,1)=(C*yyq(j,4)^3*sqrt(d))/(sqrt(Rn)); end if(flag6==0) x000=x00; flag6=1; end if(max(q)>qlim ) x0=x000; x0(5,:)=x0(5,:)+0.00174533; x0(5,:) x000=x0; flag2=1; flag3=0; flag4=0; 105 flag5=0; flag7=0; flag_control=2; else flag_prin=0; end i=1; end [tq,yyq]=ode45(@System,[1 tq],x000); for j=1:1:size(yyq,1) h=yyq(j,1)-r0; if(h<0) h=0; end h_index=round(h/paso_h2)+1; if (h_index >34) h_index=34; end h_tramo=h-d_par2(h_index,4); d=d_par2(h_index,1)+h_tramo*d_par2(h_index,3); q(j,1)=(C*yyq(j,4)^3*sqrt(d))/(sqrt(Rn)); end [t,y]=ode45(@System, [1 tf],x000); for k=0:1:size(y,1)-1 h_limit=(y(size(y,1)-k)-r0); if(h_limit>10000) t_limit=(t(size(t,1)-k)); tf=t_limit break; end end [t,y]=ode45(@System, [1 tf],x000); plot(t,y(:,2)),xlabel('tiempo(s)'),ylabel('longitud(rad)' ) plot(t,y(:,3)),xlabel('tiempo(s)'),ylabel('latitud(rad)') varlon=y(end,2)-y(1,2); varlat=y(end,3)-y(1,3); lon0=y(1,2)-varlon; lat0=y(1,3)-varlat; x000(2,:)=lon0; x000(3,:)=lat0; x=x000; 106 ANEXO C: Código de simulación clear all % DATADROP global m; global Cl; global Cd; global S; global Rn; global r0; global g0; global rho0; global C; global V0; global tcon; global rear; % simulador global tf global flag_control; global sigma; flag_control=0; % Carga de Constantes f=feval(@DataDrop) % Carga de objetivos lon=25; lat=87; gamma_0=-5; azi=270; % Generación de controlador [x0,tq,q]=Control(lon,lat,gamma_0, azi); tf=round(tf); [t,y]=ode45(@System, [1 tf],x0); ts=linspace(1,tf,tf); for s=1:1:size(ts,2) tsig=ts(1,s); if(tsig<tcon(1,1)) sig=sigma(1,1); elseif(tsig>=tcon(1,1) && tsig< tcon(1,2)) sig=sigma(1,2); elseif(tsig>=tcon(1,2) && tsig< tcon(1,3)) sig=sigma(1,3); elseif(tsig>=tcon(1,3) && tsig< tcon(1,4)) sig=sigma(1,4); 107 elseif(tsig>=tcon(1,4) && tsig< tcon(1,5)) sig=sigma(1,5); elseif(tsig>=tcon(1,5) && tsig< tcon(1,6)) sig=sigma(1,6); elseif(tsig>=tcon(1,6)) sig=sigma(1,7); end act(1,s)=sig; end plot(ts,act),xlabel('tiempo(s)'),ylabel('angulo inclinación lateral(º)') plot(tq,q),xlabel('tiempo(s)'),ylabel('flujo térmico(w/m2)') max(q) subplot(6,1,1),plot(t,y(:,1)- r0),xlabel('tiempo(s)'),ylabel('altura(m)') subplot(6,1,2),plot(t,y(:,4)),xlabel('tiempo(s)'),ylabel( 'velocidad(m/s)') subplot(6,1,3),plot(t,y(:,2)),xlabel('tiempo(s)'),ylabel( 'longitud(rad)') subplot(6,1,4),plot(t,y(:,3)),xlabel('tiempo(s)'),ylabel( 'latitud(rad)') subplot(6,1,5),plot(t,y(:,6)),xlabel('tiempo(s)'),ylabel( 'ángulo de rumbo(rad)') subplot(6,1,6),plot(t,y(:,5)),xlabel('tiempo(s)'),ylabel( 'ángulo de trayectoria(rad)') xt=y(:,1).*cos(y(:,3)).*cos(y(:,2)); yt=y(:,1).*cos(y(:,3)).*sin(y(:,2)); zt=y(:,1).*sin(y(:,3)); [x,y,z]=sphere(500); x=x*3389000; y=y*3389000; z=z*3389000; figure(); plot3(x,y,z,xt,yt,zt),xlabel('(m)'),ylabel('(m)'),zlabel( '(m)') 108 ANEXO D: Código de implementación de métodos numéricos Método Runge-Kutta de orden 4 (clásico) function[t,y]=rk4(f,t0,tfinal,y0,M) % f es la función de la ecuación diferencial % t0 y tfinal son los extremos del intervalo % y0 es la condición inicial(t0) % M es el número de pasos h=((tfinal-t0)/M); t=t0:h:tfinal; y=zeros(length(y0),M+1); y(:,1)=y0; for i=1:M k1=feval(f,t(i),y(:,i)); k2=feval(f,t(i)+(h/2),y(:,i)+(h/2)*k1); k3=feval(f,t(i)+(h/2),y(:,i)+(h/2)*k2); k4=feval(f,t(i)+h,y(:,i)+h*k3); y(:,i+1)=y(:,i)+((h/6)*(k1+2*k2+2*k3+k4)); end t=t'; y=y'; Método Predictor Corrector orden 4 function [t,y]=PECE4(f,t0,tfinal,y0,M) %Entrada - f es la funcion de la ecuacion diferencial % - t0 y tfinal son los extremos del intervalo % - y0 es la condicion inicial y(t0) % - M es el numero de pasos = (tfinal-t0)/h %Salida - [t,y] donde t es el vector de abscisas e % y es el vector de ordenadas (ambos se devuelven como % vectores columna) h=(tfinal-t0)/M; t=t0:h:tfinal; y=zeros(length(y0),M+1); y(:,1)=y0; [tt,yy]= rk4(f,t(1),t(4),y0,3); y(:,1:4)=yy'; for n=1:4 ff(:,n)=feval(f,t(n),y(:,n)); end 109 for i=1:M-3 y(:,i+4)=y(:,i+3)+(h/24)*(55*ff(:,4)- 59*ff(:,3)+37*ff(:,2)-9*ff(:,1)); ff(:,1:3)=ff(:,2:4); ff(:,4)=feval(f,t(i+4),y(:,i+4)); y(:,i+4)=y(:,i+3)+(h/24)*(9*ff(:,4)+19*ff(:,3)- 5*ff(:,2)+ff(:,1)); ff(:,4)=feval(f,t(i+4),y(:,i+4)); end t=t'; y=y';