Full text
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 Trabajo Fin de Grado Grado en Ingeniería Aeroespacial Control de actitud de las E-sail mediante la modulación del voltaje de los tethers Autor: Juan Pérez Gómez Tutor: Guillermo Pacheco Ramos y Rafael Vázquez Valenzuela Dpto. de Ingeniería Aeroespacial y Mecánica de Fluidos Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2022
Trabajo Fin de Grado Grado en Ingeniería Aeroespacial Control de actitud de las E-sail mediante la modulación del voltaje de los tethers Autor: Juan Pérez Gómez Tutor: Guillermo Pacheco Ramos y Rafael Vázquez Valenzuela Profesor Sustituto Interino y Profesor Titular Dpto. de Ingeniería Aeroespacial y Mecánica de Fluidos Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2022
Trabajo Fin de Grado: Control de actitud de las E-sail mediante la modulación del voltaje de los tethers Autor: Juan Pérez Gómez Tutor: Guillermo Pacheco Ramos y Rafael Vázquez Valenzuela 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 Tr as 4 años de duro esfuerzo y trabajo esta etapa llega a su fin. Han sido 4 años duros en los que he pasado por momentos muy difíciles pero de los he conseguido sobreponerme y sobre todo aprender. No quería cerrar esta etapa sin agradecer a todas las personas que me han acompañado en el camino. Los pilares de mi vida, mis padres, sin sus esfuerzos y sacrificios hoy no estaría en la situación en la que estoy. Por su apoyo incondicional en los peores momentos y por haberme dado la educación que tengo, ellos siempre han estado ahí. Sin ellos, hoy no sería nada. También me gustaría acordarme de mi hermana, Lucía gracias por estar ahí siempre y por sacarme una sonrisa en lo momentos de estrés, eres una persona fundamental en mi vida y sin ti no habría logrado llegar hasta aquí. A mis abuelos, mis segundos padres, por su amor y su apoyo incondicional. Abuela allá donde estés, quiero decirte que lo he conseguido y que todo el esfuerzo ha merecido la pena, sigue siendo la luz que me guía desde arriba. Quería acordarme también de mis amigos, por su constante apoyo y por animarme en los momentos duros. No quería irme sin agradecer a Rafa y a Guillermo toda su ayuda. A pesar de que el comienzo de este proyecto no fue fácil, ellos han estado siempre ahí para ofrecer cualquier ayuda. Gracias por hacer que aún me apasione más este mundo. Por todas aquellas personas que se han cruzado en mi camino durante estos años. De todas he aprendido algo y me han ayudado a crecer como persona, solo tengo palabras de agradecimiento para ellas. Juan Pérez Gómez Sevilla, 2022 I
Resumen En esta memoria se abordará el tema de control de actitud de las velas solares eléctricas, un nuevo mecanismo de generación de empuje que usa los iones procedentes del plasma del viento solar y que permite, por tanto, eliminar la masa de propulsante necesaria para viajar por el Sistema Solar. Se comenzará realizando un estudio del estado actual de las velas solares eléctricas, sus posibles aplicaciones y sus ventajas respecto a sistemas cuyo fin es el mismo. Posteriormente, se introducirán aspectos básicos de la representación y la cinemática de la actitud, para luego, basándonos en la literatura de referencia presentar el modelo dinámico del que se ha hecho uso. Además, se reproducirán los resultados obtenidos aplicando la ley de control expuesta en [ 2 ] con el objetivo de comprender mejor el comportamiento de las velas solares eléctricas. Asimismo, se propondrá una nueva ley de control basada en el Control LQR con el objetivo de mejorar los resultados anteriormente comentados. Por ultimo, se analizarán los resultados obtenidos en varios casos en los que se han modificado distintos parámetros del sistema, para ver como se comporta la E-sail con la nueva ley de control propuesta. II
Abstract Th is report will address the issue of attitude control of electric solar sails, a new oush generation mechanism that uses ions from solar wind plasma and therefore allows to eliminate the mass of propellant necessary to travel through the Solar System. A study of the current state of electric solar sails, their possible applications and their advantages over systems whose purpose is the same is realised in the first place. Subsequently, basic aspects of the representation and kinematics of attitude are introduced. Then, the dynamic model thata has been used and that has benn obtained from the reference literature is presented. In addition, the results obtained by applying the control law set out in [ 2 ] will be reproduced in order to better understand the behaviour of electric solar sails. Likewise, a new control law based on LQR Control will be proposed with the aim of improving the results mentioned above. Finally, the results obtained in several cases in which different parameters of the system have been modified will be analysed, to see how the E-sail behaves with the proposed new control law. III
2Índice de Figuras 4.6 Evolución temporal de αn al aplicar la ley de control (4.52) en la fase de tracking y posteriormente la ley de control (4.47) en la fase de regulación 53 4.7 Evolución temporal de φ y θ al aplicar la ley de control (4.52) en la fase de tracking y posteriormente la ley de control (4.47) en la fase de regulación 53 4.8 Evolución temporal de ψ al aplicar la ley de control (4.52) en la fase de tracking y posteriormente la ley de control (4.47) en la fase de regulación 54 4.9 Evolución temporal de las velocidades angulares al aplicar la ley de control (4.52) en la fase de tracking y posteriormente la ley de control (4.47) en la fase de regulación 55 4.10 Evolución temporal de las variables de control al aplicar la ley de control (4.52) en la fase de tracking y posteriormente la ley de control (4.47) en la fase de regulación. 55 4.11 Comparativa de las distintas evoluciones que se tienen para diferentes valores de referencia de αn durante la fase de tracking 56 4.12 Escalón producido en la evolución de las variables de control al pasar de la ley de control (4.52) en régimen transitorio a la ley de control (4.47) en régimen permanente 57 4.13 Evolución temporal de αn al aplicar la ley de control (4.52) en la fase de tracking y posteriormente la ley de control (4.47) en la fase de regulación, añadiendo además el filtro de primer orden 58 4.14 Evolución temporal de θ y φ al aplicar la ley de control (4.52) en la fase de tracking y posteriormente la ley de control (4.47) en la fase de regulación, añadiendo además el filtro de primer orden 58 4.15 Evolución temporal de ψ al aplicar la ley de control (4.52) en la fase de tracking y posteriormente la ley de control (4.47) en la fase de regulación, añadiendo además el filtro de primer orden 59 4.16 Evolución temporal de las velocidades angulares al aplicar la ley de control (4.52) en la fase de tracking y posteriormente la ley de control (4.47) en la fase de regulación, añadiendo además el filtro de primer orden 59 4.17 Evolución temporal de las variables de control al aplicar la ley de control (4.52) en la fase de tracking y posteriormente la ley de control (4.47) en la fase de regulación, añadiendo además el filtro de primer orden 60 4.18 Evolución de las variables de control en las proximidades de la transición de la ley de control (4.52) en la fase de tracking a la ley de control (4.47) en la fase de regulación 60 5.1 Evolución de αnvariando el número de tethers 62 5.2 Evolución de φy de θvariando el número de tethers 63 5.3 Evolución de ψvariando el número de tethers 64 5.4 Evolución de las velocidades angulares variando el número de tethers 65 5.5 Evolución de las variables de control variando el número de tethers 66 5.6 Evolución de αnvariando el valor de ¯ αn f 68 5.7 Evolución de φyθvariando el valor de ¯ αn f 69 5.8 Evolución de ψvariando el valor de ¯ αn f 70 5.9 Evolución de las velocidades angulares variando el valor de ¯ αn f 71 5.10 Evolución de las variables de control variando el valor de ¯ αn f 72 5.11 Evolución de αnvariando el valor de tf73 5.12 Evolución de φyθvariando el valor de tf74 5.13 Evolución de ψvariando el valor de tf75 5.14 Evolución de las velocidades angulares variando el valor de tf76 5.15 Evolución de las variables de control variando el valor de tf77
Índice de Tablas 2.1 Parámetros de los diferentes subsistemas, con a⊕=1mm/s210 2.2 Valores de parámetros geométricos y funcionales de una E-sail. 12 3.1 Definición de los valores de los distintos parámetros de la E-sail 29 3
1 Introducción 1.1 El concepto de las E-sail La E-sail, es una idea que surge en el año 2004, de la mano del físico espacial, astrobiólogo e inventor finlandés Pekka Janhunen [ 1 ]. Esta idea se basa en las velas magnéticas que se encuentran en una magnetosfera artificial y que fueron introducidas por Zubrin y Andrews en 1991, mediante dicha magnetosfera, el viento solar es deflectado y esto permite extraer cantidad de movimiento de él. La adaptación de estas velas magnéticas, hace que surjan las E-sail. Las velas eléctricas (E-sail), consisten en un número de cables conductores de gran longitud (tethers) dispuestos alrededor de un eje de simetría y que están cargados positivamente. Dichos cables de gran longitud son desplegados mediante el giro del vehículo espacial alrededor de su eje de simetría, que permite además que estos cables se mantengan extendidos. Esencialmente, estas velas usan la presión dinámica del viento solar para generar una pequeña, pero continua cantidad de empuje. La principal causa de por qué esta tecnología es interesante, reside en que permite eliminar la mayor parte de masa de propelente que debe cargar el vehículo espacial para realizar su misión, la que puede suponer una restricción, en cuánto a peso, a la hora de realizar ciertas misiones. Además, la potencia consumida para producir 1 N de empuje es mucho menor que en la mayoría de sistemas de propulsión existentes, cómo se verá más adelante en la sección 2.4. 1.2 Principio del funcionamiento de las E-sail Las velas solares eléctricas son un nuevo método de propulsión que, usa la presión dinámica del viento solar como fuente de generación de empuje, fundamento por el cual este sistema propulsivo es eficiente desde el punto de vista de consumo de combustible, cómo ya se ha mencionado anteriormente. La base de su funcionamiento, es que las partículas ionizadas del viento solar son repelidas, mediante la Ley de Coulomb, debido a que los tethers que componen la vela están cargados positivamente. Esto modifica las trayectorias de los protones pertenecientes al viento solar, y es dicha modificación en la trayectoria de los protones la que permite la extracción de cantidad de movimiento del viento solar. La ganancia de cantidad de movimiento por parte de los tethers será la empleada para obtener empuje, el cual se transmitirá mecánicamente al vehículo espacial y permitirá el movimiento de este por el Sistema Solar. Debido a que los tethers se encuentran a una diferencia de potencial positiva, los electrones del plasma del viento solar son atraídos por esa carga positiva, esto provoca que la diferencia de potencial antes mencionada se neutralice. Es por ello, que los electrones que quedan atrapados en los conductores se deben desechar, puesto que penalizan al empuje, magnitud cuya relación con la diferencia de potencial es directa. La eliminación de estos, se lleva a cabo con un cañón de electrones que funciona gracias a unos paneles solares u otros sistemas de generación de potencia eléctrica. En el caso en el que se quiera que la vela eléctrica deje de proporcionar empuje solo bastaría con desactivar el cañón de electrones, de tal forma que la diferencia de potencial en los tethers se neutralizaría [ 7 ]. También, se ha analizado que la ausencia de dichos electrones que quedan atrapados hace que el empuje sea cinco veces superior con respecto a la situación en la que sí hay, cómo se trata en la referencia [ 3 ]. Además, la diferencia de tensión que se crea en los cables conductores, no solo hace que se genere empuje, sino que también permite controlar la dirección en la que se mueve el vehículo propulsado por la vela solar eléctrica mediante la modulación de la caída de tensión en cada uno de los tethers que la componen. 4
1.3 Estado actual de las E-sail 5 Figura 1.1 Esquema conceptual del funcionamiento de una E-sail. Extraída de [2]. Para que sirva de referencia, a una distancia de 1 AU del Sol, la presión dinámica del viento solar es de 2 nPa. Por lo que si se tiene un vehículo espacial, equipado con este sistema de propulsión, cuya masa sea de 100 kg, con unos 100 tethers de 10 km de longitud cada uno, se puede conseguir que este sufra una aceleración de 1 mm/ s2 aproximadamente y, por tanto, un empuje de 0,1 N, tal y cómo se indica en [ 7 ]. Es por ello, que este nuevo sistema de propulsión es tan prometedor para misiones futuras. En la figura 1.1 , extraída de [ 2 ], se puede ver un esquema conceptual de este nuevo mecanismo de generación de empuje. En dicha figura se puede ver el cañón de electrones que permite la eliminación de los mismos como se ha dicho antes y el conjunto de cables que conforman la vela solar eléctrica y que se encuentran cargados. 1.3 Estado actual de las E-sail Las E-sail se trata de una tecnología que ha surgido hace relativamente poco tiempo (2004), es por ello que no se ha usado en demasiadas misiones. La misión pionera en usar este tipo de tecnología ha sido la ESTCube-1, es por ello, que es en la que nos basaremos para comentar el estado actual de las velas solares eléctricas. Se trata de un proyecto de estudiantes de la Universidad de Taru, Estonia, en colaboración con la Agencia Espacial Europea. El principal fin de dicha misión era el de colocar en órbita el primer satélite que usase una vela solar eléctrica como mecanismo de generación de empuje, y comprobar el funcionamiento de la misma. El ESTCube-1 se trata de un satélite en forma de cubo de pequeñas dimensiones (10 × 10 × 10 cm) y una masa ligeramente superior a 1 kg. Dicho satélite esta conformado por el sistema de potencia eléctrica (EPS), el sistema de comunicación (COM), el sistema de control de datos y comandos (CDHS), el sistema de determinación y control de actitud (ADCS) y la carga asociada a la E-sail. Dicho satélite fue lanzado el 7 de mayo de 2013, y durante su desarrollo, uno de los retos más difíciles de enfrentar fue la integración e instalación de los cables y mazos eléctricos de los que constaba el satélite, debido a su pequeño tamaño. Con respecto, a la vela solar que lo propulsaba, hay que decir que estaba compuesta por 15 m de cable en total, una fuente de alta tensión al final de cada tether para cargarlo eléctricamente y un tether secundario (se ve en más detalle en la Sección 2.1) que conectaba las fuentes de alta tensión de cada tether principal. El principal problema de la misión fue que al llevar a cabo el despliegue de los tethers, este o fue exitoso dado que se produjo un fallo en las medidas de diagnóstico, la razón exacta de este fallo siguen sin conocerse. Según la referencia [ 11 ] uno de las posibles mejoras para las E-sail sería la de instalar sensores que detecten el comportamiento y la situación en cada instante de los distintos componentes de la vela. Lo que nos hace ver que si bien es una tecnología muy prometedora porque permite eliminar la masa de propelente del vehículo espacial, aún no está lo suficientemente madura como para poder implementarla en misiones actuales. A pesar de ello, son muchos los estudios y ensayos en laboratorios lo que se están llevando a cabo para lograr
6Capítulo 1. Introducción Figura 1.2 Satélite ESTCube-1. Extraída de [11]. solventar los problemas encontrados en cuanto a la complejidad mecánica y estructural del sistema, dado el gran número de cables con los que cuenta la E-sail. 1.4 Objetivo y alcance del proyecto Una vez definido que son las E-sail y cuál es el mecanismo básico de su funcionamiento, el objetivo de este proyecto es el de modelar una ley de control que permita controlar la actitud de este nuevo método de propulsión mediante la modulación del voltaje al cual están sometidos los tethers. Para ello, se parte de un modelo dinámico obtenido de la literatura que se ha usado de referencia y mediante control LQR se halla esta nueva ley de control que mejore los resultados mostrados en la literatura en la que se basa el modelo dinámico que se ha tomado. Con ello se pretende realizar una maniobra de actitud y posteriormente mantener la actitud alcanzada para una órbita prefijada, sin que se tenga en cuenta el control de la órbita, es por ello que como se verá más adelante el horizonte de tiempo es pequeño. Por último, comentar que para la resolución del problema planteado se despreciarán todas las perturbaciones a las que se encuentran sometidos vela y vehículo y que únicamente se considerará la perturbación provocada por los propios tethers. 1.5 Estructura del documento Tras presentar cuál es el objetivo que se pretende alcanzar con este proyecto, se prosigue definiendo la estructura, en cuanto a contenido, que seguirá el documento. En primer lugar, en el Capítulo 2 se hace una descripción de las E-sail y de sus aplicaciones. En este estudio se presentan los principales componentes de una velar solar eléctrica y las posibles misiones a realizar con este mecanismo de propulsión. Asimismo, se muestran los distintos modelos dinámicos que han ido surgiendo a lo largo de estos años y durante el desarrollo del mismo se irán comentando las ventajas e inconvenientes que este sistema presenta con respecto a otros cuyos fines son similares. Posteriormente, en el Capítulo 3 se muestran los sistemas de referencia de los que se hace uso, así como conceptos básicos pero esenciales para representar la actitud de un vehículo espacial y conocer su cinemática. Además, se presenta el modelo por el que se ha optado a la hora de representar la dinámica del sistema. Dicho modelo incluye las fuerzas que intervienen en el problema, todas las expresiones analíticas de los pares de fuerza que actúan sobre la vela y el modelo de deformada de los cables de la E-sail. Por otro lado, se analizan los resultados obtenidos mediante simulaciones numéricas al aplicar la ley de control propuesta en la referencia [2] para el control de la actitud de una vela solar eléctrica.
1.5 Estructura del documento 7 Luego en el Capítulo 4, se presenta una nueva ley de control que mejora la propuesta por el artículo que se ha usado como referencia para este proyecto ([ 2 ]). Para ello se usa el control LQR como se mencionó en la sección 1.4, pero previamente se llevará un proceso de adimensionalización y de linealización de las variables y ecuaciones que componen el problema. Por otro lado, en el Capítulo 5 se analizarán una variedad de casos en los que se modifican distintos parámetros del problema para ver como se comportan las leyes de control obtenidas en el Capítulo 4 ante dichas modificaciones. Por último, en el Capítulo 6 se presentarán las conclusiones obtenidas tras la realización del presente trabajo y posibles líneas de trabajo futuro con las que poder ampliarlo.
2 Descripción de las E-sail y sus aplicaciones Es te capítulo tiene cómo objetivo mostrar cuál es el estado actual de las E-sail. Para ello, primeramente se definen cuáles son los componentes que constituyen una vela solar eléctrica. Por otro lado, también se realiza un análisis de los órdenes de magnitud, en cuanto a consumo y empuje que se pueden tener en éste nuevo sistema propulsivo. Asimismo, se muestran cuáles pueden ser los campos de aplicación de este novedoso mecanismo de generación de empuje y las consecuencias de su uso. Adicionalmente, se muestra la variedad de modelos dinámicos que han surgido a lo largo de estos últimos años y sus características. 2.1 Principales componentes de las E-sail Se muestra, de una manera descriptiva, cuáles son los principales subsistemas y componentes necesarios a la hora de integrar la vela solar eléctrica en un vehículo espacial y el impacto que tiene esta sobre dichos subsistemas. Nos centraremos sobre todos en las estimaciones de masa, pues son las que definen las necesidades de empuje para la puesta en órbita del vehículo espacial. Todas las deducciones que aquí se muestran han sido basadas en las que aparecen en [4]. Se muestra en la figura 2.1 un esquema simplificado, a modo descriptivo, de los principales componentes que configuran una E-sail. La imagen ha sido extraída de [4]. Comenzaremos por los tethers principales (main tethers) que se tratan de unos cables de elevada longitud y que están cargados positivamente. Cada uno de ellos está formados por cuatro cables de un material conductor como puede ser el aluminio, donde uno de ellos es recto y es el cable base, mientras que los otros tres son se disponen formando “loops”. Generalmente, sus longitudes son bastante elevadas, siendo 2 km un valor bastante habitual en la longitud de un solo tether. Su masa dependerá linealmente de la longitud que estos tengan, considerando que λt p =1,155 ×10−5 kg/m es la masa por unidad de longitud del tether principal, se tiene que su masa será mt p =λt p ·L. El carrete ensamblador del tether principal (main tether reel assembly) se trata de un mecanismo motorizado que mantiene el tether enrollado en su interior y lo despliega cuando el vehículo espacial se encuentra en órbita. Por otro lado, el tether auxiliar (auxiliary tether) se usa para conectar las RUs (unidades remotas) con el objetivo de evitar colisiones entre tethers adyacentes. Su longitud es aproximadamente igual a la longitud de una circunferencia de radio L , donde L es la longitud del tether principal. Mientras que su masa se obtendría de manera análoga al caso del tether principal, pero con la longitud y la masa por unidad de longitud correspondientes. Con respecto a dicho componente decir que algunas propuestas lo consideran mientras que otras no. Para el caso que nos ocupa con el fin de simplificar no se ha considerado dicho tether. Otro elemento a considerar es la unidad remota (RU) ya mencionada, cuya función es la de controlar el giro de la plataforma que componen los tethers. Este componente puede contener impulsores que proporcionen el impulso necesario para producir el giro inicial o para producir giros pequeños durante operaciones en vuelo. Adicionalmente, los tethers necesitarán controladores y cámaras, estas últimas se encuentran dispuestas a lo largo del perímetro del vehículo espacial y su misión es la de localizar la posición de cada RU en el tether principal. 8
2.1 Principales componentes de las E-sail 9 Figura 2.1 Componentes de la E-sail. Extraída de [4]. Cabe hablar también, del subsistema de generación de potencia. Este está formado por paneles solares, que cuentan con sus respectivos mecanismos de despliegues; una unidad encargada de procesar la potencia que produce el voltaje del bus y la gran mayoría de velas incorporan también un paquete de baterías. Su función es la de proveer de potencia suficiente tanto a la carga de pago y al sistema de telemetría, que generalmente estarán compuestos por equipo científico; como a la E-sail, pues la potencia aquí es necesaria para que el cañón de electrones pueda funcionar. En general, unos valores típicos de potencia exigida por estos sistemas suelen ser 1 W/kg para la carga de pago y el sistema de telemetría, más 10 W de potencia base para la E-sail, todo esto referido a distancias entre 0,9 AU y 4 AU. Estos requerimientos de potencia determinarán la masa que tendrá el subsistema. En general, el sistema de telemetría se considera conjuntamente con la carga de pago, esto es debido a que cada instrumento científico genera datos que deben ser transmitidos por el sistema de telemetría. Además, las E-sail requieren de un sistema de actitud (ACS) cuyo objetivo es apuntar el eje de giro hacia el Sol, iniciando así el movimiento de giro al inicio necesario para el despliegue de la E-sail. Aunque la mayor parte de la cantidad de movimiento al inicio es suministrada por los impulsores de las RUs, una pequeña parte es propocionada por el ACS gracias a que este orienta el eje de giro adecuadamente. En la tabla 2.1 se muestran las estimaciones de las masas de los distintos subsistemas y algunos parámetros de interés, para una aceleración característica ( a⊕ , que es la aceleración máxima propulsiva que se tiene a la distacia de referencia, r⊕) de 1mm/s2. Con esta tabla, se hace ver cuáles son los órdenes de magnitud con los que se trabaja y cuáles serán los requerimientos de empuje que la E-sail deberá satisfacer para poder propulsar el vehículo espacial y todos los subsistemas que lo componen. Estos resultados no son sencillos de obtener, pues hay que tener en cuenta
10 Capítulo 2. Descripción de las E-sail y sus aplicaciones la integración entre los distintos subsistemas que componen el vehículo y las implicaciones que esto conlleva a la hora de diseñarlos, ya que cada subsistema impone unos requisitos sobre los demás. Tabla 2.1 Parámetros de los diferentes subsistemas, con a⊕=1mm/s2. Carga de pago (kg) 100 200 300 Número de tethers (N) 44 62 86 Longitud del tether L (km) 15,3 19,4 20 Empuje E-sail a 1 AU (N) 0,39 0,70 1,00 Potencia E-sail (W) 409 720 1026 Potencia con carga de pago no operando (W) 12,5 25,0 37,5 Potencia con carga de pago operando (W) 125 250 375 Masa tethers principales (kg) 7,79 13,9 19,8 Masa carrete tether principal (kg) 9,21 14,8 20,8 Masa cañón de electrones (kg) 1,20 2,13 3,05 Masa cámaras y computador (kg) 1,48 1,48 1,48 Masa RUs (kg) 43,1 59,3 77,6 Masa tether aux. (kg) 26,1 32,9 34,0 Masa sist. de Potencia (kg) 32,0 61,6 91,3 Masa sist. de Telemetría (kg) 25,0 50,0 75,0 Masa ACS (kg) 13,8 24,7 35,3 Masa estructura (kg) 65,2 116 166 Masa total sin E-sail (kg) 248 490 732 Masa efectiva E-sail (kg) 143 206 264 Masa Total (kg) 391 696 996 2.2 Estimación simplificada del empuje En esta sección se muestra a efectos descriptivos los órdenes de magnitud de los principales parámetros que podemos encontrar en una E-sail. Para ello, se considerará un modelo muy simple consistente en una malla compuesta de cables conductores, que representan a los tethers, tal y como se ha hecho en [ 1 ], referencia en la que se basan todas las deducciones que aquí se muestran. Dicha malla se dispone de tal forma que el viento solar la atraviesa, como se muestra en la figura Asumiendo que r=5×10−6 m, el área de la sección transversal de los cables que componen la malla será, Atrans =πr2(2.1) Teniendo en cuenta que la malla es cuadrada para simplificar cálculos y que el espacio entre cables es △x , el número total de cables será N=L/△x(2.2) donde L es la longitud lateral de la malla, suponiendo unos valores de L=35 km y de △x=5 km se obtiene que la longitud total de cable en la estructura es 2NL =2L2/△x=490 km (2.3) y su masa será m=2L2Aρ/△x=342,512 kg (2.4) donde se ha asumido que la densidad del material de los alambres conductores es ρ=8900 km/m 3 . Una vez calculados estos parámetros básicos, se procede a calcular la fuerza que actúa sobre la malla, esta fuerza vendrá dada por la presión dinámica que el viento solar ejerce sobre la estructura,
2.2 Estimación simplificada del empuje 11 Figura 2.2 Esquema simplificado de una E-sail. Extraído de [1]. pdin =ρvsv2(2.5) siendo ρvs =n·mp la densidad del viento solar que viene dada por el producto entre la densidad de iones del viento solar ( n [m −3 ]) y la masa del protón ( mp [kg]), mientras que vvs =400 km/s, es la velocidad del viento solar que se supone constante. Todos estos datos se han obtenido asumiendo una distancia de 1 AU. Con todo esto obtenemos que pdin =2nPa. Así, se obtiene un valor de F=pdinL2=2,45 N (2.6) Aunque, el empuje que proporcionan las E-sail pueda parecer pequeño, debido a que el peso de los tethers por unidad de longitud es muy pequeño, se pueden conseguir velocidades del orden de los 50 km/s según se indica en [ 6 ], valor que es mucho mayor que las velocidades que se alcanzan en un cohete de combustible sólido ( ∼ 1 km/s) o en un cohete bipropelente ( ∼ 5 km/s). Además, una de las ventajas que tienen las E-sail con respecto a las velas solares convencionales es que en el caso de las primeras el empuje es inversamente proporcional a la distancia que separa el Sol del vehículo espacial, mientras que en las segundas el empuje es inversamente proporcional al cuadrado de la la distacia que separa el Sol del vehículo espacial, como se indica en [ 3 ]. Es por tanto, que con las E-sail se pueden realizar misiones más lejanas. Esta diferencia, en la relación entre el empuje y la distancia al Sol en velas eléctricas y convencionales, se debe a que en ambos casos al alejarnos del Sol se produce un aumento de la longitud de Debye, que se trata de la escala a través de la cual los electrones móviles generan un apantallamiento del campo eléctrico en los conductores, de tal forma que el potencial eléctrico disminuye y por tanto el empuje. Sin embargo, en las E-sail esta disminución se ve contrarrestada porque a medida que aumenta la distancia entre el vehículo y el Sol, el área efectiva de la E-sail aumenta de tamaño también y por tanto, se produce un crecimiento del diámetro del campo eléctrico y como consecuencia de ello, del empuje como se indica en [12]. Se procede ahora a ver cuál es la potencia consumida para evacuar los electrones del plasma del viento solar fuera de los cables de la malla, para ello se tiene que Ie=enveAs=17,7 mA (2.7) siendo ela carga que posee un electrón, mientras que As= (2N)2πrL =15393,8 m2(2.8) es el área total de todos los cables conductores y ve=106 m/s es la velocidad de los electrones del viento solar. Tomando como valor de voltaje uno razonable, por ejemplo V=6kV, se tiene que P=VIe=106,2 W (2.9)
18 Capítulo 3. Modelo básico del problema Figura 3.1 Sistema de referencia heliocéntrico junto con el sistema de referencia orbital. Adaptada de [9]. • Sistema de referencia ejes cuerpo: Se trata de un sistema de referencia ligado al cuerpo principal del vehículo espacial y que se mueve solidariamente con él. Este sistema está definido por B ( S ; xB , yB , zB ), cuyos vectores unitarios son ( ˆ iB , ˆ jB , ˆ kB ), donde ˆ kB se encuentra alineado con el eje de simetría de del vehículo espacial, mientras que los vectores ( ˆ iB , ˆ jB ), definen el plano nominal de la vela, como se ve en la figura 3.2, obtenida de [ 2 ]. Además, S que es el origen del sistema de referencia, se encuentra en el centro de masas del vehículo espacial. Figura 3.2 Sistema de referencia ejes cuerpo. Extraída de [2]. La vela eléctrica consistirá en un número de cables (tethers), cada uno de ellos denotado por el subíndice k , con k∈[0,1,2,...,N−1] . La posición de cada tether en el plano ( xB , yB ) viene definida por el ángulo de azimut que comienza a medirse desde el eje xB en sentido antihorario. Este ángulo se define como:
3.1 Sistemas de referencia 19 ζk≜2π Nk(3.1) Por otro lado, la orientación del vehículo espacial viene definida por [ αn , δn ]. Donde αn∈ [0, π /2] es el ángulo de inclinación (pitch angle), siendo este el ángulo entre ˆr y el eje zB , definiéndose ˆr como el vector unitario que determina la dirección Sol-vehículo y que coincide aproximadamente con la dirección del viento solar. Mientras que δn∈ [0, 2 π ), es el ángulo horario, siendo este el formado por la proyección de ˆr sobre el plano nominal de la vela y el eje xB . Basándonos en las anteriores definiciones de ambos ángulos, las expresiones que los determinan serán las mostradas en (3.2) y (3.3). αn≜arccos(ˆr·ˆ kB)(3.2) δn≜ arccosˆr·ˆ iB ∥ˆr·ˆ kB∥si (ˆr·ˆ jB)≥0 2π−arccosˆr·ˆ iB ∥ˆr·ˆ kB∥si (ˆr·ˆ jB)<0 (3.3) Estas relaciones trigonométricas se pueden obtener fácimente de la figura 3.3, extraída de [2]. Figura 3.3 Representación gráfica de los ángulos de actitud. Extraída de [2]. • Sistema de referencia inercial ligado al cuerpo: Como su propio nombre indica, se trata de un sistema de referencia inercial, cuyo origen es el centro de masas del vehículo y que viene definido por I ( S ; xI , yI , zI ) siendo la base de vectores unitarios ( ˆ iI , ˆ jI , ˆ kI ). En la situación de partida, estos vectores unitarios coinciden con los vectores unitarios de los ejes cuerpo. Dicho sistema se usa para referir la actitud que tiene el vehículo en un determinado instante con respecto a la que tenía inicialmente. • Sistema de referencia auxiliar: Se trata de un sistema de referencia no inercial. Este sistema viene denotado por A ( S ; xA , yA , zA ), la base de vectores unitarios que lo define es ( ˆ iA , ˆ jA , ˆ kA ). El vector ˆ jA coincide en todo momento con el vector ˆ jI de tal forma que tanto xA como zA se encuentran girando alrededor del eje yA hasta alcanzar ambos un ángulo αn con respecto a los ejes xI y zI respectivamente, cuando la maniobra es completada. Este sistema de referencia, definido en [ 2 ], se usa posteriormente para definir la evolución temporal de αn es por ello que se incluye aquí. Un esquema visual de dicho sistema se puede ver en la figura 3.14.
20 Capítulo 3. Modelo básico del problema 3.2 Situación de partida A la hora de definir la posición desde la que partirá el vehículo cabe mencionar que el viento solar es variable, este pierde intensidad al alejarnos del Sol y por ello la fuerza que experimenta la E-sail variará a lo largo de su trayecto, tal y como se expuso en el capítulo 2. Es por ello, que seleccionar una situación de partida u otra influye en el empuje inicial que experimente el vehículo espacial. Dicho esto, parece razonable elegir como posición de partida aquella en la que el vehículo se encuentre a una distancia de 1 AU del Sol coincidiendo ˆr con el vector ˆ kI en dirección, en la que el empuje al inicio de la misión no será demasiado pequeño. Dicha configuración inicial se muestra en la figura 3.4, extraída de [2]. Figura 3.4 Situación de partida. Extraída de [2]. 3.3 Representación de la actitud La actitud de un vehículo se define como la orientación relativa de un sistema de referencia ligado a este respecto a un sistema de referencia concreto, en nuestro caso estos sistemas serían el sistema de referencia ejes cuerpo By el sistema de referencia inercial Irespectivamente. En general, hay varias maneras de representar la actitud de un vehículo como los ángulos de Euler, los cuaterniones o la matriz de los cosenos directores. Cada representación, presenta ventajas e inconvenientes. Uno de los beneficios de que existan varias formas de representar la actitud es que la dificultad de un problema puede ser rebajada en función de si se usa un tipo de representación en lugar de otro. De esta forma, se puede trabajar con el modo de representación de la actitud que más convenga y que menos inconvenientes genere a la hora de abordar un problema. Además, es importante destacar también que posible pasar de un tipo de representación a otra. Otro aspecto interesante a comentar es que existe la posibilidad de que dos combinaciones de parámetros representen la misma actitud. En este caso se dice que la representación tiene ambigüedades. Lo que se hace en la mayoría de estas situaciones es limitar los posibles valores que pueden tomar los parámetros que definen
3.3 Representación de la actitud 21 la actitud. Algunas de dichas ambigüedades se comentarán para cada caso de representación de la actitud cuando se definan cada una de ellas. Por otro lado, hay casos en los que es necesario pasar por sistemas de referencia intermedios para alcanzar una actitud final partiendo de una inicial. En estos casos lo que se hace para llegar a esa actitud final es componer actitudes, dicha composición dependerá de la representación utilizada. 3.3.1 Matriz de cosenos directores (DCM) Dado el sistema de referencia I cuya base de vectores unitarios es ( ˆ iI , ˆ jI , ˆ kI ) y el sistema de referencia B determinado por ( ˆ iB , ˆ jB , ˆ kB ), ambos definidos en la Sección 3.1, la orientación de B con respecto a I , está completamente determinada por la matriz de cambio de base CB I , que permite expresar un vector genérico v , definido inicialmente en la base de vectores unitarios de I , en la base de vectores unitarios de B mediante: vB=CB IvI(3.4) donde CB I= c11 c12 c13 c21 c22 c23 c31 c32 c33 (3.5) Teniendo en cuenta la expresión 3.4 el vector ˆ iIexpresados en ejes cuerpo sería: ˆ iB I=CB Iˆ iI I=CB I[1,0,0]T= [c11,c21,c31]T(3.6) Usando el resultado de 3.6, si se realiza el producto escalar entre ˆ iB y ˆ iI estando ambos expresados en la base de vectores unitarios de B, se llega a que ˆ iB·ˆ iI= [1,0,0]·[c11,c21,c31]T=c11 (3.7) Por lo que el elemento c11 de la matriz DCM sería el coseno director entre ˆ iB y ˆ iI , de ahí el nombre que esta toma. Con el resto de términos se procedería igual y entonces quedaría que: CB I= ˆ iB·ˆ iIˆ iB·ˆ jIˆ iB·ˆ kI ˆ jB·ˆ iIˆ jB·ˆ jIˆ jB·ˆ kI ˆ kB·ˆ iIˆ kB·ˆ jIˆ kB·ˆ kI (3.8) Razonando de la misma forma se tendría que el cambio de base inverso vendría dado por la siguiente matriz: CI B= ˆ iB·ˆ iIˆ jB·ˆ iIˆ kB·ˆ iI ˆ iB·ˆ jIˆ jB·ˆ jIˆ kB·ˆ jI ˆ iB·ˆ kIˆ jB·ˆ kIˆ kB·ˆ kI (3.9) puesto que, CB I= (CI B)−1 , al observar las matrices (3.8) y (3.9) vemos que CI B= (CB I)T , por lo que se concluye que la DCM que define la actitud a través de nueve parámetros dependientes entre sí no puede ser cualquiera, sino que debe ser una matriz ortogonal y con determinante +1. Esta última condición implica que el sistema de referencia en cuestión será dextrógiro, es decir, girará en el sentido de las agujas del reloj. Por otro lado, la composición de actitudes en este tipo de representación viene dada por un simple producto matricial, es decir, supongamos unos sistemas de referencia genéricos S1 , S2 y S3 , y supongamos además que se tiene que la actitud de S2 respecto a S1 viene dada por CS2 S1 y que la actitud de S3 respecto a S2 viene definida por CS3 S2, la actitud de S3respecto a S1vendrá determinada por CS3 S1=CS3 S2CS2 S1. 3.3.2 Ángulos de Euler Estos ángulos permiten definir la actitud de un sistema de referencia respecto al sistema de referencia origen mediante una secuencia de tres rotaciones en ejes no consecutivos. Existen hasta 12 posibles secuencias de ángulos de Euler para representar la actitud. Esto implica que el orden en el que se realizan las rotaciones es importante, pues es posible llegar a una misma actitud pero usando diferentes ternas de rotaciones. La elección de una u otra dependerá de la aplicación. En particular, en este proyecto se seguirá la secuencia:
22 Capítulo 3. Modelo básico del problema Iψ −→ zI S′φ −→ xS′ S′′ θ −→ yS′′ B Dicha secuencia se puede explicar de la siguiente manera: se parte del sistema de referencia inercial I , primeramente se ejecuta un giro de ángulo ψ alrededor del eje zI , mediante dicho giro se pasa a unos ejes intermedios que se denotarán por S′ . Una vez se está en los ejes S′ , se realiza un giro de ángulo φ alrededor del eje xS′ y se llega a unos nuevos ejes S′′ que serán también intermedios. Partiendo de estos ejes S′′ se vuelve a realizar una rotación de ángulo θ pero esta vez alrededor del eje yS′′ llegando finalmente a los ejes cuerpo B . Cabe aclarar que los ejes S′ y S′′ no son más que sistemas auxiliares, no tienen ningún tipo de aplicación en el proyecto, es por ello que no se definen como sí que se ha hecho con los ejes cuerpo y los ejes inerciales I . La obtención de la matriz de cosenos directores a partir de los ángulos de Euler se realiza mediante la multiplicación de matrices de rotación elementales que se muestran a continuación: CB S′′(φ) = 1 0 0 0 cosφsinφ 0−sinφcosφ CS′′ S′(θ) = cosθ0−sinθ 0 1 0 sinθ0 cosθ CS′ I(ψ) = cosψsinψ0 −sinψcosψ0 0 0 1 (3.10) Para el caso que nos ocupa, que es el de la secuencia mostrada anteriormente, la DCM se obtendría mediante el producto matricial representado en la expresión (3.11), este producto entre las matrices de rotación correspondiente a cada ángulo de Euler debe seguir el orden inverso al realizado en la secuencia, al igual que se indicó en la Sección 3.3.1: CB I=CB S′′(θ)CS′′ S′(φ)CS′ I(ψ)(3.11) Una vez realizado este producto matricial, la matriz a la que se llega es la mostrada en (3.12) CB I= cosθcosψ−sinθsinψsinφcosθsinψ+sinθsinφcosψ−sinθcosφ −cosφsinψcosψcosφsinφ sinθcosψ+sinφcosθsinψsinθsinψ−sinφcosθcosψcosθcosφ (3.12) Obsérvese que analizando la matriz (3.12), se puede comprobar que la actitud mostrada por ( 180 +ψ , 180 −φ , 180 +θ ) representa la misma orientación que ( ψ , φ , θ ). Es por ello, que generalmente se suele limitar φ∈[-90º, 90º]. La gran ventaja de usar los ángulos de Euler es el significado físico que aportan. Sin embargo, operar con dichos ángulos a veces es complicado, pues en muchos casos se tienen singularidades que evitan que en ciertas situaciones la actitud esté bien definida. Además, otro problema que surge al usar este tipo de representación, es que a la hora de componer actitudes, cuando se tienen rotaciones mixtas, las matrices resultantes son demasiado complejas para ciertas aplicaciones. 3.3.3 Ángulo y eje de Euler Según el Teorema de Euler: “El movimiento más general posible de un sólido con un punto fijo es una rotación alrededor de un único eje”. A partir del Teorema anterior, se puede ver que la actitud de un sistema respecto de otro, en nuestro caso B y I respectivamente, queda totalmente definida indicando el ángulo de giro θ y el eje respecto del cual se produce ese giro, eB/I=[ex,ey,ez]T. Esto se traduce en que el eje de Euler es el eje en torno al cual debe de producirse el giro del sistema de referencia origen para obtener la actitud del sistema de referencia objeto de representación. Mientras que el ángulo de Euler es la cantidad a girar en torno al eje de Euler. Es por ello, que tiene la ventaja de poseer un sencillo significado físico. Sin embargo, presenta ambigüedades para ciertas rotaciones, por ejemplo un giro de 0ºy uno de 360ºson físicamente iguales, pero para el primero θ=0 , mientras que para el segundo no está bien definido. Esencialmente, es una representación que se usa para análisis teóricos o para ángulos pequeños, es por ello que no se entra en ella con mucho más detalle. Aunque otro posible uso sería el de obtener con cierta facilidad errores de apuntamiento, es decir, este tipo de representación permite ver con un solo número (ángulo de Euler) cuál es la desviación existente entre dos sistemas de referencia.
3.3 Representación de la actitud 23 3.3.4 Cuaterniones Los cuaterniones fueron introducidos por Hamilton en el siglo XIX como herramienta matemática necesaria para la resolución de ciertos problemas físicos. Los cuaterniones son una extensión de los números complejos a cuatro dimensiones. Un cuaternión genérico se trata de una cuaterna de número q0 , q1 , q2 , q3 , que se escribe de la siguiente manera: q=q0+iq1+jq2+kq3(3.13) Generalmente q0 se denomina parte escalar, mientras que q= [ q1 , q2 , q3 ] T es conocida como la parte vectorial de q. Algunos aspectos importantes de su álgebra se muestran a continuación: • Suma: La suma se realiza componente a componente, es decir, dados dos cuaterniones q=q0+iq1+ jq2+kq3yq′=q′ 0+iq′ 1+jq′ 2+kq′ 3, se tiene que q+q′= (q0+q′ 0)+i(q1+q′ 1)+ j(q2+q′ 2)+k(q3+q′ 3)(3.14) • Producto: El producto se da componente a componente, conociendo las siguientes reglas de multiplicación: i·i=−1,i·j=k,i·k=−j,j·i=−k,j·j=−1,j·k=i,k·i=−i,k·k=−1(3.15) Se tiene tamién la fórmula de Hamilton que implica que i·j·k = − 1. Además, como se puede deducir de las relaciones mostradas en (3.15), el producto no es conmutativo, es decir, q·q′=q′·q. El producto entre dos cuaterniones también se puede escribir en forma matricial y en forma vectorial [10]. • Conjugado: Se define el conjugado de un cuaternión de manera similar a cómo se hace en el caso de los números complejos, es decir, si se tiene que q=q0+iq1+jq2+kq3su conjugado será: q∗=q0−iq1−jq2−kq3(3.16) •Módulo: El módulo se define como: |q|2=q·q∗=q2 0+q2 1+q2 2+q2 3(3.17) •División: Se define haciendo uso del conjugado y del módulo: q′/q=q′/q·q∗/q∗= (q′·q∗)/|q|2(3.18) Una vez definida el álgebra del cuaternión, cabe decir que se trata de una representación matemática sin sentido físico, pero que ofrece ciertas facilidades a la hora de operar con ellos. Dada una actitud determinada mediante el ángulo y el eje de Euler ( e y θ ), el cuaternión que representa a esta actitud se puede definir como: q0=cosθ/2, q=sinθ/2e(3.19) Además, se observa de las expresiones (3.19) que |q|=1 , propiedad que se debe cumplir para un cuaternión que represente una actitud. Sin embargo, se puede observar que para cuaterniones también existe un problema de ambigüedad, pues la actitud representada por ( e , θ ) es equivalente a la representada por ( −e , 360 −θ ), lo cual indicaría que q=−qestableciéndose que un cuaternión y su opuesto representan la misma actitud. Por otro lado, la fórmula que se sigue para pasar de cuaterniones a DCM, C , es la de Euler-Rodrigues para cuaterniones: C= (q2 0−qTq)Id +2qqT−2q0q×(3.20)
24 Capítulo 3. Modelo básico del problema donde el operadorq×toma la siguiente forma: q×= 0−q3q2 q30−q1 −q2q1 0 (3.21) Aplicando la fórmula de Euler-Rodrigues, (3.20), se llega a: C(q) = q2 0+q2 1−q2 2−q2 32(q1q2+q0q3)2(q1q3−q0q2) 2(q1q2−q0q3)q2 0−q2 1+q2 2−q2 32(q2q3+q0q1) 2(q1q3+q0q2)2(q2q3−q0q1)q2 0−q2 1−q2 2+q2 3 (3.22) Esta matriz es la DCM en función del cuaternión asociado, que permitirá transformar los datos proporcionados en ejes inerciales a ejes cuerpo. Si comparamos término a término la matriz (3.22) con la matriz de (3.12) se pueden obtener las siguientes relaciones que son las que permiten conocer los tres ángulos de Euler a partir del cuaternión asociado. φ=arcsin(2(q2q3+q0q1)),θ=arctan2(q0q2−q1q3) q2 0−q2 1−q2 2+q2 3,ψ=arctan2(q0q3−q2q1) q2 0−q2 1+q2 2−q2 3 (3.23) donde el problema de soluciones múltiples estaría solventado al limitar el rango de variación de φ∈ [-90º, 90º], por la razón mencionada anteriormente. Una vez hecho esto ya se podría conocer en que cuadrante se encuentran θ y ψ mediante los términos, (1, 3) y (3, 3), y, (2, 1) y (2, 2), de la matriz (3.12) respectivamente. Sin embargo, las fórmulas que permiten pasar de la DCM C al cuaternión asociado son las mostradas en (3.24) q0=p1+Tr(C) 2, q×=1 4q0 (CT−C)(3.24) siendo Tr(C)la traza de la matriz DCM. 3.4 Cinemática de la actitud En esta sección, se entra a describir las ecuaciones cinemáticas de la actitud para cada una de las representaciones descritas anteriormente. En el caso de la actitud, las ecuaciones cinemáticas diferenciales (EDC) son las que relacionan la representación de la actitud con la velocidad angular. Generalmente, estas ecuaciones son no lineales. Es interesante conocer las diferentes EDC para las diferentes representaciones de la actitud, para así ver cuál es la que tiene más ventajas desde el punto de vista computacional. 3.4.1 Ecuaciones diferenciales cinemáticas para la DCM A continuación, se muestra la ecuación cinemática diferencial correspondiente a la representación mediante DCM a efectos descriptivos únicamente, ya que, no se entrará en mucho detalle, puesto que, estas ecuaciones no se usarán a lo largo del desarrollo del proyecto. Para ello, hay que tener en cuenta que el sistema de referencia B gira respecto a I con una velocidad angular ΩB B/I . Siguiendo el desarrollo que aparece en [ 10 ], se llega a la siguiente expresión: d dt CB I=−( ΩB B/I)×CB I(3.25) se trata de una ecuación diferencial matricial, que habrá que resolver componente a componente, por lo que se tendrían nueve ecuaciones diferenciales acopladas dado que dichos componentes son dependientes entre sí. Sin embargo, no es muy operativo trabajar con ella, ya que la integración numérica de las ecuaciones formuladas en términos de la DCM no garantiza que la matriz resultante sea ortogonal debido a errores de redondeo que se introducen y la corrección de la matriz mediante un método iterativo para que esta sea ortogonal hace que este tipo de representación no sea eficiente. En cambio, es bastante útil para realizar el cambio de base de un vector de un sistema de referencia a otro.
3.4 Cinemática de la actitud 25 3.4.2 Ecuaciones diferenciales cinemáticas para los ángulos de Euler Partiendo de la definición de la secuencia mostrada anteriormente, se tiene que la relación entre las velocidades angulares obtenidas de la mecánica y correspondientes a las rotaciones, expuestas en la Sección 3.3.2, de cada sistema de referencia es: ΩB B/I= ΩB B/S′′ + Ωb S′′/S′+ ΩB S′/I(3.26) donde ΩB B/I es la velocidad angular con la que B (sistema de referencia ejes cuerpo) gira respecto a I (sistema de referencia inercial). La definición del resto de velocidades angulares se realizaría de forma análoga. Como consecuencia de las rotaciones producidas durante el proceso de composición de actitud de ejes cuerpo con respecto a ejes inerciales, se tiene que ΩB B/S′′ = [ ˙ φ , 0, 0] T , ΩS′′ S′′/S′= [0, ˙ θ , 0] T y ΩS′ S′/I= [0, 0, ˙ ψ ] T , por lo que para llegar a la ecuación 3.26 tendremos que expresar cada velocidad angular en el sistema de referencia Bmediante las matrices de rotación elementales mostradas en (3.10). ΩB B/I= ΩB B/S′′ +CB S′′ ΩB S′′/S′+CB S′CS′′ S′ Ωb S′/I(3.27) Mediante la ecuación (3.27), se llega al sistema de ecuaciones mostrado por la expresión (3.28) ΩB B/I= cosθ0−cosφsinθ 0 1 sinφ sinθ0 cosφcosθ ˙ φ ˙ θ ˙ ψ (3.28) en la que, para obtener una expresión para las derivadas de los ángulos de Euler en función de ΩB B/I = [ Ω1 , Ω2,Ω3], habría que invertir la matriz llegándose al sistema (3.29). ˙ φ ˙ θ ˙ ψ =1 cosφ cosφcosθ0 cosφsinθ sinφsinθcosφ−cosθsinφ −sinθ0 cosθ Ω1 Ω2 Ω3 (3.29) Como se puede ver, el sistema de ecuaciones (3.29) está compuesto por 3 ecuaciones diferenciales no lineales con multitud de funciones trigonométricas que además posee una singularidad para φ=±90 º. Lo mismo sucede para cualquier otro conjunto de ángulos de Euler, posee singularidades. A pesar de ello, en el caso que aquí se trata estos problemas de singularidades no se tienen, ya que, los valores que toman los ángulos de Euler son pequeños. 3.4.3 Ecuaciones diferenciales cinemáticas para cuaterniones Partiendo de la formulación de las ecuaciones cinemáticas diferenciales en función del ángulo y del eje de Euler, que se muestran en las expresiones de (3.30) ˙ eB B/I=1 2(eB B/I)×+1 tanθ/2(Id −eB B/I(eB B/I)T) ΩB B/I,˙ θ= (eB B/I)T ΩB B/I(3.30) para obtener las ecuaciones diferenciales cinemáticas formuladas con cuaterniones es necesario derivar las expresiones (3.19) respecto del tiempo y sustituir ˙ eB B/Iy˙ θen ellas. Tras esto se llega a que: ˙q= 1 2qT ΩB B/I 1 2[q×+q0Id] ΩB B/I!=1 2 −q1−q2−q3 q0−q3q2 q3q0−q1 −q2q1q0 Ωx Ωy Ωz (3.31) Se tratan de cuatro ecuaciones diferenciales bilineales y sin singularidades. Estas propiedades son la razón más importante por la cual generalmente son las usadas a la hora de representar la actitud en vehículos espaciales. Se realizan los cálculos necesarios con cuaterniones y posteriormente, si es necesario visualizar la solución, se transforma esta a ángulos de Euler. Debido a la acumulación de errores a causa de imprecisiones numéricas, que provoca que el módulo del cuaternión que representa la actitud ya no sea unitario, es posible llevar a cabo un proceso de renormalización. En cualquier caso, la manera de corregir esta fuente de error es mucho más sencilla que en el caso relativo a los problemas de ortogonalidad en la DCM.
26 Capítulo 3. Modelo básico del problema 3.5 Definición de las diferentes fuerzas que intervienen en el problema En este apartado, se presentan las distintas fuerzas que intervienen en el problema, aunque no todas ellas se tendrán en cuenta, pues la contribución de algunas, como puede ser el caso de las fuerzas elásticas o de las fuerzas de amortiguación se desprecian en el modelo implementado en el presente trabajo, por lo que simplemente se comenta la posibilidad de la existencia de las mismas. Dicho esto, las fuerzas que actúan son las siguientes: • Fuerza gravitatoria: Esta fuerza aparece debido a la atracción gravitatoria que el Sol ejerce sobre el vehículo por el hecho de que este tiene masa. El modelo gravitatorio que se usará a la hora de plantear las ecuaciones dinámicas viene definido por el Problema de los dos cuerpos, y será el mostrado por la ecuación (3.32) Fg=−G·M·m r2·ur(3.32) dónde G representa la constante de gravitación universal, cuyo valor es G=6,67 ·10−11 N m 2 /kg 2 . Por otro lado, M y m son las masas de los cuerpos que interaccionan, en nuestro caso M se correspode con la masa del Sol mietras que m con la masa del vehículo espacial. Por último, r hace referencia a la distancia que separan a ambas masas y ur es el vector unitario que expresa la dirección de actuación de la fuerza. • Fuerzas inerciales: Su existencia se debe al movimiento del vehículo espacial a lo largo de su órbita. En concreto, encontramos una fuerza inercial que se opondrá al movimiento orbital del vehículo en torno al Sol, y otra fuerza inercial que se opondrá al de giro de la E-sail alrededor de su eje, definido por el eje zB. Ambos movimientos están acoplados, a causa del efecto Coriolis, como se puede ver en [8]. • Fuerza de Coulomb: Su inclusión se debe a la interacción que se produce entre los protones que se encuentran confinados en los tethers y las partículas ionizadas del viento solar. Una estimación de la misma ya ha sido realizada en la Sección 2.2. Es esta fuerza, la que permite la generación de empuje y la que nos permitirá controlar la actitud de la E-sail. Su magnitud dependerá del viento solar siendo la velocidad y energía de este variables, pues estas últimas varían en función de la distancia a la que el vehículo se encuentre del Sol. • Fuerza de ligadura: Esta fuerza surje al considerar la composición total del vehículo como un sistema multicuerpo, formado por los cables y el cuerpo principal del vehículo, que se idealiza como una masa puntual. Es por ello, que será necesario eliminar el vínculo entre ambos para poder ver la contribución que el cuerpo principal del vehículo ejerce sobre el tether y viceversa. • Fuerzas elásticas: Estas fuerzas aparecen debido a que se considera a los tethers como elementos que se pueden deformar y que tienen rigidez, por lo que dependerán del modelo que se elija a la hora de representar la deformada de los tethers. • Fuerzas de amortiguamiento: Son fuerzas que surgen debido a la consideración del tether como un elemento que posee amortiguación. Sucede igual que en el caso de las fuerzas elásticas, dependerán del modelo de deformada que se elija. 3.6 Modelo de la deformada del tether Para el presente proyecto, se considera que los tethers son elementos que se pueden deformar debido a la interacción de estos con el viento solar. A la hora de seleccionar el modelo que define la deformada de estos, se ha optado por aquel que usa la función logarítmica proporcionada por [2] y ya introducida en la Sección 2.4. La ventaja de emplear este modelo es que aporta una gran simplicidad sin dejar de representar fielmente la realidad, siempre y cuando la vela eléctrica tenga una configuración de “Sun-facing” en la que la recta Sol-vehículo espacial es perpendicular al plano nominal de la vela. En caso de no ser así el modelo pierde validez, pues la deformada no se aproximaría al arco descrito al graficar una función logarítmica.. Además, nos permite tener una expresión analítica con la que poder trabajar. Esta expresión viene dada por la ecuación (3.33).
3.7 Ecuaciones dinámicas de Euler 27 f=blLln(1+x L)(3.33) Donde f representa la distancia de un elemento diferencial de longitud del tether ( ds ) al plano ( xB , yB ), el parámetro x representa la distancia de ds al eje zB , mientras que L es la longitud del cable. Por otro lado, la variable bl que aparecen en (3.33) se trata de un coeficiente adimensional de forma cuya expresión viene dada por la ecuación (3.34) bl≜2σvvs ρ ω2L(3.34) siendo vvs la velocidad del viento solar cuyo valor ya es conocido de secciones anteriores, ρ es la densidad de masa lineal del tether ( ∽ 10 −5 kg/m), mientras que ω es la velocidad de giro de la E-sail. Adicionalmente, σ se trata de un parámetro de diseño que representa a la fuerza de Coulomb y que puede modificar su valor en función de la tensión eléctrica que tengan los tethers, sus unidades serían [kg/(m · s)]. Su expresión analítica viene dada por la ecuación (3.35) y esta tendrá validez siempre y cuando la distancia que separa al Sol y al vehículo espacial sea de 1 AU, lo cual supone una nueva limitación del modelo. σ=0,18 ·max(0,V−Vw)·pε0mpn(3.35) En la ecuación (3.35), V representa el voltaje del tether, Vw es el voltaje de las partículas ionizadas del viento solar (cuyo valor aproximado es de 1 kV, como se indica en [ 2 ]), ε0 es la permitividad del vacío siendo su valor ε0=8,8541878176 ·10−12 F/m, mientras que el resto de parámetros son ya conocidos de secciones anteriores. A continuación, en la figura 3.5 se muestra un esquema del modelo expuesto anteriormente con el objetivo de aclarar ideas. Dicha imagen ha sido extraída de [2]. Figura 3.5 Esquema del modelo de deformada del tether. Extraída de [2]. 3.7 Ecuaciones dinámicas de Euler En el presente apartado se presentan las ecuaciones de Euler para un sólido rígido, procedentes de la mecánica, que definen la evolución temporal de la velocidad angular bajo la acción de unos pares y que son integradas numéricamente junto con las ecuaciones diferenciales cinemáticas presentadas en la Sección 3.4.2 para obtener la actitud de la vela una vez transcurra un periodo de tiempo concreto. Para ello, nos basaremos en la formulación de la referencia [2].
34 Capítulo 3. Modelo básico del problema 3.8.2 Evoluciones temporales de las componentes del par de fuerzas total A continuación, se muestran cual es la evolución temporal de las distintas componentes ([ E , F , G ], definidas por las expresiones 3.45, 3.46, 3.47) que generan el momento de fuerzas total que actúa sobre la vela para el caso en el que se considera la maniobra y los parámetros de la E-sail definidos al inicio de la sección. Cabe recordar, que estas componentes incluyen tanto la contribución del momentos de fuerzas dado por la modulación del voltaje eléctrico, así como los efectos de la deformación en los tethers debido a la presión dinámica del viento solar y el par causado por la fuerza centrífuga debida al giro de la vela. Figura 3.11 Evolución temporal de la componente E del momento de fuerzas total considerando la maniobra definida y los valores de los parámetros que aparecen en la tabla 3.1. Figura 3.12 Evolución temporal de la componente F del momento de fuerzas total considerando la maniobra definida y los valores de los parámetros que aparecen en la tabla 3.1.
3.8 Maniobra de control 35 Figura 3.13 Evolución temporal de la componente G del momento de fuerzas total considerando la maniobra definida y los valores de los parámetros que aparecen en la tabla 3.1. 3.8.3 Ley de control activo Para conseguir que el vehículo espacial se desplace por el Sistema Solar es necesario la existencia de una componente transversal del empuje que será la que permita direccionar el vehículo. Para ello, el plano nominal de la vela eléctrica debe ser orientado de tal forma que se tenga un valor de αn pequeño, pero distinto de 0, con el fin de que la E-sail siga manteniendo simetría axial. Con esta sección, se pretende mostrar cómo la actitud del vehículo y, por tanto, αn , pueden ser ajustados con el par Tc presentado anteriormente. Para esto, se usarán los razonamientos deducidos en [2]. Como se dijo en la Sección 3.8.1, la presencia del plano de control divide el plano nominal de la E-sail en dos. Ambos semiplanos están sometidos a diferentes niveles de tensión. Dichos niveles de tensión vienen impuestos por la ley de control, de tal modo que el semiplano en el que se encuentra la proyección de ˆr está sometido a un nivel de tensión mayor que el semiplano restante. Con esta modulación lo que se pretende es contrarrestar el par perturbador y generar la componente transversal del empuje requerida. Antes de entrar a definir la ley de control, se establece que β=β ( t , αn ) ∈ [0, 2 π ] es el ángulo que define la orientación del plano de control. Dicho ángulo es el formado por Tc y el eje xA medido en sentido antihorario desde xA, como se muestra en la figura 3.14 adaptada de [2]. Este ángulo está definido por la ecuación (3.56), y tiende a π 2cuando Tatiende a 0. β=(arctanTd Tasi t ≤tf π 2si t >tf (3.56) El valor de σk correspondiente a un tether genérico dependerá del desplazamiento angular γk de este con respecto al plano de control, dicho desplazamiento angular es: γk≜ωt+ζk−β(3.57) Por tanto, la ley de control resultante sería σk=(σ−λTc 2si γk∈[0,π) σ+λTc 2si γk∈[π,2π)(3.58) donde λes el △σrequerido por unidad de Tdgenerado:
36 Capítulo 3. Modelo básico del problema Figura 3.14 Esquema de la maniobra de actitud con el par de control. Extríada de [2]. λ≜△σ Td =2 sin(π/N) vvs L2cosαn (3.59) Como se puede ver en (3.58), se trata de una ley de control en bucle abierto, es decir, sin realimentación, puesto que la salida de la ley de control no es remitida de nuevo a la entrada como primera información, permitiendo que esta se regule de forma independiente. Por otro lado, el incremento del parámetro de diseño en el semiplano en el cual se encuentra la proyección de ˆr viene dado por el término λTc 2 , siendo esta misma cantidad la que disminuye este mismo parámetro en el semiplano restante. Además, el valor de σ para todos los tethers que se encuentran en un mismo semiplano es el mismo, variando este únicamente con el tiempo y con αn. 3.8.4 Simulaciones numéricas Con esta sección, se muestra la evolución de αn , obtenida al integrar numéricamente, con el comando ODE45 de MATLAB, el sistema de ecuaciones diferenciales compuesto por las ecuaciones diferenciales cinemáticas (3.29) y las ecuaciones dinámicas de Euler (3.36), (3.37) y (3.38), que a continuación se presenta, de esta forma se reproducen y analizan los resultados del artículo que se ha usado como referencia ([2]). f1≡˙ φ=Ωxcosθ+Ωzsinθ f2≡˙ θ=Ωy−(Ωzcosθ−Ωxsinθ)tanφ f3≡˙ ψ= (Ωzcosθ−Ωxsinθ)secφ f4≡˙ Ωx=It−Iz ItΩyΩz+E It f5≡˙ Ωy=Iz−It ItΩxΩz+F It f6≡˙ Ωz=G Iz (3.60)
3.8 Maniobra de control 37 Para ello se toman los siguientes valores de los parámetros del sistema, siendo estos obtenidos de [ 2 ]: Iz=3It= 3000 kg m 2 , L= 2 km, ω= 0,0758 rad/s, N= 16, V= 16,5 kV, y Vw= 1 kV. Las condiciones iniciales son φ(0) = θ(0) = ψ(0) = Ωx=Ωy=0 y Ωz=ω . Por otro lado, los valores elegidos para ¯ αn f y t f son: ¯ αn f =5◦ytf=120 s. Además, esta evolución es comparada con la evolución que seguiría αn al resolver el mismo sistema de ecuaciones, pero eliminando las perturbaciones existentes e imponiendo que el único momento de fuerzas aplicado fuese el obtenido al exigir que la velocidad del sistema en las ecuaciones dinámicas de Euler sea ΩA=˙ ¯ αnˆ jA , donde ˙ ¯ αn viene dada por la ecuación (3.49). Se considera que la evolución que aquí se obtiene es la ideal, es decir, la que seguiría αn sin ningún tipo de perturbación. El proceso seguido para obtener el par anterior se muestra a continuación: It˙ Ωx+(Iz−It)ΩyΩz=Mx(3.61) It˙ Ωy+(It−Iz)ΩxΩz=My(3.62) Iz˙ Ωz+(It−It)ΩxΩy=Mz(3.63) Teniendo en cuenta que tanto Ωx como Ωy son nulas y que Ωy=˙ ¯ αn , la única componente del momento de fuerzas distinta de 0 sería My, cuya expresión es: My=It 6¯ αn f tf"1 tf−2t t2 f#[1(t)−1(t−tf)] (3.64) Por lo que las ecuaciones dinámicas de Euler a integrar junto con 3.29 quedarían: ˙ Ωx=It−Iz It ΩyΩz(3.65) ˙ Ωy=Iz−It It ΩxΩz+6¯ αn f tf"1 tf−2t t2 f#[1(t)−1(t−tf)] (3.66) ˙ Ωz=0(3.67) Figura 3.15 Evolución ideal (azul) de αnjunto con la evolución considerando perturbaciones (rojo) de αn.
38 Capítulo 3. Modelo básico del problema Figura 3.16 Evolución temporal de φ,θyψpara la maniobra correspondiente a αn f = 5◦. En la figura 3.15, se puede ver la evolución, considerando perturbaciones, de αn en comparación con la evolución ideal, refiriéndonos por evolución con perturbaciones a la evolución obtenida al resolver íntegramente el sistema de ecuaciones compuesto por (3.29), (3.36), (3.37) y (3.38). En ambos casos, se ha tomado un tiempo de simulación mayor al de 2 minutos para verificar que el comportamiento de la vela eléctrica es capaz de mantenerse en la actitud final propuesta mediante la ley de control que se indicó en la Sección 3.8.3. Como puede observarse, las dos evoluciones son prácticamente coincidentes desde que se inicia la maniobra de actitud hasta que se finaliza. Sin embargo, mientras que en la evolución ideal αn es capaz de mantenerse perfectamente en la actitud final deseada transcurrido un tiempo de 2 minutos, en la evolución normal αn tiene un comportamiento oscilatorio alrededor de un valor muy próximo al ángulo deseado que es de 5◦. Se muestran además las evoluciones temporales de las distintas variables de estado del problema en las figuras 3.16 y 3.17, y de las potencias de control en la figura 3.18 Como se puede ver la evolución de ψ es periódica en todo instante de tiempo, esto es debido a que la vela está constantemente girando alrededor de su eje de giro, siendo este el eje zB . Por otro lado, las evoluciones de φ y θ tras un periodo de tiempo, que coincide con la fase en la que se realiza la maniobra para llegar a la actitud deseada, alcanzan también una periodicidad en sus evoluciones. Es interesante resaltar que debido a la presencia del par perturbador y a la componente G del par total en el eje zB , se produce una pequeña variación de Ωz que hace que αn no se mantenga oscilando en ese valor próximo al deseado, sino que su valor decae al aumentar el tiempo de integración, como se indica en [ 2 ] y como se puede ver en la figura 3.19. Por lo que se concluye que la ley de control presentada en la referencia [ 2 ], si bien funciona razonablemente bien para intervalos de tiempo cortos, si queremos que el vehículo espacial se mantenga en una determinada actitud un periodo de tiempo mayor como puede ser una hora, esta ley de control no es muy efectiva, ya que, el comportamiento de αn sería mejor en el caso de que no hubiese ley de control aplicada, lo que se correspondería con un valor de σk = 0 y por tanto, con un valor nulo de G .
3.8 Maniobra de control 39 Figura 3.17 Evolución temporal de las velocidades angularespara la maniobra correspondiente a αn f = 5◦. Figura 3.18 Evolución temporal de las variables de control para la maniobra correspondiente a αn f = 5◦.
40 Capítulo 3. Modelo básico del problema Figura 3.19 Evolución de αncon Gnulo (rojo) y con Gno nulo (azul)..
4 Aplicación del control LQR Co n este capítulo, se pretende mejorar la ley de control propuesta en la referencia [ 2 ]. Tal y como se muestra en la Sección 3.8.4 la ley de control propuesta en [ 2 ] presenta algunos aspectos mejorables, entre ellos: la eliminación de las oscilaciones que se producen cuando se llega al régimen permanente y la estabilización de la vela en la actitud final deseada pues el valor de αn decae con el tiempo, como se ve en la figura 3.19. Para ello, primero se lleva a cabo un proceso de adimensionalización de todas las variables del problema. Posteriormente, se lleva a cabo la linealización del sistema de ecuaciones diferenciales (3.60). Por último, se presenta el control LQR tanto de horizonte infinito como de horizonte finito que se usará para obtener las ganancias que definan la nueva ley de control. 4.1 Adimensionalización del problema Para evitar la proliferación de valores muy elevados que oculten posibles errores que se cometan durante el proceso, se llevará a cabo un proceso de adimensionalización. De esta forma aseguramos que las simulaciones numéricas realizadas son correctas pues los valores de las variables que intervienen no son elevados y de esta manera no se están ocultando errores que pueda haber a lo largo del proceso. Para ello, previamente se mostrarán las variables adimensionales a utilizar y posteriormente, se presentará como queda el sistema de ecuaciones (3.60) adimensionalizado. 4.1.1 Variables adimensionales Se presentan a continuación las variables adimensionales que se usarán a lo largo del desarrollo de este capítulo, y por tanto para definir la nueva ley de control. Las variables son las siguientes: ωx=Ωx ω,ωy=Ωy ω,ωz=Ωz ω(4.1) Γk=σk σ(4.2) t∗=ω·t(4.3) siendo ω la velocidad de giro nominal de la vela solar eléctrica , como se expuso ya con anterioridad, y estando σ definido por la ecuación (3.35), mientras que σk viene dada por las expresiones (3.58). Hay que tener en cuenta que los ángulos de Euler son ya adimensionales, por lo que no se introducen variables para ellos. 4.1.2 Sistema de ecuaciones adimensionalizado Teniendo en cuenta las variables adimensionales definidas en (4.1), (4.2) y (4.3), y aplicando la regla de la cadena se llega a que: dφ dt =dφ dt∗·dt∗ dt =dφ dt∗·ω(4.4) 41
42 Capítulo 4. Aplicación del control LQR dΩx dt =dΩx dωx·dωx dt∗·dt∗ dt =ω2·dωx dt∗(4.5) Procediendo de igual manera para el resto de variables de estado, el sistema definido en (3.60) en su forma adimensional quedaría de la siguiente manera: f= f1≡dφ dt∗·ω=ω·(ωxcosθ+ωzsinθ) f2≡dθ dt∗·ω=ωy·ω−ω·(ωzcosθ−ωxsinθ)tanφ f3≡dψ dt∗·ω=ω·(ωzcosθ−ωxsinθ)secφ f4≡dωx dt∗·ω2=It−Iz It·ω2ωyωz+E It f5≡dωy dt∗·ω2=Iz−It It·ω2ωxωz+F It f6≡dωz dt∗·ω2=G Iz (4.6) donde las componentes del par de fuerzas total, mostradas en las expresiones (3.45), (3.46) y (3.47); quedarian de la siguiente forma introduciendo la variable adimensional de control: Ek=uL2Γkσblsinαn[cos(δn−ζk)sinζk(ln8 −2)−sinδn(ln4 −1)]+ cosαnsinζk 2(4.7) Fk=uL2Γkσblsinαn[cosδn(ln4 −1)−cos(δn−ζk)cosζk(ln8 −2)]−cosαncosζk 2(4.8) Gk=uL2Γkσsin(δn−ζk)sinαn 2(4.9) 4.2 División en fases del problema. Linealización En el presente apartado se presentan las ecuaciones resultantes de linealizar el sistema de ecuaciones diferenciales adimensional (4.6). Este proceso de linealización se realiza con el objetivo de poder aplicar el método LQR que nos permite obtener la ley del control óptimo, y así mejorar los resultados mostrados en la referencia [2]. La linealización se hace en torno a la trayectoria de equilibrio en cada una de las fases de las que se compone el problema: una primera fase de tracking, donde el sistema dinámico se encuentra realizando la maniobra que le permite pasar de una actitud inicial a una final y una segunda fase de regulación donde las variables del problema llegan a régimen permanente y se mantienen en los valores que se desean tener al finalizar la maniobra de actitud. En el presente proyecto, primeramente se ha resuelto la fase de regulación y posteriormente la de tracking, con el objetivo de comprobar primeramente si la vela solar eléctrica era capaz de mantenerse en la actitud deseada durante una cantidad considerable de tiempo. 4.2.1 Valores de referencia Valores de referencia en la fase de regulación Los valores de referencia aquí mostrados, se obtienen del equilibrio del sistema, siendo este un equilibrio dinámico pues la vela al final de la maniobra sigue girando alrededor de su eje. Estos valores son los valores que tienen que tener las distintas variables del problema para que la actitud al final de la maniobra sea la deseada. Dichos valores, además, deben verificar las ecuaciones diferenciales que componen el sistema de ecuaciones (4.6). Es por ello, que el valor de referencia de algunas de las variables se obtiene a partir de dichas ecuaciones. Teniendo en cuenta la situación de partida definida en la Sección 3.2; la expresión (3.12) de la matriz DCM, que permite expresar un vector definido en ejes inerciales I en ejes cuerpo B y la definición de los parámetros αnyδnde la Sección 3.1, se llega a las expresiones:
4.2 División en fases del problema. Linealización 43 cosαn=cosθcosφ(4.10) sinδnsinαn=sinφ(4.11) cosδnsinαn=−sinθcosφ(4.12) Por lo que si fijamos que αre f n = 5 ◦ , que es la actitud que queremos mantener al final de la maniobra y tenemos en cuenta la expresión (4.10), podríamos fijar que φre f = 0 ◦ y, por tanto, obtendríamos que θre f =αre f n . El valor de referencia de δn se podría obtener fácilmente de la expresión (4.11) o de la (4.12), y este sería δre f n = 180 ◦ . Una vez hecho esto, considerando la expresión f2 del sistema (4.6) y teniendo en cuenta las imposiciones anteriores con respecto a φre f y a θre f se deduce que ωyre f = 0. Por otro lado, si imponemos que ωzre f = 1 de la ecuación f1 del sistema (4.6) obtendríamos que el valor de referencia de la velocidad angular en el eje xBadimensionalizada es ωxre f =−tanθre f . Por último, para obtener el valor de referencia de la variable de estado ψ se recurre a la tercera ecuación del sistema (4.6), en la que sustituyendo los valores de referencia de las variables anteriores se llega a que ˙ ψre f =ω/cosθre f , por lo que integrando se tendría que ψre f =ω·t/cosθre f . A modo de resumen, se muestran a continuación los valores de referencia de las distintas variables de estado adimensionalizadas Xre f = φre f θre f ψre f ωxre f ωyre f ωzre f = 0 αre f n ω·t/cosθre f −tanθre f 0 1 (4.13) Con el procedimiento seguido para el cálculo de los valores de referencia en la fase de regulación es posible asegurar que dichos valores son válidos para el proceso de linealización, pues satisfacen las ecuaciones del problema (4.6). Por otro lado, para obtener los valores de referencia de las variables de control, Γre f , se minimiza la siguiente función f( Γ) = Γ2 0+Γ2 1+Γ2 2+...+Γ2 N−1(4.14) pero teniendo en cuenta que, dichos valores tienen que satisfacer las ecuaciones dinámicas de Euler, siendo estas las expresiones f4 , f5 y f6 del sistema (4.6) ya adimensionalizado, lo que supondrán restricciones a la hora de resolver el sistema de ecuaciones resultante. Además, dichas ecuaciones deberán ser particularizadas para los valores de referencia ya calculados anteriormente, por lo que, los términos correspondientes a las derivadas de las velocidades angulares en las ecuaciones dinámicas de Euler serían nulos. Por tanto, las restricciones a las que estaría sujeto el sistema serían ω× re f ·I· ωre f =D· Γre f (4.15) donde I sería la matriz de inercias del sistema, mientras que, D sería la matriz de coeficientes de las variables de control, estando esta particularizada para los valores de referencia. Teniendo en cuenta esto el sistema (4.15) un poco más desarrollado sería 0−ωzre f ωyre f ωzre f 0−ωxre f −ωyre f ωxre f 0 · It0 0 0It0 0 0 Iz · ωxre f ωyre f ωzre f =D· Γre f (4.16) La matriz Dtendría la siguiente forma D= d11 ... d1N d21 ... d2N d31 ... d3N (4.17)
50 Capítulo 4. Aplicación del control LQR δ Γ=K(t)·δ X(4.52) desarrollándola un poco más se tiene que Γ= Γre f +K(t)( X− Xre f )(4.53) donde la ganancia K(t)se encuentra de la siguiente manera K(t) = −R−1(t)B(t)TP(t)(4.54) siendo P(t) la matriz resultante de resolver la Ecuación Diferencial de Ricatti que a continuación se muestra. Dicha ecuación es siempre resoluble pero no se puede resolver por los métodos convencionales del cálculo diferencial. Además, no se puede resolver en tiempo real, sino que hay que resolverla desde el futuro hacia el presente pues la condición de contorno está dada en el instante final (actitud final). −˙ P=AT(t)P(t)+P(t)A(t)−P(t)B(t)R−1(t)BT(t)P(t)+Q(t),P(tf) = Qfin (4.55) La ecuación (4.52) define la ley de control que se aplica en la fase transitoria, durante la realización de la maniobra. Al igual que antes en la fase de regulación dicha ley de control está compuesta por una referencia a seguir y una realimentación que pretende minimizar las desviaciones que se produzcan respecto a los valores de referencia. La principal diferencia ahora en relación con el régimen permanente es que se tiene una dependencia con el tiempo tanto de la ganancia como de los valores de referencia. 4.4 Resultados al aplicar el Control LQR 4.4.1 Resultados al aplicar el control LQR de horizonte infinito A continuación, se muestran los resultados obtenidos, con respecto a las variables del problema sin adimensionalizar, en las simulaciones numéricas al aplicar la ley de control (4.47) durante todo instante de tiempo. Figura 4.1 Evolución temporal de αnaplicando en todo instante de tiempo la ley de control (4.47). Como se puede ver en la figura 4.1 tras un transitorio la variable αn alcanza el valor deseado de 5 ◦ , que se fijo al inicio. Se ha logrado el primer objetivo que era el de eliminar las oscilaciones que se producían en
4.4 Resultados al aplicar el Control LQR 51 Figura 4.2 Evolución temporal de φ (azul) y θ (rojo) aplicando en todo instante de tiempo la ley de control (4.47). Figura 4.3 Evolución de temporal de ψaplicando en todo instante de tiempo la ley de control (4.47). dicha variable al alcanzar el régimen permanente. Además se puede ver como αn se mantiene en dicho valor, no decae con el tiempo, algo que sí ocurría al aplicar la ley de control de la referencia [2]. En la figura 4.2, se puede ver como también las variables φ y θ se estabilizan en torno a sus valores de referencia tras el transcurso de un cierto periodo de tiempo. Además, se puede observar como la evolución de θ es prácticamente la misma que la de αn , algo lógico pues se había fijado que durante la fase de regulación ambas variables tomasen el mismo valor de referencia. Por otro lado, en la figura 4.3 se ve como la evolución de ψ es periódica debido al valor de referencia de esta variable, que se tomó en la fase de regulación. Dicho valor era dependiente del tiempo y es esa dependencia la que produce la periodicidad que se observa. Como pasa con el resto de variables de estado, las velocidades angulares alcanzan los valores de referencia
52 Capítulo 4. Aplicación del control LQR Figura 4.4 Evolución temporal de las velocidades angulares aplicando en todo instante de tiempo la ley de control (4.47). fijados en la fase de regulación tras un periodo transitorio, siendo estos valores una vez se deshace la adimensionalización ωzre f =ω= 0,0758 rad/s (valor que se fijó en la Sección 3.8.4), ωyre f = 0 y ωxre f = −ω·tanθre f como bien se observa en la figura 4.4. Por ultimo, en la figura 4.5 se puede ver la evolución de todas las variables de control para el caso en el que se tienen 16 tethers. Como se puede observar el control ha quedado relativamente repartido con equidad entre los distintos tethers. Además, dichas variables también se estabilizan en sus valores de referencia cuando las variables de estado lo hacen, pues si se observa la ley de control (4.47) el término correspondiente a la diferencia entre los valores de las variables de estado y sus valores de referencia se anularía, y por tanto, las variables de control serían iguales a sus valores de referencia. Figura 4.5 Evolución temporal de las variables de control aplicando en todo instante de tiempo la ley de control (4.47).
4.4 Resultados al aplicar el Control LQR 53 4.4.2 Resultados al aplicar Control LQR de horizonte finito en la fase de tracking y Control LQR de horizonte infinito en la fase de regulación A continuación, se muestran los resultados obtenidos al aplicar, como se ha dicho al inicio de la sección, la ley de control (4.52) en la fase de tracking y posteriormente la ley de control (4.47) en la fase de regulación. Figura 4.6 Evolución temporal de αn al aplicar la ley de control (4.52) en la fase de tracking y posteriormente la ley de control (4.47) en la fase de regulación. Figura 4.7 Evolución temporal de φyθal aplicar la ley de control (4.52) en la fase de tracking y posteriormente la ley de control (4.47) en la fase de regulación. Como se puede ver en las figuras 4.6 y 4.7 la evolución en la fase transitoria de las variables aquí representadas ha mejorado con respecto a las que se veían en las figuras 4.1 y 4.2, puesto que ya no se producen tantas oscilaciones. Sin embargo, es cierto que en las figuras 4.6 y 4.7 se puede apreciar una cierta
54 Capítulo 4. Aplicación del control LQR oscilación en las variables αn y θ , justo antes de llegar al régimen permanente, coincidiendo esto con el instante en el que se pasa de aplicar una ley de control a aplicar la otra. A pesar de ello, dicha oscilación no es muy elevada y tanto en régimen transitoria (obviando la oscilación) como en régimen permanente las evoluciones se ajustan bastantes a las referencias (representadas en línea negra rayada). En esta sección, se estudia la mejoría de la sobreoscilación comentada a través de la consideración de un polinomio cúbico como valor de referencia en la variable αn durante la fase transitoria, pues al ser su derivada un polinomio de segundo grado la transición entre leyes de control podría ser más suave. Por otro lado, la evolución de ψ mostrada en la figura 4.8 es idéntica a la que aparecía en la figura 4.3, algo que tiene lógica pues los valores de referencia obtenidos en la fase de tracking y en la de regulación para dicha variable eran idénticos. Figura 4.8 Evolución temporal de ψ al aplicar la ley de control (4.52) en la fase de tracking y posteriormente la ley de control (4.47) en la fase de regulación. Con respecto a las velocidades angulares puede visualizarse en la imagen 4.9 que la evolución en la fase transitoria tanto de Ωx como de Ωy se ha suavizado bastante con respecto a las evoluciones mostradas en la figuras 4.4. Sin embargo, Ωz permanece prácticamente igual que en la fase de regulación, pues se han tomado los mismos valores de referencia en ambos casos. La componente Ωx se ajusta bastante bien a la evolución de su referencia, mientras que en la componente Ωy el ajuste ya no es tan exacto, hay diferencias entre los valores que toma la variable y los que toma su valor de referencia en la fase de tracking. Por otro lado, la evolución de Ωz difiere un poco de la evolución de su valor de referencia, esto se puede deber a las perturbaciones existentes, sin embargo, esta diferencia no es muy grande pues la escala que se muestra en la tercera imagen de la figura 4.9 es muy pequeña. Por último, en la figura 4.10 puede verse como las oscilaciones que se tenían en la fase transitoria al aplicar únicamente la ley de control (4.47) se han eliminado y ahora se tienen evoluciones mucho más suaves en dicha fase, aunque en ciertas variables de control, se puede ver como se producen pequeñas oscilaciones al darse el cambio en la aplicación de una ley de control a la otra. Además, se sigue teniendo el control relativamente repartido con equidad entre los distintos tethers que componen la vela.
4.5 Posibles extensiones y mejoras a realizar 55 Figura 4.9 Evolución temporal de las velocidades angulares al aplicar la ley de control (4.52) en la fase de tracking y posteriormente la ley de control (4.47) en la fase de regulación. Figura 4.10 Evolución temporal de las variables de control al aplicar la ley de control (4.52) en la fase de tracking y posteriormente la ley de control (4.47) en la fase de regulación.. 4.5 Posibles extensiones y mejoras a realizar El objetivo aquí es el de presentar diferentes alternativas con el fin de conseguir eliminar la oscilación que se produce al pasar de aplicar la ley de control (4.52) en la fase de tracking a la ley de control (4.47) en la fase de regulación. Primeramente, se ha modificado el valor de referencia que se ha tomado en la fase de tracking para la variable αn , con el objetivo de suavizar la transición entre una ley de control y la otra. Este nuevo valor propuesto es
56 Capítulo 4. Aplicación del control LQR αre f n=Ztf 0 K·t2·(tf−t)2(4.56) siendo K una constante a determinar para que se cumpla que la integral de K·t2·(tf−t)2 entre 0 y tf sea igual a ¯ αn f . Imponiendo esto el valor de Kobtenido es K=30 ¯ αn f t5 f (4.57) Por lo que el nuevo valor de referencia en el régimen transitorio de la variable αnes αre f n=30 ¯ αn f t5 f t3 3t2 f+t5 5−t4 2tf!(4.58) con las implicaciones que esto tendría en el resto de valores de referencia que depende de dicha variable. Con el valor de la expresión (4.58), se tendría que la aceleración ¨ αre f n sería nula tanto al inicio como al llegar al final de la maniobra, es por ello que la transición comentada anteriormente podría realizarse de una manera más suave de forma que se omita la oscilación que se observa en las figuras 4.6 y 4.7 con respecto a la evolución de las variables αnyθ. La comparativa se muestra en la figura 4.11 donde también se muestra los resultados obtenidos para el polinomio K·t3·(tf−t)3donde aquí Kse determina de la misma forma que se hizo anteriormente. Figura 4.11 Comparativa de las distintas evoluciones que se tienen para diferentes valores de referencia de αndurante la fase de tracking. Como se puede ver no hay mejora, la oscilación se sigue produciendo a pesar de modificar αre f n en la fase de tracking. Por lo que se considera una alternativa diferente. Al hacer una ampliación en la figura 4.10 en la zona donde se da la transición entre las leyes de control (en torno a 120 s), se observa que al cambiar la ley de control entre el transitorio y el régimen permanente se produce un escalón en las variables de control, que se puede observar en la figura 4.12. Por tanto, la solución adoptada para este caso, es realizar una pequeña transición que permita realizar el paso de una ley de control a la otra de una manera más suave, para así eliminar la oscilación que se genera en las evoluciones de αnyθ. Para ello, se introduce un filtro de primer orden, de tal forma que el valor de las variables de control en la transición es
4.5 Posibles extensiones y mejoras a realizar 57 Figura 4.12 Escalón producido en la evolución de las variables de control al pasar de la ley de control (4.52) en régimen transitorio a la ley de control (4.47) en régimen permanente. σtrans N=σtrack N(t)·e−c·t+σregul N·(1−e−c·t)(4.59) donde σtrack N(t) es el último valor de las variables de control durante la fase de tracking al aplicar la ley de control (4.52), mientras que σregul N se trata del primer valor que toman las variables de control tras la transición, una vez se alcanza el régimen permanente, donde se aplica la ley de control (4.47). Por otro lado, c es un parámetro a elegir que cuanto mayor sea, más rápido se producirá la transición. Es por ello, que se ha decidido tomar c = 1000. Por último, decir que el tiempo en la expresión (4.59) empezaría a contar cuando de comienzo la transición. Como se puede observar en la figura 4.13, al añadir esta pequeña transición al modelo la oscilación se ha eliminado, lo que supone una considerable mejora. Ya se consigue llegar a la actitud deseada o final en el tiempo de maniobra fijado. Además, se puede apreciar una pequeña diferencia entre la evolución de αn y su referencia, sin embargo, esta diferencia es pequeña, de tal forma que αn se ajusta bastante bien a los valores que toma su valor de referencia. Dicha oscilación en la evolución de la variable θ también se ha eliminado como se puede ver en la figura 4.14. Con respecto a φ , su evolución es muy similar a la mostrada en la figura 4.7 y se ajusta bastante bien a su referencia. Por las razones ya comentadas, la evolución temporal de ψ no sufre ninguna variación al introducir el filtro de primer orden que elimina las oscilaciones visualizadas en las variables αn y θ , como se puede ver en la figura 4.15. Por otro lado, con respecto a las velocidades angulares, en la figura 4.16 no se aprecia ninguna diferencia significativa con respecto a la figura 4.9 Para finalizar, en la figura 4.17 se puede ver que aparentemente al introducir el filtro de primer orden las evoluciones de las variables de control son iguales a las mostradas en la figura 4.10. Sin embargo, si nos centramos en las proximidades del instante de tiempo (120 s) en el cual se produce la transición de una ley de control a la otra (figura 4.18), se puede ver como el escalón observado en la figura 4.12 se ha eliminado, gracias a la transición introducida. Dicho escalón es el que originaba la oscilación en las evoluciones de las variables αn y θ , y su eliminación ha conllevado la omisión de las oscilaciones comentadas, por lo que el objetivo del proyecto se ha conseguido de manera satisfactoria.
58 Capítulo 4. Aplicación del control LQR Figura 4.13 Evolución temporal de αn al aplicar la ley de control (4.52) en la fase de tracking y posteriormente la ley de control (4.47) en la fase de regulación, añadiendo además el filtro de primer orden. Figura 4.14 Evolución temporal de θ y φ al aplicar la ley de control (4.52) en la fase de tracking y posteriormente la ley de control (4.47) en la fase de regulación, añadiendo además el filtro de primer orden.
4.5 Posibles extensiones y mejoras a realizar 59 Figura 4.15 Evolución temporal de ψ al aplicar la ley de control (4.52) en la fase de tracking y posteriormente la ley de control (4.47) en la fase de regulación, añadiendo además el filtro de primer orden. Figura 4.16 Evolución temporal de las velocidades angulares al aplicar la ley de control (4.52) en la fase de tracking y posteriormente la ley de control (4.47) en la fase de regulación, añadiendo además el filtro de primer orden.
66 Capítulo 5. Análisis paramétrico de los resultados obtenidos con la ley de control óptimo (a) 32 tethers. (b) 8 tethers. Figura 5.5 Evolución de las variables de control variando el número de tethers. En la figura 5.5 se puede ver como las evoluciones de las potencias de control en ambos casos son similares. La única salvedad es el número de variables de control que se tienen en un caso y en otro, estando dicho número directamente relacionado con el número de tethers que tiene la vela (una variable de control por cada theter que tiene la vela). Además, los valores que adquieren las variables de control en el caso en el que se
5.2 Análisis de la influencia de la actitud final deseada 67 tienen 32 tethers es menor pues al tener un mayor número de tethers el control está más repartido. Por otro lado, se sigue teniendo ese escalón al pasar de una ley de control a la otra, pues como se ha dicho ya este análisis se está realizando sin incluir el filtro de primer orden. 5.2 Análisis de la influencia de la actitud final deseada Aquí se estudia el comportamiento de las distintas variables del problema en función de la actitud final que se quiera tener que viene definida, como en ocasiones anteriores, por el parámetro αn. Al igual que antes, el resto de parámetros se consideran iguales a los que se tenían en la situación nominal, y sólo se varía el valor de ¯ αn f . En concreto se analizan dos casos, uno en el que ¯ αn f = 2,5 ◦ y otro en el que ¯ αn f = 10 ◦ Con respecto a αn , se puede visualizar en la figura 5.6 que no se aprecian cambios significativos en la forma de la evolución. Se sigue teniendo dicha oscilación, pues se está realizando el presente análisis para el caso en el que se tienen 16 tethers, si este número aumentase la amplitud de la oscilación se vería disminuida como se mostró en la Sección 5.1. En ambos casos αn consigue alcanzar el valor final deseado de una forma bastante exacta al igual que sucede en el caso en el que αn f = 5 ◦ . Sin embargo, dicho ajuste si es bastante preciso en la fase de tracking. Al igual que para el caso nominal, en los dos casos bajo estudio la evolución de θ es prácticamente igual a la de αn , mientras que φ se mantiene prácticamente en 0 ◦ durante la totalidad del tiempo de simulación, por las razones ya comentadas en el Capítulo 4. Sin embargo, no se aprecian diferencias significativas con respecto a la evolución de la variable ψ en los dos casos bajo estudio y en la situación nominal mostrada en el Capítulo 4, así como tampoco se observan en las diferentes componentes de la velocidad angular de la vela.
68 Capítulo 5. Análisis paramétrico de los resultados obtenidos con la ley de control óptimo (a) ¯ αn f = 2,5 ◦. (b) ¯ αn f = 10 ◦. Figura 5.6 Evolución de αnvariando el valor de ¯ αn f .
5.2 Análisis de la influencia de la actitud final deseada 69 (a) ¯ αn f = 2,5 ◦. (b) ¯ αn f = 10 ◦. Figura 5.7 Evolución de φyθvariando el valor de ¯ αn f .
70 Capítulo 5. Análisis paramétrico de los resultados obtenidos con la ley de control óptimo (a) ¯ αn f = 2,5 ◦. (b) ¯ αn f = 10 ◦. Figura 5.8 Evolución de ψvariando el valor de ¯ αn f . Como se puede ver en la figura 5.10 al aumentar el valor de αn que se quiere tener al final de la maniobra, el módulo de los valores de las variables de control aumenta, algo lógico pues mayor debe ser la diferencia de potencial eléctrico entre los distintos cables para alcanzar un valor de αn mayor. Además, se ha comprobado, para el caso en el que ¯ αn f = 10 ◦ , que las potencias de control no alcanzan su valor de saturación. Dicho valor de saturación se ha definido como 1,15 ·σ y se alcanza para ¯ αn f = 15 ◦ aproximadamente, a partir de dicho valor, las variables del problema dejan de comportarse como lo hacen en la mayoría de casos aquí mostrados y no alcanzan los valores que se desea que alcancen en régimen permanente.
5.3 Análisis de la influencia del tiempo de maniobra 71 (a) ¯ αn f = 2,5 ◦. (b) ¯ αn f = 10 ◦. Figura 5.9 Evolución de las velocidades angulares variando el valor de ¯ αn f . 5.3 Análisis de la influencia del tiempo de maniobra Se estudia a continuación la influencia que tiene el tiempo de realización de la maniobra de actitud, tf sobre la evolución de las distintas variables del problema. Igual que en casos anteriores, únicamente se modifica tf , el resto de parámetros serán iguales al caso nominal. Más concretamente se hace el análisis para los casos en los que se tiene que tf = 60 s y tf = 240 s. Como se puede ver en las figuras 5.11 y 5.12 al disminuir el tiempo de maniobra, es decir, al hacer la maniobra de actitud más rápida, mayor es la amplitud de las oscilaciones que se producen en las variables αn y θ , al pasar de la ley de control (4.52), en la fase de tracking, a la (4.47) en la fase de regulación, lo que
72 Capítulo 5. Análisis paramétrico de los resultados obtenidos con la ley de control óptimo (a) ¯ αn f = 2,5 ◦. (b) ¯ αn f = 10 ◦. Figura 5.10 Evolución de las variables de control variando el valor de ¯ αn f . no permite llegar a la actitud final deseada en el tiempo indicado. Sin embargo, al aumentar el tiempo de maniobra, dicha oscilación disminuye en amplitud, pudiendo incluso desaparecer como se muestra en la figura 5.11b. Por otro lado, al igual que en los casos analizados anteriormente, la modificación del tiempo de maniobra no provoca variaciones en la evolución de ψ. Con respecto a las velocidades angulares en la figura 5.14 se puede visualizar como al aumentar el tiempo de maniobra la oscilación que se produce en Ωy se suaviza disminuyendo la amplitud, además, el ajuste con respecto a su referencia ya no es tan exacto como lo es en el caso en el que el tiempo de maniobra es de 60 s. Mientras que con el resto de componentes no hay cambios significativos.
5.3 Análisis de la influencia del tiempo de maniobra 73 (a) tf= 60 s. (b) tf= 240 s. Figura 5.11 Evolución de αnvariando el valor de tf. Por último, con respecto a las variables de control, decir que los órdenes de magnitud de las mismas no cambian de un caso a otro, pues se quiere llegar en ambos a la misma actitud final. Sin embargo, la pendiente de las potencias de control durante la fase de tracking es mayor en el caso en el que tf = 60 s que cuando tf = 240 s, algo evidente pues en el primero de los casos se desea llegar a la actitud final en un intervalo de tiempo menor y por tanto, las potencias de control tienen que alcanzar su valor de régimen permanente mucho más rápido. Además en la figura 5.15b se puede ver como el escalón que se produce en el tránsito de una ley de control a la otra ha disminuido, esta es la razón por la que en las figuras 5.11b y 5.12b a penas se observa la oscilación en las evoluciones de αnyθ.
74 Capítulo 5. Análisis paramétrico de los resultados obtenidos con la ley de control óptimo (a) tf= 60 s. (b) tf= 240 s. Figura 5.12 Evolución de φyθvariando el valor de tf.
5.3 Análisis de la influencia del tiempo de maniobra 75 (a) tf= 60 s. (b) tf= 240 s. Figura 5.13 Evolución de ψvariando el valor de tf.