Modelado y diseño del sistema de control de reentrada en la atmósfera de un vehículo de lanzamiento

Proyecto Fin de Carrera Ingeniería de Telecomunicación Formato de Publicación de la Escuela Técnica Superior de Ingeniería Autor: F. Javier Payán Somet Tutor: Juan José Murillo Fuentes Dep. Teoría de la Señal y Comunicaciones Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2013 Proyecto Fin de Grado Ingeniería Aeroespacial Modelado y diseño del sistema de control de reentrada en la atmósfera de un vehículo de lanzamiento Autor: Alberto López Platero Tutor: Carlos Bordons Alba Dpto. Ingeniería de Sistemas y Automática Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2022 Proyecto Fin de Grado Ingeniería Aeroespacial Modelado y diseño del sistema de control de reentrada en la atmósfera de un vehículo de lanzamiento Autor: Alberto López Platero Tutor: Carlos Bordons Alba Catedrático de Universidad Dpto. Ingeniería de Sistemas y Automática Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2022 Proyecto Fin de Grado: Modelado y diseño del sistema de control de reentrada en la atmósfera de un vehículo de lanzamiento Autor: Alberto López Platero Tutor: Carlos Bordons Alba El tribunal nombrado para juzgar el trabajo arriba indicado, compuesto por los siguientes profesores: Presidente: Vocal/es: Secretario: acuerdan otorgarle la calificación de: El Secretario del Tribunal Fecha: Agradecimientos La llegada al final de un camino en ocasiones pedregoso, arduo y empinado es la conclusión de un peregrinaje hacia el conocimiento y la felicidad. Este hito, similar a coronar una montaña de más de 4000 metros, no es más que una etapa en el viaje de la vida. Pese a la persecución de objetivos que se antojan alcanzables únicamente por deidades, siempre he tenido una filosofía clara: «si tienes un ojo puesto en el destino que esperas alcanzar, sólo te queda otro para que te guíe en el viaje». En primer lugar, agradecer a mis padres Antonio y Loreto y a mi hermana Elena, que han sido un elemento fundamental para mantener ese ojo en el camino y disfrutarlo, por haberme ayudado a luchar contra las adversidades que me he encontrado y celebrar las pequeñas victorias que he conseguido. A mis amigos más cercanos, por hacerme feliz y mantenerme motivado. En el ámbito académico, agradecer la confianza puesta en mí desde el colegio, pasando por el instituto, hasta la universidad. Las instituciones públicas educativas han conseguido el reto más complicado al que se enfrentan, la formación a nivel de conocimiento y de raciocinio. También, agradecer a Carlos por su implicación en este trabajo, sin la cual no hubiera salido adelante. Por último, desarrollar este trabajo me ha permitido progresar personalmente gracias a la búsqueda constante de una versión mejor de mí, a la disciplina a lo largo de estos años y a la motivación que me han permitido seguir adelante aún cuando no se vislumbraba atisbo alguno de luz en la más profunda inmensidad. Gracias a todos los que han participado en esta etapa del viaje, sin vosotros no hubiera sido posible. Alberto López Platero Sevilla, 2022 I Resumen La irrupción en el sector espacial del concepto de sistemas de lanzamiento reutilizables ha conseguido hacer del espacio un lugar más accesible. Este trabajo se centra en la maniobra de descenso para la recuperación de un vehículo de lanzamiento. En él, se analizan las características del vehículo Falcon 9 en su descenso a tierra mediante el modelado aerodinámico, atmosférico y propulsivo. Con la finalidad de completar un aterrizaje satisfactorio, una vez validado el modelo, se diseñan los sistemas de control necesarios, el control de empuje y el de actitud. Por último se exploran alternativas para optimizar y mejorar tanto modelos como controladores. III XÍndice 4.1.1 Sistema Inercial Topocéntrico 19 4.1.2 Sistema de ejes horizonte local 19 4.1.3 Sistema de ejes viento 19 4.1.4 Sistema de ejes cuerpo 21 4.2 Ecuaciones del movimiento 23 4.2.1 Ecuaciones dinámicas 23 Tensores de inercia 24 4.2.2 Ecuaciones cinemáticas de la posición 25 4.2.3 Ecuaciones cinemáticas de la actitud 25 Ecuaciones cinemáticas de la actitud con ángulos de Euler 25 Ecuaciones cinemáticas de la actitud con cuaterniones 26 4.2.4 Definición de los dispositivos de control. Fuerzas y momentos 28 Hipótesis 30 Ecuaciones del sistema 30 5 Modelado en Simulink 33 5.1 Estructura 33 5.1.1 Fuerzas y Torques 33 5.1.2 Ecuaciones del movimiento 34 5.1.3 Posición y actitud 35 5.2 Comportamiento del sistema 36 5.2.1 Caída libre vertical 37 5.2.2 Encendido de planta propulsora 38 5.2.3 Deflexión de superficie de control 1 39 5.2.4 Deflexión de superficie de control 4 41 5.2.5 Acoplamiento en actitud 43 6 Sistema de control del empuje 45 6.1 Objetivo 45 6.2 Estructura 45 6.3 Resultados 46 7 Sistema de control de la actitud 49 7.1 Objetivo 49 7.2 Estructura 49 7.2.1 Diseño del controlador para el ángulo de asiento 50 7.2.2 Diseño del controlador para el ángulo de trayectoria 52 7.2.3 Diseño de trayectoria para el aterrizaje 54 7.3 Resultados 55 7.3.1 Simulación dentro del rango operacional 55 7.3.2 Simulación fuera del rango operacional 57 8 Conclusiones 61 8.1 Trabajo futuro 62 Índice de Figuras 63 Índice de Tablas 65 Bibliografía 67 Notación hAltura hp0Altura de referencia para presión en modelo atmosférico hρ0Altura de referencia para densidad en modelo atmosférico ρDensidad del aire ρSL Densidad del aire al nivel del mar pPresión del aire pSL Presión del aire al nivel del mar TTemperatura del aire αTGradiente térmico DResistencia VVelocidad del aire SSuperficie CdCoeficiente de resistencia aerodinámica αÁngulo de ataque M o Mach Número de Mach Re Número de Reynolds CD,bCoeficiente de resistencia de la base CD,fCoeficiente de resistencia del grid fin YFuerza lateral LMomento de balance NMomento de guiñada MMomento de cabeceo lDistancia del grid fin al eje de simetría CY0Coeficiente de fuerza lateral sin ángulo de resbalamiento CYβ Coeficiente de fuerza lateral en función del ángulo de resbalamiento CYδf Coeficiente de fuerza lateral en función del ángulo de deflexión del grid fin Cl0 Coeficiente de momento de balance sin ángulo de resbalamiento ClβÍndice de estabilidad estática lateral Clδf Coeficiente de momento de balance en función del ángulo de deflexión del grid fin Cn0 Coeficiente de momento de guiñada sin ángulo de resbalamiento CnβÍndice de estabilidad estática direccional Cnδf Potencia de control en guiñada Cm0 Coeficiente de momento de cabeceo para ángulo de ataque nulo CmαÍndice de estabilidad estática longitudinal XI XII Notación Cmδf Potencia de control en cabeceo EEmpuje ˙mGasto másico VsVelocidad de los gases de salida de la tobera AsÁrea de salida psPresión de salida paPresión atmosférica T0Temperatura de remanso p0Presión de remanso a0Velocidad del sonido de remanso ρ0Densidad de remanso γCoeficiente de dilatación adiabática αGradiente térmico MsMach de salida γÁngulo de asiento de la velocidad ξÁngulo de guiñada de la velocidad µÁngulo de balance de la velocidad βÁngulo de resbalamiento θÁngulo de asiento ψÁngulo de guiñada φÁngulo de balance  ωVector velocidad angular del vehículo ΩTensor asociado al vector velocidad angular gGravedad uVelocidad en el eje xb vVelocidad en el eje yb wVelocidad en el eje zb pVelocidad angular alrededor del eje xb qVelocidad angular alrededor del eje yb rVelocidad angular alrededor del eje zb  HVector momento angular  TVector de momentos 1 Introducción Desde tiempos inmemoriales, el ser humano ha asociado la bóveda celeste y los astros que en ella se encuentran, con divinidades y dioses siempre rodeados de misticismo y perfección. Son muchas las civilizaciones antiguas que adoraban a los astros. En la mitología egipcia, el dios Ra es el dios del Sol y del origen de la vida, y de la mitología griega surgieron las constelaciones que conocemos hoy en día. Así, el ser humano ha intentando durante milenios darle forma y vida a un mundo totalmente desconocido e inmenso. Con el paso de los siglos y la aportación de maestros como Aristótles, Galileo o Kepler, se fue forjando una concepción cada vez más realista del Universo que nos rodea. Incluso Julio Verne [1], el popular escritor visionario, en 1865, se atrevió a imaginar una expedición hacia los cielos, lugar de las divinidades de las mitologías antiguas, en concreto hacia la Luna en su libro De la Tierra a la Luna. 1.1 Historia de la exploración espacial No fue hasta la Segunda Guerra Mundial, donde un mayor desarrollo tecnológico y armamentístico podían inclinar la balanza hacia la victoria, que se dieron los primeros pasos hacia la carrera espacial. El principal artífice de esta visión fue Wernher Von Braun [2], amante de los relatos de aquel visionario escritor francés Julio Verne. Gracias a él, los nazis concibieron el misil balístico V-2, Vergeltungswaffe-2 (Arma de la venganza-2), destinado a bombardear y realizar el mayor daño posible al enemigo. El V-2 alemán fue el primer vehículo en abandonar la atmósfera terrestre el 20 de junio de 1944, superando la línea de Karman a más de 100 km de altura. A partir de este hito, unos años más tarde, comenzó la carrera espacial motivada por razones geopolíticas que propició el auge de las tecnologías espaciales. 1.1.1 Guerra Fría La Guerra Fría fue un enfrentamiento entre el bloque occidental-capitalista, liderado por Estados Unidos (EEUU) y el bloque oriental-comunista liderado por la Unión Soviética (URSS). Esta tensión surgió después de la Segunda Guerra Mundial por motivos económicos, políticos, sociales e ideológicos lo que llevó a ambas potencias a participar e influir en conflictos alrededor del mundo. Dicha intrusión por parte de tanto EEUU como la URSS desembocó en guerras como en Vietnam o Afganistán. La tensión constante entre ambos bloques llevó consigo una inversión armamentística y militar sin precedentes, incluida la exploración espacial. En esta época se produce el nacimiento de la carrera espacial [3], la cual merece que se analice cada potencia por separado. •Estados Unidos Estados Unidos fue pionero en poner en órbita los primeros satélites de comunicaciones. Programas como el Project Score, que en 1958 fue el primer satélite de comunicaciones puesto en órbita, el Syncom 2, primer satélite en órbita geoestacionaria y el Syncom 3, que se usó para retransmitir en directo en EEUU las Olimpiadas de 1964 de Tokyo. Tras las sucesivas conquistas soviéticas en la carrera espacial, Estados Unidos lanzó numerosas misiones para explorar la Luna, Pioneer, Ranger, Lunar Orbiter oApollo. El hito más conocido de la historia de la exploración espacial fue la llegada del hombre a la Luna el 16 de julio de 1969 a bordo del Apollo 11. 1 2Capítulo 1. Introducción •Unión Soviética La URSS logró el 4 de octubre de 1957 lanzar al espacio el primer satélite que orbitó la Tierra, el Sputnik 1. El satélite tenía el tamaño de un balón de baloncesto, pesaba 83 kilogramos y duró 20 días en órbita. Unos meses más tarde, el 3 de noviembre de 1957, los rusos consiguieron poner en órbita al primer ser vivo a bordo del Sputnik 2. Lamentablemente, la perra Laika no sobrevivó debido a un fallo en el control de temperatura del vehículo. Estos dos hitos al conmocionaron de la Guerra Fría conmocionó a la opinión pública estadounidense. En la carrera por llegar a la Luna, se desarrolló el programa de exploración Luna con el cual consiguieron acercarse a la inmediaciones selenocéntricas e incluso impactar en la Luna. Más adelante en 1961, la Unión Soviética fue de nuevo partícipe de la puesta en órbita del primer ser humano, Yuri Gagarin, a bordo de la Vostok 1. Años más tarde, en 1963 la primera mujer llegó al espacio. Valentina Tereshkova formó parte de una misión que duró 3 días, completando 48 vueltas alrededor de la Tierra. Por último, tras el gran hito norteamericano de alcanzar la Luna en 1969, los rusos lideraron numerosas misiones interplanetarias entre las que destaca el programa Venera que logró alcanzar Venus en 1970. El fin de la Guerra Fría tuvo lugar en 1975, simbolizado por el encuentro Apollo-Soyuz, una misión conjunta entre las agencias espaciales de ambos potencias. Figura 1.1 Cronograma carrera espacial. 1.1.2 Décadas de los 80-00 Una vez acabada la carrera espacial se llevaron a cabo misiones científicas y fueron puestas en órbita numerosas constelaciones de satélites para establecer sistemas de comunicaciones y navegación. En estos años organizaciones gubernamentales como la NASA, la ESA o la agencia espacial rusa ROSCOSMOS tenían todo el protagonismo en el sector espacial. Las misiones científicas tenían como objetivo el estudio y escudriño del sistema solar y sus astros. Ejemplos [4] como Galileo, alcanzando la órbita joviana en 1995; SOHO, principal fuente de información del estado del Sol que se encuentra en una órbita Halo en torno al punto de Lagrange L1; Cassini, que llegó a Saturno en 2004; o la misión New Horizons con destino Plutón siendo el vehículo más rápido jamás lanzado desde la Tierra. También, se llevó a cabo la construcción de la Estación Espacial Internacional en 1998 en una órbita de 350 kilómetros de altura. En ella se han llevado a cabo estudios en muchos ámbitos de la ciencia en un entorno de ingravidez. Más allá de las misiones científicas, se constituyeron las primeras constelaciones de satélites. Entre ellas destaca la constelación GPS (Global Positioning System), constituida por 31 satélites fue un proyecto del ejército de EE.UU. que permite determinar la posición y velocidad en cualquier punto del globo terráqueo. Por otro lado, Rusia lanzó su propia constelación, GLONASS para geolocalización. En el campo de las comunicaciones, constelaciones como Iridium, Inmarsar oGlobalstar, con cobertura global, supusieron una mayor seguridad y fiabilidad en sectores como el marítimo y el aéreo,. Con el paso de los años el número de satélites en órbita baja (LEO) ha aumentado enormemente y consigo la basura espacial obligando a los satélites a realizar maniobras para evitar colisiones. 1.2 Objetivos 3 1.1.3 Últimos años En las últimas décadas numerosas agencias espaciales se han sumado al sector espacial. Entre ellas la agencia espacial china (CNSA) o la agencia espacial india (ISRO) con constelaciones propias de satélites de navegación, BeiDou [5] y IRNSS [6]. Por otro lado, la ESA está en proceso de constituir la constelación Galileo en vista a sustituir al GPS. Sin embargo, la revolución que ha experimentado el sector espacial ha sido causada por la irrupción del sector privado, con empresas como SpaceX, Blue Origin oVirgin. Los programas espaciales siempre han sido financiados por presupuestos multimillonarios únicamente alcanzables por las agencias gubernamentales. Sin embargo, la reducción del coste en el lanzamiento mediante el uso de vehículos reusables y el aumento de la eficiencia de los mismos ha permitido la explotación privada del espacio. Además, se están realizando satélites cada vez más pequeños denominados nanosatélites cuyo objetivo es reducir la carga de pago enviada al espacio y por lo tanto abaratar el coste del lanzamiento. Figura 1.2 Coste de lanzamiento. Fuente (editada):[7]. En capítulos posteriores se describirá el estado actual de esta tecnología así como el potencial que tiene. La accesibilidad al espacio está aumentando y, con ello, la creación de nuevas empresas e incluso el turismo espacial. 1.2 Objetivos El presente trabajo tiene como objetivo principal modelar y diseñar el sistema de control para la reentrada de un vehículo de lanzamiento, con el fin de recuperar íntegramente el vehículo y así reutilizarlo en misiones posteriores. La línea de trabajo estará enfocada en simplificar el problema de tal manera que se pueda abordar académicamente sin abandonar la cercanía a la realidad que se pretende con este estudio. Se tomarán hipótesis justificadamente para diseñar los modelos con un nivel de profundidad suficiente para el correcto transcurso del problema. Por otro lado, aún siendo un sistema de 6 grados de libertad no lineales el problema de control se va a acometer con controladores sencillos y sin optimizar como una primera aproximación a un problema mucho más complejo. Parte de la finalidad de este proyecto es también actuar como punto de partida para trabajos futuros a raíz de los resultados o modelos propuestos. 4Capítulo 1. Introducción 1.3 Estructura Previamente al planteamiento del problema, se define el estado del arte de los sistemas de lanzamiento, tanto de agencias gubernamentales como de empresas privadas, llegando al por qué de la elaboración de este trabajo y su utilidad futura. Una vez especificado el motivo del estudio, se analizan las características del vehículo y se elaboran los modelos aerodinámico, propulsivo y de atmósfera con el suficiente grado de profundidad. Acto seguido, para finalizar la construcción del problema se exponen las ecuaciones que rigen el movimiento del vehículo en su descenso hacia la zona de aterrizaje. Un total de 6 grados de libertad gobiernan el vehículo de los cuales algunos son controlados y otros se restringen por simplificación. Seguidamente, se comprueba la validez del modelo mediante simulaciones donde se analiza el comportamiento frente a las deflexiones de dispositivos de control o la ignición de la planta propulsora. En estas simulaciones es de especial interés cómo se comporta el sistema trabajando en torno a una indeterminación. Una vez validado el modelo, es hora de diseñar los controladores para su gobierno. En primer lugar, se diseña el control por empuje con un control por relé y se prueba su viabilidad. En segundo lugar el control de actitud, más exigente que el anterior, se diseña mediante controladores proporcionales derivativos y de la misma manera se estudia su comportamiento analizando sus deficiencias. Por último el modelo con sus sistemas de control incorporados es simulado y analizado para distintos puntos de aterrizaje, lo que permitirá llegar a las conclusiones y trabajos futuros en vista a mejorar el proyecto. 2 Lanzadores El presente trabajo tiene como objetivo conseguir aterrizar un vehículo de lanzamiento. Esta tendencia comenzó recientemente con el desarrollo del Falcon 9 y se ha extendido rápidamente. Con este capítulo se pretende dar una visión de mercado de esta tecnología y su potencial. En él, se describirán los numerosos sistemas de lanzamiento utilizados para un amplio rango de misiones, desde misiones interplanetarias al turismo espacial. Además, se analizará el potencial que estos pueden tener en un futuro y la ventaja de reutilizarlos. 2.1 Organizaciones gubernamentales Desde los inicios de la era espacial, las agencias espaciales de algunos países han sido partícipes de los hitos espaciales más relevantes de la humanidad. Actualmente, más países se han unido al progreso espacial como China (CNSA) o la India (ISRO). En esta sección se van a analizar los sistemas de lanzamiento de las organizaciones gubernamentales más importantes. •NASA La NASA tiene previsto con el Space Launch System (SLS) [8] sustituir al transbordador STS. Con un coste estimado de 10000 millones de dólares, una altura de 100 m y una capacidad de carga de hasta 130 toneladas para alcanzar una órbita LEO, será el vehículo con mayor carga de pago que existe. Además, será posible reducir a la mitad o más la duración de misiones interplanetarias hacia las inmediaciones del sistema solar. •ESA La ESA está desarrollando el lanzador Vega-C [9], un cohete de 35 m de altura capaz de transportar una masa de 2200 kg a una órbita polar de 700 km de altura. Vega-C será el primer cohete que integre el vehículo reusable Space Rider[10] que una vez acabada la misión reentrará en la atmósfera terrestre de manera controlada. Este sistema de lanzamiento estará destinado a poner en órbita desde Cubesats a grandes satélites. •ROSCOSMOS Soyuz-2 [11] es un nuevo vehículo de lanzamiento que permitirá en el futuro reemplazar los vehículos de lanzamiento Soyuz-U, Soyuz-FG y Molniya-M. Con una altura de 46,3 m y una masa de 312 000 kg , permitirá lanzar naves espaciales en varios tipos de órbitas: baja, media, altamente elíptica, heliosincrónica, geotransicional y geoestacionaria. 2.2 SpaceX SpaceX, fundada por Elon Musk, es sin duda la empresa con proyectos más ambiciosos del sector. SpaceX ha ganado la atención mundial por una serie de hitos históricos. Es la única empresa privada capaz de recuperar una nave espacial desde la órbita terrestre baja, y en 2012 el vehículo Dragon se convirtió en la primera nave espacial privada en entregar carga en la Estación Espacial Internacional. También en 2020, se convirtieron en la primera empresa privada en llevar humanos al espacio. 5 6Capítulo 2. Lanzadores Tal es el éxito, que tiene acuerdos con la NASA y ha visitado 29 veces, hasta la fecha, la Estación Espacial Internacional con la Dragon Spacecraft. También tiene previsto participar en el programa Artemis de la NASA para volver a la Luna con el vehículo Starship. 2.2.1 Lanzadores Falcon El Falcon 9 es el lanzador escogido para su modelado en este trabajo y sus características serán definidas en detalle en el siguiente capítulo. En cuanto al Falcon Heavy [12], es un vehículo formado por 3 lanzadores Falcon 9 cuyos 27 motores Merlin generan más de 22000 kN de fuerza al despegue, equivalente a dieciocho aviones Boeing-747. Es capaz de transportar una carga de pago de 16800 kg a Marte. Tiene una altura de 70 m y una masa de 1,42 t. 2.2.2 Dragon La nave espacial Dragon [13] es capaz de transportar hasta 7 pasajeros hacia la órbita terrestre y más allá. Es la única nave espacial actualmente en vuelo que es capaz de devolver cantidades significativas de carga a la Tierra, y es la primera nave espacial privada que lleva humanos a la estación espacial. 2.2.3 Starship El vehículo espacial Starship [14] y el cohete Super Heavy de SpaceX (denominados colectivamente como Starship) representan un sistema de transporte totalmente reutilizable diseñado para transportar tanto a la tripulación como a la carga a la órbita terrestre, la Luna, Marte y más allá. Starship será el vehículo de lanzamiento más poderoso del mundo jamás desarrollado, con la capacidad de transportar más de 100 toneladas métricas a la órbita terrestre y con una altura de 100 m. Starship está diseñado para llevar satélites más lejos y a un costo marginal por lanzamiento más bajo que los vehículos Falcon actuales. Su capacidad de carga es tal que crea posibilidades para nuevas misiones, incluidos telescopios espaciales más grandes que el James Webb. Además, se tiene previsto que Starship participe en misiones hacia la Luna y Marte. De hecho, en 2023, el empresario japonés Yusaku Maezawa y la tripulación de dearMoon [15] se convertirán en los primeros pasajeros civiles en una misión lunar del Starship, con un sobrevuelo de la Luna durante su viaje de una semana. Este vuelo es un paso importante para permitir el acceso de las personas que sueñan con viajar al espacio. 2.3 Blue Origin Blue Origin fue fundada por Jeff Bezos con la visión de permitir un futuro donde millones de personas estén viviendo y trabajando en el espacio con el objetivo de beneficiar a la Tierra. Al igual que Space X, esta empresa defiende el reúso de sus vehículos de lanzamiento. Además, no descarta explotar los recursos que el Sistema Solar ofrece. 2.3.1 New Shephard Nombrado en honor al astronauta de Mercury Alan Shepard, el primer estadounidense en ir al espacio, New Shepard [16] es un sistemade lanzamiento suborbital reutilizable diseñado para llevar a los astronautas y cargas útiles de investigación más allá de la línea de Kármán, el límite del espacio reconocido internacionalmente. Con New Shephard, Blue Origin tiene como objetivo hacer del espacio un lugar más accesible para las personas y explotar el turismo espacial. 2.3.2 New Glenn Nombrado en honor al astronauta pionero John Glenn, New Glenn [17] es un vehículo de lanzamiento de carga pesada de configuración única capaz de transportar personas y cargas útiles de forma rutinaria a la órbita terrestre y más allá. Tiene una capacidad de transportar 13 toneladas métricas a una órbita estacionaria. Después de la separación de la primera etapa, esta vuela de regreso a la Tierra y aterriza en un barco en movimiento incluso con estados de mar gruesa. 2.4 Virgin 7 2.3.3 Blue Moon Blue Moon [18] es un módulo de aterrizaje flexible que entrega una amplia variedad de cargas útiles pequeñas, medianas y grandes a la superficie lunar. Ya sea carga o tripulación, su capacidad para proporcionar aterrizajes suaves y precisos permitirá una presencia humana sostenida en la Luna. Para llevar a cabo esta misión será necesario el diseño del sistema de control para un correcto y seguro aterrizaje del módulo de aterrizaje. Es por eso que el estudio de la entrada de vehículos y su aterrizaje en otros planetas como Marte o la Luna tiene un gran potencial en futuras misiones. Figura 2.1 Lanzadores actuales. Fuente:[17]. 2.4 Virgin Virgin Group es un conglomerado multinacional fundado por Richard Branson. Este conglomerado tiene dos filiares llamadas Virgin Orbit y Virgin Galactic, con una visión prometedora hacia el sector aeroespacial. 2.4.1 Virgin Orbit Más capaces que nunca, los satélites pequeños necesitan proveedores de lanzamiento que sean tan ágiles como las propias naves espaciales. Para satisfacer las necesidades de los fabricantes y operadores de satélites de hoy en día, Virgin Orbit ha desarrollado LauncherOne, el servicio de lanzamiento más flexible del mundo para satélites comerciales y construidos por el gobierno. Desde pequeñas empresas hasta grupos universitarios, Virgin Orbit está eliminando barreras al hacer que el espacio sea más accesible. Estos servicios "low-cost" están en resonancia con la tendencia actual de reducir el tamaño de satélites a cubos de 10 cm de lado permitiendo utilizar lanzadores más pequeños. Virgin utiliza un Boeing 747 modificado para transportar el sistema de lanzamiento en sus alas para así una vez en el aire, realizar el lanzamiento. 2.4.2 Virgin Galactic La filial Virgin Galactic está destinada a proporcionar vuelos suborbitales tripulados tanto para misiones científicas como para turismo espacial. El sistema de lanzamiento utilizado difiere de los actualmente vistos ya que se ha diseñado una aeronave para realizar el vuelo suborbital. Esta aeronave esta constituida a su vez por dos etapas, la primera etapa se asemeja a un avión convencional y realiza el vuelo en la atmósfera baja hasta una cierta altitud. Una vez alcanzada dicha altitud se desprende la segunda etapa, una aeronave supersónica que proyecta a los pasajeros a una altitud de 82 km para realizar un vuelo suborbital. 2.5 "Cubesats" Tendemos a pensar en los satélites como enormes naves espaciales que se elevan sobre los ingenieros que los construyen. Considerando Envisat, SOHO yGOCE: todos los satélites de la ESA son similares en tamaño 14 Capítulo 3. Modelización •Hipótesis 4: CD,bconstante para Mach >3 •Hipótesis 5: CD,fconstante para Mach >2 Estas hipótesis se han tomado para simplificar el problema aerodinámico en cierto grado sin obviar el régimen transónico y en base a resultados experimentales. Una vez tomadas dichas hipótesis queda modelar la transición del régimen subsónico al supersónico donde se produce un aumento puntual de la resistencia. 1. Base: De [29] se obtienen los siguientes resultados en función del número de Mach. Figura 3.6 Coeficiente de resistencia de la base en función de Mach. 2. Grid fins: Si se observa la figura 3.2 existe una gran variedad de grid fins con mallado distinto. Una mayor densidad de mallado proporciona un coeficiente de momentos mayor (mayor potencia de control) y un aumento de resistencia. En este problema interesa que ambos valores sean altos. No obstante, un coeficiente de momentos mayor, requiere que el actuador trabaje bajo momentos más exigentes. También, una mayor resistencia pone en juego la integridad estructural del sistema de control. El mallado escogido es el del propio vehículo Falcon 9 [21]. Para definir la función del coeficiente de resistencia del grid fin respecto al Mach se parte de [30], con la geometría estudiada representada en la figura 3.7a. Sin embargo, véase la diferencia entre el grid fin del estudio y el que se utiliza en el vehículo en la figura 3.7. La razón por la que se ha utilizado estos resultados es que al coeficiente de resistencia es adimensional con la superficie, por lo que la diferencias de superficie no afectan. Lo único que afecta es la densidad del mallado, el cual guarda una gran similitud con la superficie de control de referencia. (a) Geometría del grid fin del estudio. Fuente [30]. (b) Geometría del grid fin del vehículo. Figura 3.7 3.4 Modelo aerodinámico 15 La evolución del coeficiente de resistencia del grid fin sigue la misma tendencia, como cabe esperar, que la base. La gráfica representada en la figura 3.8 es relativa a una única superficie de control. Figura 3.8 Coeficiente de resistencia del grid fin en función de Mach. 3.4.1 Derivadas de estabilidad Una vez definido el coeficiente de resistencia del vehículo, el parámetro más importe del problema aerodinámico, es necesario definir los coeficientes aerodinámicos relativos a fuerzas transversales y momentos del vehículo. Estos coeficientes son denominados derivadas de estabilidad y por lo general se utilizan potentes softwares de CFD. No obstante, se ajustarán los valores empíricamente, en función a la evolución del sistema en la simulación o en base a ciertos estudios. Obsérvese que las derivadas de estabilidad dependen fuertemente de la geometría e inercia del vehículo. Figura 3.9 Momentos en ejes cuerpo. Sean Y una fuerza transversal y L, M y N los momentos respectos al eje x, y, z respectivamente (figura 3.9), se obtienen los siguientes coeficientes adimensionalizados, CY=Y 1 2ρV2SCl=L 1 2ρV2Sl Cn=N 1 2ρV2Sl Cm=M 1 2ρV2Sl (3.6) Si expresamos las ecuaciones de fuerzas y momentos, (CY)A=CY0+CYββ+CYδf δf(3.7) (Cl)A=Cl0+Clββ+Clδf δf(3.8) 16 Capítulo 3. Modelización (Cn)A=Cn0+Cnββ+Cnδf δf(3.9) (Cm)A=Cm0+Cmαα+Cmδf δf(3.10) El estudio de todas las derivadas de estabilidad se sale fuera de los límites de este trabajo, a continuación se van a definir las derivadas más importantes, CY0=0CYβ=1.5CYδf=0 Cl0=0Clβ=0Clδf=0 Cn0=0Cnβ=0.0012 Cnδf=0.0951 Cm0=0Cmα=−0.0012 Cmδf=−0.0951 En primer lugar, al ser un vehículo simétrico: CY0=Cl0=Cn0=0 . CYβ se obtiene de [31] donde dada la geometría del vehículo se ha basado este coeficiente a la resistencia que ofrece un cilindro. Por otro lado en lo respecto al momento L no se ha contemplando ningún acoplamiento respecto a las superficies de control o ángulo de resbalamiento para este trabajo. En cuanto al índice de estabilidad estática direccional y longitudinal, Cnβ y Cmα , se han tomado los valores empíricamente en función de la evolución del sistema durante la simulación. Para que el vehículo sea estable, en todo momento, Cnβ>0yCmα<0. Por último, las derivadas de control Cnδf y Cmδf se obtienen de la relación Sfin/S . No obstante, a la hora de la definición de las ecuaciones se verá con más detalle las derivadas de control y su sentido. 3.5 Modelo propulsivo El vehículo está compuesto por 9 motores Merlin 1D cuya distribución en planta se puede observar en la figura (3.10) y usan como combustible oxígeno líquido y queroseno RP-1 [21]. Figura 3.10 Distribución de la planta motora. La ecuación que define el empuje proporcionado por una tobera es E=˙mVs+As(ps−pa)(3.11) 1 Estos valores no se utilizan en las ecuaciones del problema pero son cualitativos para analizar la estabilidad el vehículo y su potencia de control 3.5 Modelo propulsivo 17 Dada esta ecuación, el empuje del motor cohete depende en mayor medida de las condiciones de los gases en la cámara de combustión siendo indiferente la velocidad de vuelo. Si bien, la presión atmosférica afecta negativamente al empuje. A continuación, se han tomado una serie de hipótesis con el objetivo de simplificar el modelo. •Hipótesis 6: Gasto másico constante El gasto másico se define como ˙m=ρ(s)V(s)A(s) donde s es el parámetro que define la sección de la tobera. •Hipótesis 7: Evolución isentrópica dentro de la tobera Dada una evolución isentrópica se cumple: T0 T=1+γ−1 2M2 s(3.12) p0 p=1+γ−1 2M2 sγ γ−1 (3.13) ρ0 ρ=1+γ−1 2M2 s1 γ−1 (3.14) Como consecuencia de la isentropía, las magnitudes de remanso son constantes en toda la tobera. Reformulando la ecuación del gasto másico y utilizando las ecuaciones (3.12) y (3.14): ˙m=ρ0a0M(s)A(s)1+γ−1 2M(s)2−γ+1 2(γ−1)(3.15) Si tomamos como referencia la sección de salida, con el gasto másico, el coeficiente de dilatación adiabático y las magnitudes de remanso constantes, en la ecuación (3.15), da lugar a un Mach de salida constante. A su vez, según las ecuaciones (3.12), (3.13) y (3.14), la temperatura, presión y densidad a la salida también serán constantes. Dadas las hipótesis 8 y 9, si las condiciones en la cámara de combustión son constantes, entonces el Mach, temperatura, presión y densidad de salida son constantes. •Hipótesis 8: Tobera adaptada a nivel del mar (ps=pa|SL) Tomadas estas hipótesis es posible definir los parámetros que caracterizan al motor cohete a partir de los datos [21]. Además, de la ecuación (3.11) se deduce que el único parámetro que varía el empuje es la presión atmosférica que a su vez dependerá de la altitud. En la figura (3.11) el empuje varía desde 914 kN a nivel del mar a 982 kN a medida que la altitud es mayor. Tabla 3.3 Parámetros del motor cohete. Parámetro Símbolo Valor Unidad Empuje a nivel del mar E|SL 914 kN Empuje en vacío E|V982 kN Impulso específico en vacío Isp|V310 s Área de salida As0.6799 m2 Gasto másico ˙m300.5857 kg/s Velocidad de salida Vs3041.0998 m/s Presión de salida ps101325 Pa 18 Capítulo 3. Modelización Figura 3.11 Empuje-Altitud. 4 Leyes del movimiento Es hora de plantear las leyes del movimiento que caracterizan al sistema. Sin embargo, antes de comenzar con su planteamiento, es necesario establecer los sistemas de referencia que se utilizarán en el problema: el sistema inercial topocéntrico, el sistema de ejes horizonte local, el sistema de ejes viento y el sistema de ejes cuerpo. Todo ello, para determinar la actitud del vehículo y su movimiento a lo largo de su trayectoria. Una vez definidos los sistemas de referencia, se plantearán las ecuaciones que rigen el movimiento del vehículo: ecuaciones dinámicas, ecuaciones cinemáticas de la posición y ecuaciones cinemáticas de la actitud. Un total de 6 ecuaciones, con 6 incógnitas: las velocidades en cada eje y los ángulos de Euler. 4.1 Sistema de referencia Para abordar las ecuaciones que afectan al vehículo es necesario establecer los distintos sistemas de referencia que van a participar en el problema. En primer lugar se definen los sistema de referencia inerciales donde a partir de los cuales se definirá el sistema de referencia relativo al viento y al vehículo. 4.1.1 Sistema Inercial Topocéntrico Para la definición de este sistema de referencia T(O, x, y, z) se toma la hipótesis de Tierra plana. Con ello, la superficie terrestre es plana y cualquier sistema de referencia fijado en ella (topocéntrico) es inercial, los ejes se definen como: •O (origen): cualquier punto de la superficie. •Eje x: Dirección Norte. •Eje y: Dirección Este. •Eje z: Completa un triedro a derechas (dirigido hacia abajo). 4.1.2 Sistema de ejes horizonte local El Sistema de ejes horizonte local H( Oh , xh , yh , zh ) surge de la translación del Sistema Inercial Topocéntrico al centro de masas del vehículo. Los ejes de horizonte local también forman un sistema inercial ya que su orientación coincide en todo momento con los ejes T(O, x, y, z). •Oh(origen): Centro de masas del vehículo. •Eje xh: Dirección Norte. •Eje yh: Dirección Este. •Eje zh: Completa un triedro a derechas (dirigido hacia abajo). 4.1.3 Sistema de ejes viento El sistema de ejes viento W(Ow,xw,yw,zw) es solidario con el movimiento del vehículo respecto al aire. •Ow(origen): Centro de masas del vehículo. 19 20 Capítulo 4. Leyes del movimiento •Eje xw: Tiene la misma dirección y sentido que el vector velocidad aerodinámica V. • Eje zw : Contenido en el plano de simetría del vehículo y dirigido hacia abajo para un ángulo de asiento (γ) nulo. •Eje yw: Completa un triedro a derechas. La orientación de los ejes viento con respecto a los ejes de horizonte local se puede definir a través de 3 ángulos: 1. Ángulo de asiento de velocidad ( γ ): es el ángulo formado por el vector velocidad aerodinámica  V y el plano del horizonte local. Se considera positivo cuando el vector apunta hacia abajo, a diferencia del criterio de signos utilizado en la literatura [32], [33], debido a que, por conveniencia, es más intuitivo trabajar en este problema con un ángulo de asiento apuntando hacia abajo (véase figura 4.1). 2. Ángulo de guiñada de velocidad ( χ ): es el ángulo formado por la proyección del vector velocidad aerodinámica  V sobre el plano horizonte local con la dirección norte. Se considera positivo hacia el este (sentido horario). 3. Ángulo de balance o alabeo de velocidad (µ): es el ángulo que forma el eje ywcon el plano horizonte local. Para el criterio de signo ver figura (4.2) Figura 4.1 Sistema de ejes horizonte local y ejes viento I. Figura 4.2 Sistema de ejes horizonte local y ejes viento II. Para la formulación de las ecuaciones será necesaria una matriz de giro que permita transformar magnitudes de un sistema a otro. La trasformación H−→ W , se obtiene a partir de los ejes horizonte local mediante una rotación de ángulo χ alrededor del eje zh , seguida de una rotación de ángulo γ alrededor del eje y′ , seguida de otra rotación de ángulo µalrededor del eje x′′. 4.1 Sistema de referencia 21 [T]W H = [T]WX′′ ·[T]X′′X′·[T]X′H=  1 0 0 0cosµsinµ 0−sinµcosµ ·  cosγ0sinγ 0 1 0 −sinγ0cosγ ·  cosχsinχ0 −sinχcosχ0 0 0 1 = =  cosγcosχcosγsinχsinγ −cosµsinχ−sinµsinγcosχcosµcosχ−sinµsinγsinχsinµcosγ sinµsinχ−cosµsinγcosχ−sinµcosχ−cosµsinγsinχcosµcosγ (4.1) 4.1.4 Sistema de ejes cuerpo El sistema de ejes cuerpo B(Ob,xb,yb,zb) se mueve solidariamente con el vehículo, se define como: •Ob(origen): Centro de masas del vehículo. •Eje xb: En la misma dirección que el eje de simetría del vehículo. • Eje zb : Contenido en el plano de simetría del vehículo y dirigido hacia abajo para un ángulo de asiento (γ) nulo. •Eje yb: Completa un triedro a derechas. Figura 4.3 Vista de planta del vehículo (a escala). La orientación de los ejes cuerpo con los ejes viento se puede definir mediante dos ángulos. 1. Ángulo de resbalamiento ( β ): es el ángulo formado por el vector velocidad aerodinámica  V con el plano xbzb, positivo cuando el aire entra por la derecha. 2. Ángulo de ataque ( α ): es el ángulo formado por la proyección del vector velocidad aerodinámica  V sobre el plano xbzby el eje xb La trasformación W−→ B , se obtiene a partir de los ejes viento mediante una rotación de ángulo β alrededor del eje zw , seguida de una rotación de ángulo α alrededor del eje yb , quedando la matriz de transformación, 22 Capítulo 4. Leyes del movimiento [T]BW =  cosαsinβ−cosαsinβ−sinα sinβcosβ0 sinαcosβ−sinαsinβcosα (4.2) Figura 4.4 Vista de alzado del vehículo (a escala). Figura 4.5 Sistema de ejes viento y ejes cuerpo. Por otro lado la orientación de los ejes cuerpo con respecto a los ejes horizonte local está dada por tres ángulos. 1. Ángulo de asiento, o pitch angle ( θ ): ángulo que forma el eje xb con el plano horizonte local, positivo hacia abajo. De la misma manera que anteriormente este criterio de signos, diferente de la literatura se ha escogido para un mayor sentido físico. Se cumple la relación θ=α+γ 2. Ángulo de guiñada, o yaw angle ( ψ ): ángulo que forma la proyección de xb sobre el plano horizonte local con la dirección norte, positivo hacia el este. Se cumple la relación ψ=β+χ 3. Ángulo de balance, o roll angle ( φ ): es el ángulo formado entre yb y el plano horizonte local, positivo hacia la derecha. 4.2 Ecuaciones del movimiento 23 La trasformación H−→ B , se obtiene a partir de los ejes horizonte local mediante una rotación de ángulo ψ alrededor del eje zh , seguida de una rotación de ángulo θ alrededor del eje y′ , seguida de otra rotación de ángulo φalrededor del eje x′′. [T]BH =  cosθcosψcosθsinψsinθ −cosφsinψ−sinφsinθcosψcosφcosψ−sinφsinθsinψsinφcosθ sinφsinψ−cosφsinθcosψ−sinφcosψ−cosφsinθsinψcosφcosθ (4.3) 4.2 Ecuaciones del movimiento A la hora de obtener las ecuaciones escalares del movimiento es necesario proyectar las ecuaciones vectoriales en forma matricial en un sistema de ejes determinado. Más adelante, para la determinación de las ecuaciones dinámicas, resultará necesario proyectar vectores definidos por un producto vectorial. Así, dado un vector ω y un tensor asociado Ω,  ω=  p q r Ω=  0−r q r0−p −q p 0 (4.4) que permite escribir  A= ω×a=Ω·a Además, se han tomado ciertas hipótesis para formular las leyes del movimiento que rigen sobre el vehículo. •Hipótesis 9: Tierra plana En primer lugar, se ha mencionado de manera discreta esta hipótesis a la hora de definir el sistema topocéntrico inercial. Además, el problema comienza a menos de 60 km de altura y se estudia exclusivamente el vuelo atmosférico, desatendiendo la trayectoria orbital que ha seguido antes de llegar al punto de inicio de este problema. De hecho, esta etapa de descenso es en mayor parte vertical, por tanto, la hipótesis de Tierra plana es viable. •Hipótesis 10: Gravedad constante La gravedad solamente disminuye un 1.55 % a una altura de 60 km por lo que también es una buena aproximación inicial. •Hipótesis 11: Posición del centro de gravedad constante Si bien se tiene en cuenta la disminución del peso debida a la combustión del combustible, a la hora de plantear las ecuaciones dinámicas no se tendrá en cuenta el cambio de posición del centro de gravedad. 4.2.1 Ecuaciones dinámicas En primer lugar se va llevar a cabo la formulación del equilibrio de fuerzas del sistema. Se parte de la formulación matricial, m d V dt !I = ( FA,T)B+m(g)I(4.5) Para obtener las ecuaciones dinámicas escalares se hará uso del siguiente resultado de la Mecánica del Sólido Rígido. Sea un sistema de referencia fijo F y otro móvil M, la relación entre las derivadas temporales viene dada por d A dt !F = d A dt !M + ωMF × A(4.6) Este resultado aplicado a la velocidad en el sistema de ejes cuerpo y el sistema topocéntrico se obtiene d V dt !I = d V dt !B + ωBI ×( V)B(4.7) y sustituyendo en (4.5) 30 Capítulo 4. Leyes del movimiento λ1>0−→ |λ1|=µ4 λ1<0−→ |λ1|=µ2 λ2>0−→ |λ2|=µ3 λ2<0−→ |λ2|=µ1 A su vez, aplicando las propiedades del coseno a la ecuación 4.47 se obtiene, FA,T x=−T−1 2ρV2SbaseCdb −ρV2Sf inCd f (cosλ1+cosλ2)u V(4.56) Por otro lado, el consumo de combustible debido al uso de la tobera se cuantifica según la ecuación: ˙m=−300.5857 (4.57) donde el valor del gasto másico fue definido en el modelo propulsivo. Para finalizar se van a resumir todas las hipótesis y las ecuaciones alcanzadas en este capítulo. Hipótesis •Hipótesis 1: Atmósfera en calma •Hipótesis 2: Atmósfera exponencial •Hipótesis 3: Cdb yCd f constentes para Mach∈[0,0.6] •Hipótesis 4: Cdb constante para Mach>3 •Hipótesis 5: Cd f constante para Mach>2 •Hipótesis 6: Gasto másico constante •Hipótesis 7: Evolución isentrópica dentro de la tobera •Hipótesis 8: Tobera adaptada a nivel del mar (ps=pa|SL •Hipótesis 9: Tierra plana •Hipótesis 10: Gravedad constante •Hipótesis 11: Posición del centro de gravedad constante Ecuaciones del sistema 1. Ecuaciones dinámicas   ˙u ˙v ˙w =1 m  FA,T+mg ·sinθ FA,T+mg ·sinφcosθ FA,T+mg ·cosφcosθ +  r·v−q·w p·w−r·u q·u−p·v    ˙p ˙q ˙r =  (L−qr(Izz −Iyy))/Ixx (M−pr(Ixx −Izz))/Iyy (N−pq(Iyy −Ixx))/Izz   4.2 Ecuaciones del movimiento 31 2. Ecuaciones cinemáticas de la posición dx dt = (cosθcosψ)u+(−cosφsinψ−sinφsinθcosψ)v+(sinφsinψ−cosφsinθcosψ)w dy dt = (cosθsinψ)u+(cosφcosψ−sinφsinθsinψ)v+(−sinφcosψ−cosφsinθsinψ)w dh dt =−[(sinθ)u+(sinφcosθ)v+(cosφcosθ)w] 3. Ecuaciones cinemáticas de la actitud     ˙q0 ˙q1 ˙q2 ˙q3     =1 2    −q1−q2−q3 q0−q3q2 q3q0−q1 −q2q1q0      p q r  4. Ecuaciones de fuerzas y momentos FA,T x=−T−1 2ρV2SbaseCdb −ρV2Sf inCd f (cosλ1+cosλ2)u V FA,T y=1 2ρV2SlatCYβ v V FA,T z=1 2ρV2SlatCYβ w V L=0 M=1 2ρV2SfinCd f sinλ1 u V+1 2ρV2SfinCmαα N=1 2ρV2SfinCd f sinλ2 u V+1 2ρV2SfinCmββ ˙m=−300.5857 5 Modelado en Simulink Una vez establecidas las ecuaciones que rigen el movimiento del vehículo, es hora de implementarlo en Simulink. En este capítulo, se va a describir el procedimiento seguido para implementarlo en la herramienta. Además se analizará el comportamiento del sistema para distintas configuraciones antes de abordar el problema de control. 5.1 Estructura Dadas las ecuaciones desarrolladas en el capítulo anterior, para su implementación en Simulink, se va a seguir la estructura definida en el diagrama de bloques de la figura 5.1. Dentro de los bloques, la estructura es similar, integrando las ecuaciones definidas para cada uno. A continuación, se describirá el contenido de cada bloque y su función. Figura 5.1 Diagrama de bloques. Por ahora el bloque Dispositivos de control no está realimentado por la actitud del vehículo porque todavía no se han definido el control de estas variables. Se ha incluido de esta manera para poder realizar pruebas con valores determinados de las variables de control para analizar cómo responde el sistema y su comportamiento. 5.1.1 Fuerzas y Torques Los parámetros que recibe este bloque son las componentes de la velocidad en ejes cuerpo u, v, w y las variables atmosféricas ρ , p, T, M. También, en este bloque actuarán los dispositivos de control del empuje y 33 34 Capítulo 5. Modelado en Simulink de la actitud, modificando las fuerzas y momentos del vehículo y por tanto las velocidades y ángulos del mismo. La función de este bloque será por tanto calcular las fuerzas y momentos del cohete según las ecuaciones 4.47, 4.48, 4.49, 4.50, 4.54, 4.55, y 4.57. Figura 5.2 Fuerzas y torques. 5.1.2 Ecuaciones del movimiento Los parámetros recibidos por este bloque son las velocidades lineales y angulares en ejes cuerpo u, v, w, p, q, r, la actitud del vehículo a través de los ángulos de Euler θ , ψ , φ y las fuerzas y momentos Fx , Fy , Fz , L, M, N, m. Siguiendo las ecuaciones dinámicas 4.10, 4.14 se obtienen ˙u,˙v,˙w,˙p,˙q,˙r. A continuación se integra con el bloque integrator que tiene Simulink y se obtienen los valores u, v, w, p, q, r. Figura 5.3 Ecuaciones del movimiento. 5.1 Estructura 35 5.1.3 Posición y actitud Una vez calculadas las componentes de la velocidad y velocidad angular en ejes cuerpo u, v, w, p, q, r queda determinar la posición en ejes horizonte local x, y, z y la actitud del vehículo θ , ψ , φ . No obstante, en este problema es necesario trabajar con cuaterniones por lo que habrá que realizar conversiones de ángulos de Euler (θ,ψ,φ) a cuaterniones y viceversa. Figura 5.4 Posición y actitud. Para la conversión de ángulos de Euler a cuaterniones se construye la matriz de giro siguiendo la estructura de 4.3. Una vez definida la matriz de giro se aplica la ecuación 4.38 de donde se obtienen los cuatro componentes del cuaternión q0,q1,q2,q3. Por otro lado, el método inverso no es tan trivial debido al carácter dual de las funciones trigonométricas. Siguiendo la fórmula de Euler-Rodrigues para el paso de cuaterniones a la matriz de giro 4.39, se obtiene la matriz que tiene como estructura la definida en 4.3. Además, es necesario tener en cuenta el caso de θ=90o donde se produce la pérdida de un grado de libertad, razón principal por la que se ha escogido trabajar con cuaterniones. Figura 5.5 Diagrama de flujo para conversión a ángulos de Euler. 36 Capítulo 5. Modelado en Simulink Si se denota la as filas y columnas de la matriz de giro 4.3 como, C=  A B C D E F G H I  (5.1) es inmediato deducir el valor de θ, θ=arcsin(C)(5.2) Una vez obtenido el valor de θ, siguiendo el esquema 5.5, pueden existir dos casos •θ=90o se pueden obtener ψyφde forma inmediata como, ψ1=arcsin(B) cosθ(5.3) φ2=arcsin(F) cosθ(5.4) Sin embargo para el arcoseno de un ángulo existen 2 soluciones por lo que habrá que calcular las soluciones restantes, ψ2=π−ψ1(5.5) φ2=π−φ1(5.6) Para conocer cual de las soluciones es la correcta, hace falta construir una matriz de giro para cada combinación de las soluciones posibles y compararla con la matriz de giro C calculada mediante la fórmula de Euler-Rodrigues. Por tanto se construirán C1(θ,ψ1,φ1) , C2(θ,ψ1,φ2) , C1(θ,ψ2,φ1) , C1(θ,ψ2,φ1) . Una vez construidas las matrices, aquella que cumpla C−Ci=0 será la solución y por consiguiente los ángulos de Euler (θ,ψ,φ). Sin embargo, la existencia de error numérico a lo largo del proceso hace que C−Ci=0 se cumpla con menos frecuencia. Por tanto, es más fiable comparar las normas de las diferencias y escoger la menor de ellas, es decir, min(C−Ci). •θ=90o Con esta configuración no se pueden calcular los ángulos ψ y φ según 5.3 y 5.4. La pérdida de un grado de libertad conlleva la existencia de una relación entre ψ y φ . Si desarrollamos el término D de la matriz de giro C para θ=90o, se obtiene, D=−cosφsinψ−sinφcosψ=−sin(φ+ψ)(5.7) Desarrollando igualmente los términos E, G, y H de la matriz de giro C se obtiene el mismo resultado. Tomando φ=0, ψ=arcsin(−D)(5.8) 5.2 Comportamiento del sistema Una vez modelado el sistema en Simulink, se va a analizar el comportamiento del vehículo para distintas configuraciones para verificar que efectivamente el modelo tiene sentido físico y está bien implementado. Además, también servirá para entender su evolución a la hora de la implementación del control automático que más tarde se va a llevar a cabo. Las condiciones iniciales se deducen de [35]. La posición inicial es ( x0=0m , y0=0m , z0=54000 m ), la actitud inicial ( θ0=90o , ψ0=0o , φ0=0o ) y la velocidad inicial (u0=1389 m/s,v0=0m/s,w0=0m/s). 5.2 Comportamiento del sistema 37 5.2.1 Caída libre vertical En primer lugar, se simulará el caso más simple, donde los grid fins no son accionados ni tampoco la planta propulsora. El objetivo de esta simulación es analizar cómo evoluciona la altura y la presión dinámica ( 1/2ρ·U2 ) a la que se ve sometida el vehículo. Además, es necesario comprobar que en ausencia de perturbaciones, el sistema partiendo de su posición de equilibrio no sale de ella y se descontrola a lo largo de la trayectoria de descenso hacia tierra. Será interesante observar cómo varía la velocidad sin ignición alguna en la planta propulsora para analizar el papel de la atmósfera en la cinemática del vehículo. (a) Trayectoria. (b) Velocidad. Figura 5.6 (a) Altura. (b) Presión dinámica. Figura 5.7 Evolución en caída libre. El parámetro más crítico se encuentra en la figura 5.7b donde la presión dinámica ( q=1/2ρ·U2 ) toma un valor máximo de 292000 kg/(m·s) a una altura de 12000 m . Este valor se alcanza en la tropopausa donde la presencia de la atmósfera se hace evidente. Para un análisis más completo del tramo de descenso es necesario realizar un estudio estructural del vehículo para conocer los esfuerzos, las temperaturas o las vibraciones límite a las que se puede ver sometido el vehículo. No obstante, este campo de estudio se sale fuera de los límites de este trabajo. Dada la ausencia del análisis estructural se va a tomar la evolución seguida en una operación real con este vehículo de [35] que se especificará en la siguiente subsección. Por otro lado, como cabe esperar, los ángulos que definen la actitud no sufren modificación alguna. Este resultado es una señal de la robustez del modelo pues con una ligera modificación de θ , existiría un salto en las variables ψyφcomo se verá en análisis posteriores. 38 Capítulo 5. Modelado en Simulink (a) θ. (b) ψ. (c) φ. Figura 5.8 Actitud. Por último, obsérvese en la figura 5.6b que aumenta ligeramente la velocidad del vehículo dada la escasez de la atmósfera inicial. Conforme sigue el descenso la velocidad disminuye considerablemente ligado a una mayor presencia de la atmósfera. En el tramo final, atraviesa el régimen transónico para impactar contra el suelo a una velocidad de 960 km/h. 5.2.2 Encendido de planta propulsora Uno de los objetivos principales durante el descenso es reducir la presión dinámica a la que se ve sometido el vehículo. La manera de reducirla es disminuir la velocidad por lo que será necesario activar la planta propulsora durante un intervalo de tiempo. Para ello se debe llevar a cabo un proceso de optimización para conseguir una trayectoria con la mínima presión dinámica. Sin embargo, en lugar de eso se ha estimado de [35] la duración del encendido inicial 19.3s . Durante el encendido de la planta propulsora sólo se activará la tobera 9 (obsérvese la figura 3.10). En primer lugar remarcar la similitud del modelo con la realidad. En la operación real del vehículo [35] tras el apagado del motor el cohete se encuentra a una altitud de 33500 m y a una velocidad de 2702 km/h . En el modelo, tras el apagado del motor el cohete se encuentra a una altitud de 33012 m y a una velocidad de 2700 km/h , una variación media del 0.76 % . Si se analiza más adelante, en el sistema real a una altitud de 8000 m se tiene una velocidad de 1428 km/h , mientras que para la simulación se tiene una velocidad de 1470 km/h, una variación del 2.94 %. La consecuencia más evidente del encendido de la planta propulsora durante casi 20 segundos en la fase inicial del descenso se encuentra en la presión dinámica. El pico máximo se reduce a 95260 kg/(m·s) dado a una altitud de 12390 m , un 67.4 % menos. Obsérvese la figura 5.9b para apreciar la reducción de magnitud. 5.2 Comportamiento del sistema 39 (a) Comparación de alturas. (b) Comparación de presiones dinámicas. (c) Comparación de velocidades. (d) Evolución de la masa. Figura 5.9 Comparación. Por otro lado, la duración del aterrizaje aumenta hasta superar los 80 segundos y la evolución de la velocidad no es tan abrupta, por lo que las aceleraciones que experimentará el cohete serán menores. También tener en cuenta que debido al encendido de la planta propulsora la masa disminuye y por tanto también los momentos de inercia. Esta disminución se traduce en una mayor facilidad para producir un momento en el vehículo mediante los grid fins. Por último es necesario un último encendido al final del descenso para aterrizar de forma segura pero ello se analizará en el siguiente capítulo. La deflexión de los grid fins se realizará una vez se produzca el primer apagado a 33000 mde altura. 5.2.3 Deflexión de superficie de control 1 Una vez realizada la maniobra de frenado del vehículo a 33000 m se lleva a acabo el guiado aerodinámico del cohete hacia el punto de aterrizaje. Para analizar el comportamiento de las superficies aerodinámicas o grid fins, se va a realizar una simulación partiendo de las mismas condiciones iniciales impuestas al principio. La prueba consistirá en desplegar 30o el grid fin 1, ver figura 5.10, durante 5 segundos una vez que la etapa de encendido inicial ha finalizado. El objetivo de la prueba es analizar cómo responde el modelo al salir de θ=90oy su evolución. En primer lugar, al deflectar la superficie de control 1 se produce un momento alrededor del eje zb , en sentido positivo (véase la figura 5.10). Una vez que el cohete abandona la configuración vertical y se inclina, genera una fuerza de resistencia lateral de manera que el vehículo se desplaza a lo largo de yb , en sentido contrario. obsérvese que con una deflexión de tan sólo 5 segundos, se ha conseguido un desplazamiento de poco más de 50 metros (figuras 5.11a y 5.11b). Este resultado da a entender que el rango de operación de esta misión no debería de sobrepasar los 60 metros. 46 Capítulo 6. Sistema de control del empuje del sistema de control será pues seguir dicha evolución de la velocidad. Si se observa la figura 6.2, aparecen dos "picos". Estos, son los puntos donde se encuentra la transición entre tobera encendida y apagada. En primer lugar, entre los 54000 m y 33000 m se produce el primer encendido descrito en el capítulo anterior. A continuación el vehículo evoluciona libremente hasta que a una altitud de aproximadamente 2100 m se produce el último encendido. Para ello, el controlador del relé recibe la diferencia entre la velocidad deseada u∗ y la velocidad del sistema u . Si la diferencia es menor que 0, es decir, la velocidad del vehículo es mayor que la deseada, la salida del relé adquiere el valor 1. Si por el contrario la diferencia es mayor que 5, es decir, la velocidad del vehículo es ligeramente menor que la deseada, la salida del relé será entonces de valor 0 (véase la figura 6.3). Acto seguido el controlador de la tobera recibe la salida del control por relé, siendo el valor 1 encendido y 0 apagado. Figura 6.2 Diagrama de bloques del control de empuje. Figura 6.3 Arquitectura del relé. La complejidad de este sistema de control no es elevada, no obstante, se ha priorizado que la ignición de la planta propulsora se produzca una sola vez antes del contacto con el suelo. Para llevarlo a cabo es necesario que la velocidad deseada u∗ y el margen establecido en el controlador tengan una precisión tal que se consiga la evolución deseada. De lo contrario, el relé estaría continuamente encendiéndose y apagándose, algo válido y quizás más preciso en la práctica, pero que se aleja de un comportamiento real del cohete. En el sistema real cuantos más encendidos y apagados existan la probabilidad de fallo aumenta y la precisión disminuye porque existe un período de transición entre la ignición y el punto de operación de la tobera, hecho que no se contempla en el control por relé. 6.3 Resultados A continuación se van a exponer los resultados obtenidos, partiendo de unas condiciones iniciales tales que la posición inicial es ( x0=0m , y0=0m , z0=54000 m ), la actitud inicial ( θ0=90o , ψ0=0o , φ0=0o ) y la 6.3 Resultados 47 velocidad inicial ( u0=1389 m/s , v0=0m/s , w0=0m/s ). Dichas condiciones no difieren de las definidas en el capítulo anterior. Siguiendo las especificaciones de los motores de [21], los cuales tienen capacidad de regular el empuje desde un 100 % hasta un 43 % del empuje máxima operacional de la tobera, se ha decido trabajar en el encendido final al 43 % de su capacidad máxima. La razón principal se debe a que con un empuje menor se consigue una desaceleración del vehículo menor. La ventaja que proporciona es la existencia de un mayor margen en el momento de ignición de la tobera en el caso de que se produjera un adelanto o retraso en la ignición. Si la desaceleración fuera elevada, un retraso de 0.5 segundos podría suponer que la velocidad de contacto con el suelo fuera 50 km/h más elevada de la esperada porque el sistema evoluciona más rápidamente y por tanto la pendiente de la velocidad es más acuciada. Otra razón por la cuál se ha disminuido la potencia en el tramo final es para seguir fielmente el comportamiento del sistema real de [35]. (a) Empuje. (b) Masa. Figura 6.4 Control del empuje. (a) Velocidad. (b) Velocidad en instantes finales. Figura 6.5 Evolución de la velocidad. En primer lugar, el empuje y la masa del vehículo representados en las figuras 6.4a y 6.4b reflejan lo que se ha mencionado anteriormente. En el primer encendido, se trabaja con el empuje máximo operacional, mientras que en el segundo encendido al 43 % de la potencia máxima. También se aprecia con el combustible consumido, donde en el segundo encendido se aprecia una variación de la masa frente al tiempo más suave que en la primera ignición. Por otro lado, al modelizar la planta propulsora, se tuvo en cuenta la presencia de atmósfera y la disminución del empuje disponible con el aumento de la presión atmosférica, no obstante, esta no afecta en gran medida al empuje que, a grandes escalas, parece constante. Además, la masa restante es cercana a la masa estructural del vehículo de 19000 kg por lo que se utiliza el combustible justo y necesario para la maniobra de descenso. En el diseño del sistema de control se ha priorizado la existencia de una sola etapa de ignición final. Si 48 Capítulo 6. Sistema de control del empuje bien se aprecia un encendido y apagado en el instante final, la evolución del empuje es satisfactoria. No obstante comparando con el caso real de [35], la ignición final se produce a 4100 m de altura, mientras que en la simulación la ignición se produce a 2100 m . Esto puede ser debido a que el modelo propulsivo escogido no se comporta fielmente como lo hace en el sistema real, ya que hasta antes de la ignición las velocidades de tanto el real como el simulado son muy parecidas. Este retraso en la ignición durante la simulación indica que la potencia de que se dispone es mayor que en el vehículo real. En cuanto a la evolución de la velocidad, obsérvese los dos "picos" que aparecen en la figura 6.5a que se relacionan con el encendido y apagado de la planta propulsora. También véase en la figura 6.5b cómo disminuye la velocidad de 900 km/h a 3km/h en 18 segundos. Además, es apreciable el apagado y la pequeña ignición final en la velocidad. Tras varias pruebas intentando evitar este rápido encendido y apagado no ha sido posible mejorar la evolución del sistema. Los instantes finales son críticos y este resultado es en el que se ha conseguido una velocidad de impacto mínima. Por último, una vez que se produce el encendido de la tobera, la presión dinámica disminuye rápidamente por la desaceleración a la que se ve sometida el vehículo y la curva de la altura frente al tiempo se suaviza. Por tanto, el control por relé implementado para llevar a cabo la tarea de aterrizar el vehículo de forma segura es satisfactorio, considerando una velocidad de contacto con el suelo de 10 km/h como segura. Es un buen resultado teniendo en cuenta la sencillez del controlador y la complejidad del sistema ya que un segundo de adelante o de retraso pueden suponer velocidades de impacto considerables o que el vehículo no llegue a tocar el suelo y vuelva a ascender. 7 Sistema de control de la actitud Para completar la maniobra de descenso queda diseñar el piloto automático para la actitud. A diferencia del control por relé, este sistema de control estará formado por dos controladores proporcionales derivativos PD en cascada. A lo largo del capítulo se explicará la estructura del sistema, los valores para su sintonización y su comportamiento. Asimismo, se diseñará una trayectoria adecuada para guiar al vehículo en su descenso. Por último se analizarán los resultados de las simulaciones comentando las limitaciones, deficiencias y los logros. 7.1 Objetivo El objetivo principal es conseguir el seguimiento de una trayectoria diseñada con el fin de llegar al punto de aterrizaje. Dada la complejidad del sistema se han centrado los esfuerzos en el estudio de un sistema de control para el problema bidimensional. Reduciendo el problema a dos dimensiones, se utilizan dos de los cuatro grid fins mediante una variable de control ( λ1 o λ2 ), modificando así el ángulo de asiento y consigo el ángulo de trayectoria. No obstante, el vehículo no puede alcanzar el punto de aterrizaje de cualquier manera. Además, hay que tener en cuenta que la ignición de la planta propulsora con un ángulo de asiento distinto de 90o puede causar un aumento indeseado de la velocidad horizontal y afectar notablemente a la precisión del sistema de aterrizaje. 7.2 Estructura Antes de analizar la arquitectura que se va a seguir para el diseño del controlador es necesario profundizar en el concepto de ángulo de asiento θ y ángulo de trayectoria γ . En primer lugar el ángulo de asiento es el ángulo que forma el eje de simetría del vehículo con el plano horizontal y el ángulo de trayectoria se define como el ángulo formado por el vector velocidad con el plano horizontal, véase la figura 7.1. La arquitectura seguida para el diseño del piloto automático se caracteriza por la implementación de dos controladores PD en cascada, el controlador para el ángulo de asiento ( PD1 ) y el controlador para el ángulo de trayectoria ( PD2 ). Para seguir la trayectoria definida, el sistema recibe el ángulo de trayectoria deseado el controlador PD2 calcula el ángulo de asiento adecuado y el controlador PD1 calcula la deflexión necesaria de la superficie de control, véase la figura 7.2. Para la construcción de la arquitectura propuesta, se va a diseñar, en primer lugar, el controlador para el ángulo de asiento y acto seguido el controlador para el ángulo de trayectoria que completa el piloto automático. 49 50 Capítulo 7. Sistema de control de la actitud Figura 7.1 Ángulo de trayectoria γy ángulo de asiento θ. Figura 7.2 Arquitectura del piloto automático. 7.2.1 Diseño del controlador para el ángulo de asiento El objetivo de este controlador es alcanzar el ángulo de asiento deseado con la menor sobreoscilación y el menor error posibles. Para el diseño del controlador se trabajará en el plano xz, es decir, en la orientación Norte-Sur aunque es igualmente aplicable para el plano yz, orientación Este-Oeste. En primer lugar, es necesario conocer cómo están definidos los ángulos de Euler. El ángulo de asiento varía entre 0o y 90o , mientras que el ángulo de guiñada ψ toma valores de 0o si el vehículo apunta hacia el Norte o 180o si apunta hacia el Sur. Sin embargo, el controlador solo trabaja con la variable θ y para los mismo valores de θ existen dos casos, que el vehículo esté apuntando hacia el Norte o hacia el Sur, véase la figura 7.3a. Para hacer posible el control de esta variable se va a redefinir el ángulo de asiento. Este tomará valores entre 0o y 180o con el criterio de signos tomado en la figura 7.3b. Al tomar este nueva definición del ángulo de asiento, si bien pierde algo de sentido físico, permite al controlador trabajar en el plano xz prescindiendo de la variable ψy aumentando el intervalo de trabajo de θ. Una vez realizados los cambios necesarios, es hora de sintonizar el controlador. En primer lugar, este sistema se caracteriza por evolucionar lentamente, además de tener una gran inercia, sobretodo en las capas altas atmosféricas (por encima de los 20000 m) debido a la escasez de aire. Por tanto el controlador no debe tener sobreoscilación y debe de tener un carácter predictivo que le permita anticiparse a los momentos de inercia del vehículo. Es por esta razón por la que se va a trabajar con un controlador proporcional derivativo (PD1) que se caracteriza por: •Incorporar predicción de la tendencia del error • Aumentar el valor de Td puede aumentar la estabilidad estática relativa (disminuir la sobreoscilación) y permitir aumentar la velocidad de respuesta del sistema •Valores grandes de Tdamplifican el ruido con efectos negativos sobre los actuadores • Aumentar el valor proporcional P aumenta la velocidad de la respuesta, aumenta la sobreoscilación y disminuye el error en permanente 7.2 Estructura 51 (a) Definición para ángulos de Euler. (b) Definición para entrada en controlador. Figura 7.3 Criterio de signos. Por lo tanto, a la hora de sintonizar el controlador para el ángulo de asiento se tendrá el compromiso de aumentar la velocidad de la respuesta con la menor sobreoscilación y ruido posibles. Los valores obtenidos tras la sintonización del controlador son, P=170 D=750 (7.1) Además a fin de no exceder las limitaciones mecánicas de los actuadores se ha impuesto una saturación al controlador de ±π/4. (a) Entrada del controlador. (b) Evolución de θ. Figura 7.4 El sistema se encuentra en equilibrio para un ángulo de asiento de 90o . Para una entrada en escalón como la de la figura 7.4a, se obtiene una respuesta rápida, si bien con algo de sobreoscilación y un pequeño error en régimen permanente. Es inmediato intuir que un aumento en valor derivativo del controlador disminuirá la sobreoscilación y por tanto se alcanzará antes el régimen permanente. Sin embargo, se ha preferido no aumentarlo ya que ello puede provocar la aparición de ruido en el controlador. Por otro lado obsérvese la respuesta del controlador en la figura 7.5a. En ella no aparece ruido, aún siendo el valor derivativo del controlador elevado. También señalar que aunque se esté controlando una única variable, λ1 , a la hora de trasladar los resultados al movimiento de las superficies del vehículo los grif fins 2 y 4 girarán según las gráficas 7.5b y 7.5c. Remarcar que se han tomado que todos los actuadores giran en el mismo sentido. 52 Capítulo 7. Sistema de control de la actitud (a) Salida del controlador. (b) Actuador 2. (c) Actuador 4. Figura 7.5 Respuesta del controlador para el ángulo de asiento. Aunque la salida del controlador sea un única variable λ1 , es posible traducir esta respuesta al sistema diseñado en el problema donde se controla la actitud en un plano con dos actuadores. Tal y como está definido el problema, la parte positiva de la gráfica corresponderá al actuador 2 y la parte negativa al actuador 4. A partir de este punto, se representarán las gráficas referentes a los actuadores 2 y 4 para un mayor sentido físico. Tener en cuenta que en todo momento la salida del controlador es solo una, λ1. 7.2.2 Diseño del controlador para el ángulo de trayectoria Una vez el ángulo de asiento θ está bajo control, queda diseñar el controlador para el ángulo de trayectoria que completa el piloto automático (figura 7.2). La evolución del ángulo de trayectoria en este sistema tiene unas características muy similares la evolución antes estudiada del ángulo de asiento. Es necesario incorporar una predicción de la tendencia del error, además de una respuesta rápida y precisa. Es por eso, que se va a utilizar un controlador proporcional derivativo PD2, cuyos valores tras la sintonización son, P=10 D=30 (7.2) Además se ha introducido una saturación de ±0.6o ya que dada la evolución lenta del sistema, una variación superior del ángulo de asiento aumentaría el tiempo de respuesta hasta situaciones en las que el tiempo de simulación sea inferior al tiempo de respuesta. Al igual que con el ángulo de asiento, el ángulo de trayectoria se encuentra en equilibrio para γ=90o , por lo que la función escalón de entrada comienza en ese valor hasta el salto a γ=90,75o . Puede parecer un ángulo pequeño pero es debido a la gran velocidad de caída del vehículo. Si se quisieran ángulos de trayectoria mayores, el ángulo de asiento θ se alejaría aún más de su posición nominal de 90o . Ello complicaría finalizar la maniobra con γ=θ=90o pues habría más error en régimen permanente y el tiempo de respuesta aumentaría. 7.2 Estructura 53 Obsérvese en la figura 7.6b la evolución de γ . Experimenta una pequeña sobreoscilación y un error en régimen permanente que se podrían mitigar con valores mayores tanto del valor proporcional como del derivativo, no obstante, a fin de evitar la aparición de ruido se han mantenido los valores de 7.2. (a) Entrada al piloto automático. (b) Evolución de γ. Figura 7.6 Respuesta del controlador para el ángulo de trayectoria. (a) Salida de PD2(θre f ). (b) Evolución de θ. Figura 7.7 Evolución de θ. (a) Actuador 2. (b) Actuador 4. Figura 7.8 Actuadores. 54 Capítulo 7. Sistema de control de la actitud En las figuras 7.7a y 7.7b es posible ver cómo el vehículo no responde muy fielmente a cambios bruscos en el ángulo de asiento. Sin embargo para cambios progresivos en θre f el vehículo puede seguir correctamente su evolución. Por otro lado, es apreciable la aparición de algo de ruido en la respuesta del controlador para el ángulo de asiento ( PD1 ) tanto en el actuador 2 como en el actuador 4. Es un ruido que no aparece a la hora de sintonizar el controlador con la entrada en escalón de θre f pero que, sin embargo, si que aparece una vez que el vehículo debe de seguir una evolución específica de θre f . Existen varias opciones para reducir el ruido, disminuyendo el valor del término derivativo o aumentando el valor del término integral. No obstante, la disminución del término derivativo empeora notablemente la respuesta del controlador al igual que introduciendo un valor integral que inhibe la predicción de la tendencia del error que se pretende conseguir. Por tanto, la señal de salida del PD1debe ser tratada adecuadamente. Además, un valor más alto tanto en el término proporcional como en el derivativo del controlador del ángulo de trayectoria PD2 mejoraría el error en régimen permanente y reduciría la sobreoscilación respectivamente. Sin embargo, incrementando el término derivativo de PD2 aumenta también el ruido a la salida de PD1 es decir, las figuras 7.7a y 7.7b tendrían más ruido. 7.2.3 Diseño de trayectoria para el aterrizaje Para la construcción de la trayectoria es esencial tener en cuenta varios factores: •El sistema evoluciona lentamente por lo que debe ser una trayectoria poco brusca y suave • Antes de la última ignición el vehículo debe estar vertical ( θ=90o ) para que cuando se produzca el encendido este no se desvíe del punto de aterrizaje y pierda precisión •Antes de la última ignición la velocidad horizontal debe de ser nula, es decir, γ=90o Tomando estos factores para la elaboración de la trayectoria esta tendrá lugar entre los 32000m y los 4000m de altura. La última ignición se activa alrededor de los 2100m, por lo que el sistema tiene un margen para alcanzar el régimen estacionario y alcanzar los 2100m con θ=γ=90o justo encima del punto de aterrizaje. Si se observa la figura 7.9, está compuesta por dos parábolas cuyos máximos se encuentran en los 32000m y 4000m respectivamente. Figura 7.9 Trayectoria a seguir (para una distancia de aterrizaje de 7m). Si definimos la altura inicial h0 , la altura final hf y la distancia horizontal final xf podemos definir la trayectoria de manera que, siendo la altura el eje z y la distancia horizontal, el eje x, x1(z) = az2+bz +c[x′ 1(h0) = 0x1(h0) = 0x1(hm) = xm](7.3) x2(z) = dz2+f z +g[x′ 2(hf) = 0x2(hf) = xfx2(hm) = xm](7.4) siendo hm y xm la altura y la distancia media respectivamente. Planteando un sistema de ecuaciones para calcular los términos de las parábolas se obtienen, 7.3 Resultados 55   2h01 0 h2 0h01 h2 mhm1   a b c =  0 0 xm (7.5)   2hf1 0 h2 fhf1 h2 mhm1   d f g =  0 xf xm  (7.6) Una vez obtenidos los coeficientes, el ángulo de trayectoria se calcula de la siguiente forma. γ=arctan(2az +b)h0>z>hm(7.7) γ=arctan(2dz +f)hm>z>hf(7.8) Para el caso de h0=32000 m , hf=4000 m y xf=7m el ángulo de trayectoria necesario será el representado en la figura 7.10. En las sucesivas simulaciones h0yhfse mantendrán con esos valores. Figura 7.10 Ángulo de trayectoria a seguir (para una distancia de aterrizaje de 7m). 7.3 Resultados Constituido el sistema de piloto automático con los controladores del ángulo de asiento y del ángulo de trayectoria y calculada la trayectoria en función de la distancia que se quiera alcanzar, es hora de llevar a cabo la simulación de la maniobra de aterrizaje. El objetivo de las simulaciones que se van a llevar a cabo es comprobar que el vehículo aterriza en el punto indicado a una velocidad inferior a 5km/h . Para ello en todo el recorrido estará activado el control de empuje y por debajo de los 33000m estará activado el control de actitud. La razón por la que no se activa antes es debida a que a mayores altitudes la atmósfera es muy tenue por lo que el guiado aerodinámico es poco efectivo. 7.3.1 Simulación dentro del rango operacional Se ha establecido una distancia límite por encima de la cual no se garantiza que la precisión con la que llega al punto de aterrizaje sea lo suficientemente baja. Dicha distancia límite es 14 m. La limitación es debida al hecho de que si antes de la última ignición el vehículo no cumple con θ=γ=90o (el vehículo está vertical y la velocidad horizontal es nula), cuando se encienda la planta propulsora este se desviará del punto deseado. Por lo tanto si las distancias no son muy grandes el vehículo podrá alcanzar el régimen permanente más fácilmente y por tanto tener la precisión suficiente. Para la simulación que se va a llevar a cabo, a posición inicial es ( x0=0m , y0=0m , z0=54000 m ), la actitud inicial ( θ0=90o , ψ0=0o , φ0=0o ) y la velocidad inicial ( u0=1389 m/s , v0=0m/s , w0=0m/s ). Son las mismas condiciones iniciales que se han utilizado anteriormente. 62 Capítulo 8. Conclusiones se caracteriza por tener una evolución lenta con inercias elevadas y una potencia de control limitada es necesario un control con carácter anticipativo, sin sobreoscilaciones. No obstante, en un afán de perseguir esas características se entra en conflicto con otra variable, el ruido. Así pues el ruido ha sido el factor limitante en la sintonización de los controladores. Pese a las dificultades encontradas, ha sido posible conseguir el aterrizaje de forma controlada del vehículo llegando a la plataforma de aterrizaje con un margen de error de alrededor del 1 % a una velocidad que ronda entre los 2km/h y los 6km/h dentro del rango operacional de 14mde distancia horizontal. Finalmente, este trabajo ha sacado a la luz la necesidad de concurrencia de muchos campos de la ingeniería. Por un lado las características aerodinámicas, propulsivas o estructurales, la mecánica del sólido rígido, mecánica orbital y por supuesto la ingeniería de sistemas. Para que en proyecto de esta magnitud salga adelante es esencial la concurrencia de todos los departamentos. En este estudio existirán carencias en algunos campos más que en otros e incluso la inexistencia de los mismos. A lo largo del trabajo se ha intentado abarcar los campos esenciales con la mayor profundidad posible dentro de las limitaciones propias de un trabajo de fin de carrera. En vista de los resultados, la búsqueda de un trabajo en línea con los sucesos reales ha sido llevada a cabo con éxito, sin dejar de lado su margen de mejora. 8.1 Trabajo futuro La versatilidad que proporciona la estructura del modelo abre muchas puertas para la implementación de nuevos estudios aerodinámicos, propulsivos o controles más óptimos y precisos. En primer lugar, en buena parte de la maniobra de aterrizaje, el vehículo trabaja en régimen supersónico lo cual tiene consecuencias aerodinámicas evidentes. Las ondas de choque generadas por el vuelo supersónico del vehículo generan una variación de presión, temperatura, densidad y velocidad del aire. Además, a la resistencia existente se le añade la resistencia de onda. Por otro lado la velocidad aerodinámica que experimentan los actuadores podría ser distinta a la del vehículo por las razones expuestas anteriormente. Un estudio aerodinámico futuro con CFDs implementando la geometría del vehículo y sus dispositivos de control proporcionará un modelo aerodinámico más preciso. Al igual que para el modelo aerodinámico, en vista al modelo propulsivo, existe un campo de estudio relativo al diseño de un motor cohete adecuado para la misión a realizar. Sin embargo, más allá del margen de mejora de los modelos, es interesante ver de qué manera se puede optimizar la maniobra de reentrada del vehículo hasta su punto de aterrizaje. De todas las técnicas globales de optimización, los métodos basados en inteligencia de enjambre (swarm intelligence, SI) se están volviendo cada vez más populares debido a su velocidad y precisión [38]. Este método está inspirado por fenómenos naturales como el comportamiento de una colonia de hormigas o una bandada de pájaros. El Particle Swarm Optimization (PSO) evita los cálculos que precisan las aproximaciones analíticas convencionales usando una solución existente para un problema específico e intentando encontrar otras posibles trayectorias para otros objetivos de interés. Este campo proporciona una gran variedad de métodos para la optimización del problema, la cual ha sido omitida en este trabajo. Incluso, con la integración de un nuevo modelo aerodinámico se podría simular el lanzamiento de un misil o de un cohete que daría lugar al estudio de trayectorias óptimas para problemas multifase [39], como lo puede ser el lanzamiento de un cohete. Asimismo, existen muchas opciones para el diseño del piloto automático. Si bien ha sido posible su control mediante controles sencillos como relé o proporcionales derivativos, cabe mencionar que el estudio de controles LPV (Linear parameter-varying control) [40], LQR (Linear-quadratic-regulator) de horizonte finito[37] o el control predictivo [36] es muy interesante dada la dinámica del problema. El uso de un controlador LPV para el control de la posición y actitud de un vehículo espacial con 6 grados de libertad no lineales presenta ventajas frente a un sistema lineal invariante con el tiempo. Un sistema LPV está definido como un sistema lineal cuya dinámica depende de una variable externa cuyos valores son desconocidos a priori pero pueden ser calculados a medida que el sistema avanza. Por otro lado, existen campos con los que no se ha contado por su complejidad o extensión. Un ejemplo de ello es el estudio estructural y de vibraciones. Los esfuerzos debido a la resistencia del aire, la diferencia de presión y temperatura y las aceleraciones a las que se ve sometido el vehículo a lo largo de la trayectoria, son parámetros a tener en cuenta a la hora tanto del diseño del cohete como la selección de materiales. Además el paso por el régimen transónico trae consigo un aumento de las vibraciones del vehículo por lo que habría que tener en cuenta que le frecuencia natural del cohete no coincide con la del sistema. Índice de Figuras 1.1 Cronograma carrera espacial 2 1.2 Coste de lanzamiento. Fuente (editada):[7] 3 2.1 Lanzadores actuales. Fuente:[17] 7 2.2 CubeSats. Fuente:[19] 8 3.1 Falcon 9 [21] 9 3.2 Ejemplos de distintos mallados en grid fins. Fuente [22] 10 3.3 Trayectoria de reentrada del Falcon 9. Fuente: [21] 11 3.4 Modelo de atmósfera exponencial 12 3.5 Temperatura-Altitud 13 3.6 Coeficiente de resistencia de la base en función de Mach 14 3.7 14 3.8 Coeficiente de resistencia del grid fin en función de Mach 15 3.9 Momentos en ejes cuerpo 15 3.10 Distribución de la planta motora 16 3.11 Empuje-Altitud 18 4.1 Sistema de ejes horizonte local y ejes viento I 20 4.2 Sistema de ejes horizonte local y ejes viento II 20 4.3 Vista de planta del vehículo (a escala) 21 4.4 Vista de alzado del vehículo (a escala) 22 4.5 Sistema de ejes viento y ejes cuerpo 22 4.6 Funcionamiento de los grid fins (I) 28 4.7 Funcionamiento de los grid fins (II) 29 5.1 Diagrama de bloques 33 5.2 Fuerzas y torques 34 5.3 Ecuaciones del movimiento 34 5.4 Posición y actitud 35 5.5 Diagrama de flujo para conversión a ángulos de Euler 35 5.6 37 5.7 Evolución en caída libre 37 5.8 Actitud 38 5.9 Comparación 39 5.10 Deflexión de la superficie de control 1 40 5.11 Evolución 40 5.12 Ángulos 41 5.13 Deflexión de la superficie de control 4 41 5.14 Evolución 42 5.15 Ángulos 42 5.16 Evolución 43 63 64 Índice de Figuras 5.17 Actitud 43 6.1 Diagrama de bloques del control de empuje 45 6.2 Diagrama de bloques del control de empuje 46 6.3 Arquitectura del relé 46 6.4 Control del empuje 47 6.5 Evolución de la velocidad 47 7.1 Ángulo de trayectoria γy ángulo de asiento θ50 7.2 Arquitectura del piloto automático 50 7.3 Criterio de signos 51 7.4 51 7.5 Respuesta del controlador para el ángulo de asiento 52 7.6 Respuesta del controlador para el ángulo de trayectoria 53 7.7 Evolución de θ53 7.8 Actuadores 53 7.9 Trayectoria a seguir (para una distancia de aterrizaje de 7m) 54 7.10 Ángulo de trayectoria a seguir (para una distancia de aterrizaje de 7m) 55 7.11 Resultados de la simulación 56 7.12 Actuadores 56 7.13 Comportamiento de PD257 7.14 Comportamiento de PD157 7.15 Resultados de la simulación 58 7.16 Actuadores 58 7.17 Comportamiento de PD259 7.18 Comportamiento de PD159 Índice de Tablas 3.1 Parámetros del motor cohete 11 3.2 Atmósfera Estándar Internacional 13 3.3 Parámetros del motor cohete 17 65 Bibliografía [1] P. G. Martín, «Julio Verne, el escritor visionario,» National Geographic, 2019. dirección: https : // historia.nationalgeographic.com.es/ a/ julio-verne-escritor-visionario_13488. [2] J. M. Sadurní, «Wernher Von Braun, el artífice de la conquista a la Luna,» National Geographic, 2020. dirección: https:// historia.nationalgeographic.com.es/ a/ wernher-von-braun-artifice-conquista- luna_15168. [3] H. Rodríguez, «La carrera espacial paso a paso,» National Geographic, 2019. dirección: https:// www. nationalgeographic.com.es/ llegada-del-hombre-a-la-luna/ carrera-espacial-paso-a-paso_14369. [4] R. V. Valenzuela, Transparencias de Mecánica Orbital y Vehículos Espaciales. Escuela Técnica Superior de Ingeniería, 2021. [5] «BEIDOU naviagation and positioning satellites,» China National Space Administration (CNSA), dirección: http:// www.cnsa.gov.cn/ english/ n6465715/ n6465716/ c6480954/ content.html. [6] «Indian Regional Navigation Satellite System (IRNSS) : NavIC,» Indian Space Research Organisation (ISRO), dirección: https:// www.isro.gov.in/ irnss-programme. [7] B. Venditti. «The Cost of Space Flight Before and After Space X.» (), dirección: https:/ / www. visualcapitalist.com/ the-cost-of-space-flight/ . [8] «Space Launch System,» NASA, dirección: https:// ntrs.nasa.gov/ api/ citations/ 20160014576/ downloads/ 20160014576.pdf . [9] «Vega-C,» ESA, dirección: https:// www.esa.int/ Enabling_Support/ Space_Transportation/ Launch_ vehicles/ Vega-C. [10] «Space Rider,» ESA, dirección: https:// www.esa.int/ Enabling_Support/ Space_Transportation/ Space_Rider. [11] «Soyuz-2,» ROSCOSMOS, dirección: http:// en.roscosmos.ru/ 468/ . [12] «Falcon Heavy,» Space X, dirección: https:// www.spacex.com/ vehicles/ falcon-heavy/ . [13] «Dragon,» Space X, dirección: https:// www.spacex.com/ vehicles/ dragon/ . [14] «Starship,» Space X, dirección: https:// www.spacex.com/ vehicles/ starship/ . [15] «dearMoon,» dirección: https:// dearmoon.earth/ . [16] «New Shepard,» Blue Origin, dirección: https:// www.blueorigin.com/ new-shepard/ . [17] «New Glenn,» Blue Origin, dirección: https:// www.blueorigin.com/ new-glenn. [18] «Blue Moon,» Blue Origin, dirección: https:// www.blueorigin.com/ blue-moon/ lunar-transport. [19] «CubeSats,» ESA, dirección: https:// www.esa.int/ Enabling_Support/ Preparing_for_the_Future/ Discovery_and_Preparation/ CubeSats. [20] «Payload Mass to Orbit,» BryceTech, [21] Falcon user’s guide, SpaceX, 2020. dirección: https:// www.spacex.com/ vehicles/ falcon-9/ . [22] K. S. Orthner, «Aerodynamic Analysis of Lattice Grid Fins in Transonic Flow,» AIR FORCE INSTITUTE OF TECHNOLOGY, 2004. dirección: https:// scholar.afit.edu/ etd/ 3916. 67 68 Bibliografía [23] «Preguntas frecuentes sobre la reentrada de Tiangong 1,» ESA, dirección: https:// www.esa.int/ Space_ in_Member_States/ Spain/ Preguntas_frecuentes_sobre_la_reentrada_de_Tiangong_1. [24] P. Koteswaram, «The Easterly Jet Stream in the Tropics,» Tellus, págs. 43 - 57, 1958. dirección: https: // doi.org/ 10.3402/ tellusa.v10i1.9220. [25] A. Tewari, Atmospheric and Space Flight Dynamics. Modeling and Simulation with MATLAB and Simulink. Birkhauser, 2007, págs. 223-224, isbn: 0-8176-4437-7. [26] J. C. McDowell, «The edge of space: Revisiting the Karman Line,» Actra Astronáutica, vol. 151, págs. 668 - 677, 2018. dirección: https:// www.sciencedirect.com/ science/ article/ pii/ S0094576518308221. [27] R. Q. Campos, «Planificación y guiado óptimos de la trayectoria de ascenso a órbita de un vehículo lanzador sometida a perturbaciones y errores de modelado,» 2021, págs. 15 - 16. dirección: https: // idus.us.es/ handle/ 11441/ 128180. [28] «Earth Atmosphere Model - Metric Units,» NASA, 2021. dirección: https:// www.grc.nasa.gov/ www/ k- 12/ airplane/ atmosmet.html. [29] J. S. R. George E. Hansche, «Air Drag on Cubes at Mach Numbers 0.5 to 3.5,» American Institute of Aeronautics and Astronautics, págs. 83-84, dirección: https:// arc.aiaa.org/ doi/ 10.2514/ 8.2166. [30] M. D. Yan Zeng Jinsheng Cai, «Numerical Study on Drag Reduction for Grid-Fin Configurations,» American Institute of Aeronautics and Astronautics, 2009. dirección: https:// arc.aiaa.org/ doi/ abs/ 10. 2514/ 6.2009-1105. [31] J. M. Macha, «Drag of Circular Cylinders at Transonic Mach Numbers,» American Institute of Aeronautics and Astronautics, vol. 14, n.o6, pág. 606, 1977. [32] D. R. Rivas, Mecánica de Vuelo. Escuela Técnica Superior de Ingeniería, 2016, págs. 5-14. [33] R. V. Valenzuela, Fundamentos de Navegación Aérea, ép. Sistemas de referencia. Actitud. Ecuaciones de la navegación. Escuela Técnica Superior de Ingeniería, 2020, págs. 23-49. [34] A. Markelov, «Uso de Cuaterniones para Representar Rotaciones,» Repositorio Institucional de la Universidad Tecnológica de Panamá, págs. 1 - 2, 2022. dirección: http:// revistas.utp.ac.pa/ index.php/ tecnologia-hoy/ article/ view/ 628. [35] «SpaceX Launch: OTV5 / X-37B : SpaceX Landing of first stage Falcon 9.» (2018), dirección: https:// www.youtube.com/ watch?v=_eEPOYWD4WU&t=559s. [36] M. Morari, C. E. Garcia y D. M. Prett, «Model predictive control: Theory and practice,» IFAC Proceedings Volumes, vol. 21, n.o4, págs. 1-12, 1988. [37] A. Grancharova y T. A. Johansen, «Survey of explicit approaches to constrained optimal control,» Lecture Notes in Computer Science, vol. 3355, págs. 47-97, 2005. [38] A. Rahimi, K. D. Kumar y H. Alighanbari, «Particle Swarm Optimization Applied to Spacecraft Reentry Trajectory,» Journal of Guidance, Control and Dynamics, vol. 36, n. o 1, 2013. dirección: https:// arc.aiaa.org/ doi/ 10.2514/ 1.56387. [39] M. A. Patterson y A. V. Rao, «A MATLAB Software for Solving Multiple-Phase Optimal Control Problems Using hp-Adaptive Gaussian Quadrature Collocation Methods and Sparse Nonlinear Programming,» University of Florida, vol. 41, n. o 1, 2014. dirección: http:// dx.doi.org/ 10.1145/ 2558904. [40] A. S. Ghersin y R. S. S. Peña, «LPV Control of a 6-DOF Vehicle,» Universidad de Sevilla, 2002. dirección: https:// ieeexplore.ieee.org/ document/ 1058058. [41] «Updates,» Space X, dirección: https:// www.spacex.com/ updates/ #starship-moon-announcement. [42] «About Blue Origin,» Blue Origin, dirección: https:// www.blueorigin.com/ about-blue. [43] «Opening Space for good,» Virgin Orbit, dirección: https:// www.virgin.com/ virgin-companies/ virginorbit. [44] «Virgin Galactic,» Virgin Galactic, dirección: https:// www.virgingalactic.com/ . [45] J. S. Shamma y J. R. Cloutier, «Gain-Scheduled Missile Autopilot Design Using Linear Parameter Varying Transformations,» Journal of Guidance, Control and Dynamics, vol. 16, n. o 2, 1993. dirección: https:// arc.aiaa.org/ doi/ abs/ 10.2514/ 3.20997.