Análisis modal y simulación numérica de un modelo plano de E-sail
Abstract
En este proyecto se ha desarrollado un modelo dinámico de E-sail, consistente en un cuerpo rígido central al que se articulan un número variable de tetheres, cuerpos flexibles muy esbeltos en los que se han tenido en cuenta deformaciones axiales y a flexión, incluyendo disipación y actuación. Se presentará un análisis modal del mismo y un análisis de estabilidad en torno a una trayectoria de referencia operativa y realista, variando el número de tetheres del modelo para observar las posibles diferencias en sus dinámicas. Tras mostrar el análisis de estabilidad que el sistema presenta una inherente inestabilidad en torno al equilibrio escogido, se realizaron las correspondientes simulaciones numéricas del modelo no lineal, integrando sus ecuaciones de movimiento y comprobando que efectivamente el sistema no podía seguir la referencia por sí mismo debido a la inestabilidad que presentaba. Por último, con el objetivo de modificar la dinámica del sistema y conseguir la estabilidad, se aplicaron los conceptos de control LQR en el modelo, logrando que el sistema se mantuviera en la trayectoria de referencia. Además, se discutió la robustez del controlador ante perturbaciones y se calcularon las magnitudes asociadas al coste del control.
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 Análisis modal y simulación numérica de un modelo plano de E-sail Autor: Juan Manuel Huertas Nogales Tutores: Daniel García Vallejo y Guillermo Pacheco Ramos Dpto. Ingeniería Mecánica y Fabricación Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2024
Trabajo Fin de Grado Grado en Ingeniería Aeroespacial Análisis modal y simulación numérica de un modelo plano de E-sail Autor: Juan Manuel Huertas Nogales Tutores: Daniel García Vallejo y Guillermo Pacheco Ramos Catedrático de Univerdad y Profesor Sustituto Interino Dpto. Ingeniería Mecánica y Fabricación Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2024
Trabajo Fin de Grado: Análisis modal y simulación numérica de un modelo plano de E-sail Autor: Juan Manuel Huertas Nogales Tutor: Daniel García Vallejo y Guillermo Pacheco Ramos 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 C omo no podría ser de otra forma, en el primer lugar de estos agradecimientos deben de aparecer los nombres de mis padres, Manolo y Puri. Gracias no solo por brindarme la oportunidad de estudiar, con todo el sacrificio y esfuerzo que eso os conlleva, si no también por saber entenderme durante el proceso, apoyarme en los momentos más difíciles y por confiar en mí cuando ni yo mismo lo hacía. Gracias por la educación que me habéis dado, no solo por inculcarme que el trabajo duro es la única manera de conseguir nuestras metas, si no también por mostrarme que los estudios no lo son todo si no se es buena persona, generoso, simpático, honesto y leal. Como padres, habéis cumplido vuestra misión con creces, y yo como hijo espero estar haciéndolo la mitad de bien. Os lo debo todo y como mis logros son los vuestros, este título en tan mío como vuestro. También, tengo que mencionar a mis hermanas. A mi hermana mayor, Noemi, por ser la mejor referente que un niño pudo tener en su día, y que hoy por hoy sigue enseñándome tantas cosas, pues eres de las personas más buenas, generosas y trabajadoras que conozco, y no lo digo yo por ser tu hermano, pues solo hace falta escuchar a todo aquel que te conoce hablar de ti. Y a mi hermana menor, Sofía, por impulsarme a ser un ejemplo igual de bueno que el que yo tuve, pues no te mereces menos. Nunca olvides que estamos muy orgullosos de ti y sabemos que vas a llegar a grandes cosas, aunque por mucho que crezcas siempre vas a seguir siendo nuestra niña. A mi abuelo Melchor, por motivarme a seguir en cada momento y por el orgullo que siente por sus nietos. Y a los que ya no están, sé que desde donde estén me mirarán con mucho orgullo. Por supuesto, gracias también al resto de mi familia. Se me hace también necesario mencionar a mis amigos. A los de mi pueblo, los de toda la vida, los que se sienten orgulloso de mí, que se han alegrado tanto como yo por cada logro y los que tantas veces han tenido que escuchar un: "no puedo, tengo que estudiar" y al final han acabado convenciéndome para despejarme un ratito. Ellos saben quienes y que les estoy infinitamente agradecido. Y también a los nuevos que estos cuatro años y esta carrera me han regalado. Gracias por compartir conmigo cada clase, cada momento de agobio, cada llanto, pero también cada risa, cada plan, cada fiesta... Me lleva a mucha gente y muy buena, y sé que esto será para toda la vida. En especial, gracias a tí María. Quiero agradecer también a mis tutores el haberme ayudado con este proyecto. A Dani, gracias por sacar siempre el tiempo para las infinitas tutorías que te he pedido y por ayudarme con todos y cada uno de los problemas que han surgido, sin duda uno de los mejores profesores con los que me he topado en esta escuela. Y a Guillermo, por no dudar en incorporarse al trabajo y enriquecerlo aún más. Por último, quiero agradecer a aquel chaval de 18 años que decidió empezar esta carrera. Gracias por no rendirte a pesar de que los inicios no fueron de color de rosa, por aguantar y esforzarte. Hoy, salgo con 22 años hecho adulto y como una persona más paciente, luchadora y perseverante gracias a lo que esta carrera me ha enseñado, a pesar de todo lo que la vida tiene aún por enseñarme. I
Resumen E n este proyecto se ha desarrollado un modelo dinámico de E-sail, consistente en un cuerpo rígido central al que se articulan un número variable de tetheres, cuerpos flexibles muy esbeltos en los que se han tenido en cuenta deformaciones axiales y a flexión, incluyendo disipación y actuación. Se presentará un análisis modal del mismo y un análisis de estabilidad en torno a una trayectoria de referencia operativa y realista, variando el número de tetheres del modelo para observar las posibles diferencias en sus dinámicas. Tras mostrar el análisis de estabilidad que el sistema presenta una inherente inestabilidad en torno al equilibrio escogido, se realizaron las correspondientes simulaciones numéricas del modelo no lineal, integrando sus ecuaciones de movimiento y comprobando que efectivamente el sistema no podía seguir la referencia por sí mismo debido a la inestabilidad que presentaba. Por último, con el objetivo de modificar la dinámica del sistema y conseguir la estabilidad, se aplicaron los conceptos de control LQR en el modelo, logrando que el sistema se mantuviera en la trayectoria de referencia. Además, se discutió la robustez del controlador ante perturbaciones y se calcularon las magnitudes asociadas al coste del control. III
XÍndice 4.1.2 Caso Nt=645 4.1.3 Caso Nt=12 47 4.2 Simulaciones con disipación 48 4.2.1 Caso Nt=250 4.2.2 Caso Nt=655 5 Control de la E-sail 59 5.1 Concepto de control LQR (de horizonte infinito) 59 5.2 Aplicación del control LQR 61 5.3 Simulaciones numéricas sin disipación 67 5.3.1 Caso Nt=2 67 5.3.2 Caso Nt=671 5.3.3 Caso Nt=12 73 5.4 Simulaciones numéricas con disipación. 73 5.5 Robustez y coste del control 75 5.5.1 Robustez del controlador 75 5.5.2 Coste del control 79 6 Conclusiones finales y líneas futuras 85 6.1 Conclusiones finales 85 6.2 Líneas futuras y de mejora 86 Apéndice A Variación de los autovalores según los parámetros 87 Apéndice B Respuestas con control 107 Índice de Figuras 111 Índice de Tablas 117 Bibliografía 119
Notación qCoordenadas generalizadas del modelo ˙ qVelocidades generalizadas ¨ qAceleraciones generalizadas x Desplazamiento horizontal del centro de gravedad del vehículo central y Desplazamiento vertical del centro de gravedad del vehículo central θGiro del vehículo central βi Giro del tether i respecto de su dirección radial indeformada uiExtensión del extremo del tether i viFlexión del extremo del tether i φi Coordenada que localiza cada tether en el sistema de referencia solidario al vehículo central NtNúmero de tetheres NcNúmero de coordenadas del sistema Tvc Energía cinética del vehículo central mvc Masa del vehículo central Ivc Inercia del vehículo central ρRadio del vehículo central mvc Matriz de masa del vehículo central ri Vector de posición del punto de articulación del tether i R Vector de posición del centro de gravedad del vehículo central AθMatriz de giro de la coordenada θ ρi Vector de posición del punto de articulación del tether i en el sistema de referencia ligado al vehículo central rP iVector de posición de cualquier punto P del tether i Aβi Matriz de giro para pasar del sistema solidario al tether i al sistema global rP 0 Vector de posición de cualquier punto P del tether i en el sistema de referencia ligado al tether S(η)Matriz de funciones de forma saFunción de forma para la extensión sfFunción de forma para la flexión XI
XII Notación qiCoordenadas del sistema local ligado al tether i TiEnergía cinética del tether i mtiMasa del tether i LiLongitud del tether i mpiMasa puntual del tether i Bq Matriz boolena que permite obtener la coordenada qi de q ˙ rP iVelocidad de cualquier punto del tether i JP i(η,q)Matriz jacobiana del tether i miMatriz de masa del tether i TEnergía cinética del sistema completo mMatriz de masa del sistema completo ViEnergía potencial del tether i EMódulo elástico del material del tether AÁrea tranversal del tether IInercia del tether VEnergía potencial del sistema completo ψi aVector de funciones de forma axiales del tether i ψi fVector de funciones de forma a flexión del tether i kiMatriz de rigidez del tether i kMatriz de rigidez del sistema completo LLagrangiano fFuerzas del sistema QVector de fuerzas generalizadas q0Solución de referencia ωVelocidad angular del vehículo central tTiempo δWTrabajo virtual δrDesplazamientos virtuales ∆l Extensión de los extremos debido a la fuerza centrífuga ˜ q Vector de pequeñas perturbaciones en posición alrededor del movimiento de referencia ˙ ˜ q Vector de pequeñas perturbaciones en velocidad alrededor del movimiento de referencia ¨ ˜ q Vector de pequeñas perturbaciones en aceleración alrededor del movimiento de referencia R Matriz de rotación para pasar al sistema solidario al vehículo central cuando rota a ωconstante q Vector de posición en el sistema local solidario al vehículo central ˙ q Vector de velocidades en el sistema local solidario al vehículo central ¨ q Vector de aceleraciones en el sistema local solidario al vehículo central m0Matriz de masa linealizada c0Matriz de amortiguamiento linealizada k0Matriz de rigidez linealizada m Matriz de masa en el sistema local solidario al vehículo central fFuerzas en el sistema local solidario al vehículo central
Notación XIII m0 Matriz de masa linealizada en el sistema local solidario al vehículo central c0 Matriz de amortiguamiento linealizada en el sistema local solidario al vehículo central k0 Matriz de rigidez linealizada en el sistema local solidario al vehículo central x Vector de estados del sistema dinámico en primer orden A Matriz de estados del sistema dinámico en primer orden φAutovectores λAutovalores MParámetro que representa la masa del vehículo central a Parámetro adimensional que representa la masa de los tetheres b Parámetro adimensional que representa las masas puntuales c Parámetro adimensional que representa la inercia del vehículo central d Parámetro adimensional que representa el radio del disco eParámetro adimensional que representa la rigidez axial z Parámetro adimensional que representa la rigidez a flexión TPeriodo de la E-sail FRFunción de disipación de Rayleigh cfCoeficiente de amortiguamiento a flexión caCoeficiente de amortiguamiento axial ζfRelación de amortiguamiento a flexión ζaRelación de amortiguamiento axial k∗ fRigidez a flexión generalizada k∗ aRigidez axial generalizada m∗ tMasa del tether generalizada BMatriz de control del sistema linealizado uVector de control CMatriz de controlabilidad J(t)Función de coste cuadrática QMatriz de pesos de los estados del sistema RMatriz de pesos de las entradas del sistema KrMatriz de ganancia del controlador óptimo τVector de actuación Fx Módulo de la fuerza que actúa en el vehículo central según x Fy Módulo de la fuerza que actúa en el vehículo central según y Tvc Par en el vehículo central FPiMódulo de la fuerza aplicada en el extremo del tether i ti Vector que proporciona la dirección normal al extremo del tether ien cada instante ti Vector que proporciona la dirección normal al extremo del tether i en cada instante en el sistema solidario al vehículo central αtGiro del extremo de los tether
XIV Notación τ0Vector de actuación de referencia Ru Matriz de giro para pasar τ al sistema solidario al vehículo central τ Vector de actuación en el sistema local solidario al vehículo central ˜ τVector de actuación en el sistema local perturbado τcont Ley de control BxMatriz Ben el sistema de segundo orden βpert Vector de perturbaciones en βi upert Vector de perturbaciones en ui vpert Vector de perturbaciones en vi Pcont(t)Potencia instantánea asociada al control Wcont Trabajo asociado al control max(W) Trabajo máximo desarrollado por la actuación para un barrido de perturbaciones min(W) Trabajo mínimo desarrollado por la actuación para un barrido de perturbaciones max(Pmax) Máximo de la potencia máxima desarrollada por los actuadores para un barrido de perturbaciones min(Pmax) Mínimo de la potencia máxima desarrollada por los actuadores para un barrido de perturbaciones max(Fmax) Máximo de la fuerza máxima desarrollada por los actuadores para un barrido de perturbaciones min(Fmax) Mínimo de la fuerza máxima desarrollada por los actuadores para un barrido de perturbaciones
1 Introducción E n este primer capítulo introductorio, el objetivo principal es presentar el concepto de E-sail de manera comprensible, explicando tanto su funcionamiento básico como sus aplicaciones prácticas. Además, se detallarán los componentes fundamentales que forman parte de esta tecnología. Este conocimiento preliminar es importante de cara a una mejor comprensión posterior acerca del estudio de su dinámica. 1.1 Concepto de E-sail y funcionamiento El concepto de E-sail fue introducido por Pekka Jahunnen en 2004. La E-sail, o vela eléctrica, es un tipo de propulsión espacial que propone utilizar una malla compuesta por alambres conductores delgados (tetheres) dispuestos perpendicularmente al flujo del viento solar. Esta malla se mantiene a un potencial positivo, lo que genera un campo eléctrico que interacciona con los protones del viento solar. Los protones son repelidos por esta barrera eléctrica, ejerciendo una presión que puede ser aprovechada para propulsar la nave espacial. El funcionamiento básico de la E-sail es relativamente simple: al mantener los alambres a un potencial positivo adecuado, se crea un campo eléctrico que afecta a los protones del viento solar, desviándolos y generando una fuerza de empuje. Esto permite acelerar la nave espacial sin necesidad de llevar propelente adicional, lo cual es una ventaja significativa para misiones de larga duración dentro del sistema solar. Para mantener la eficiencia de la E-sail, es crucial expulsar continuamente los electrones que impactan la estructura de la malla. Esto se logra mediante la conexión de la malla a una fuente de voltaje y un emisor de electrones. Aunque el concepto es prometedor para misiones interplanetarias y puede reducir significativamente la masa total requerida en comparación con sistemas propulsores tradicionales, todavía enfrenta desafíos técnicos significativos, como el despliegue y la estabilidad de la malla en el espacio exterior.[4] Es por esto que en este proyecto se ha definido y estudiado un modelo dinámico de E-sail, para contribuir al estudio de esta novedosa tecnología con el fin de que en un futuro pueda llevarse a la realidad con éxito. 1.2 Arquitectura de la E-sail Como ya se ha comentado, la E-sail genera el empuje a partir de la interacción electrostática entre los tetheres cargados positivamente y los iones en el plasma del viento solar. Lo que distingue a una vela solar eléctrica del resto de vehículos espaciales es su ligereza y la capacidad de orientar el vector empuje (dentro de ciertos límites). 1
2Capítulo 1. Introducción Dependiendo de los requerimientos de la misión, la E-sail puede variar su orientación en el espacio para ajustar la dirección del empuje. Esto se consigue modulando el potencial de cada tehter de forma sincronizada con la rotación de la vela. Sin embargo, al realizar los cambios en la actitud de la vela es necesario tener en cuenta otros fenómenos que intervienen en el proceso como pueden ser los cambios en la velocidad angular de los tetheres a causa de dicho cambio, las fuerzas de Coriolis... Por tanto, es necesario disponer de una configuración de E-sail adecuada para evitar así colisiones entre los componentes. La disposición típica de una vela solar eléctrica consiste en una serie de correas delgadas, largas y conductoras (los tetheres ya mencionados) y un cuerpo central en que se encuentra la nave espacial. Cada extremo de los tetheres dispone de una unidad remota, que generalmente está equipada con un pequeño propulsor para actuar sobre el sistema. Las unidades remotas se mantienen unidas entre sí mediante tetheres auxiliares de tal manera que forman un círculo centrado en el cuerpo principal de la nave espacial.[1] Esta configuración se muestra en la figura 1.1, extraida de [1]. Figura 1.1 Arquitectura de la E-sail. En el presente texto, se ha desarrollado un modelo dinámico que incluye el vehículo central de la figura 1.1 así como los tetheres.
1.3 Aplicaciones de la E-sail 3 1.3 Aplicaciones de la E-sail Una valiosa característica de la E-sail es que puede proporcionar empuje en casi cualquier parte del sistema solar fuera de la magnetosfera. Las únicas restricciones son que la dirección del empuje no puede variarse más de aproximadamente 30 grados y que en la magnetosfera de los planetas gigantes se necesitan consideraciones especiales. Se tienen 5 principales aplicaciones en las que la E-sail puede dar resultados satisfactorios. En primer lugar, el acceso a asteroides y planetas cercanos a la Tierra. La E-sail necesita el viento solar u otra corriente de plasma para funcionar, por lo que no se podrá usar dentro de la magnetosfera terrestre, donde el plasma no fluye rápidamente. Sin embargo, planetas como Marte y Venus no disponen de magnetosfera, lo que convierte a la E-sail en una potencial tecnología para la llegada a estos. Las ventajas que puede ofrecer la E-sail respecto a otros métodos de propulsión es principalmente un ahorro económico y en algunos casos de tiempo. En segundo lugar, su uso en órbitas no Keplerianas. Las ventajas que ofrece la vela solar sobre las naves con propulsión tradicional es la capacidad para mantenerse en posiciones no gravitacionalmente estables con un uso mínimo o nulo de propulsante. En tercer lugar, su aplicación para misiones cercanas al Sol. Estas misiones son sobre todo de observación y la E-sail puede proporcionar empuje constante para mantener una órbita cerca del Sol o realizar maniobras de acercamiento o alejamiento repetidamente. En cuarto lugar, su uso para misiones de espacio profundo. El empuje de la E-sail decrece más lentamente que el de otras formas de propulsión, por lo que puede ser una buena opción para misiones largas que necesitan alcanzar grandes distancias sin consumir mucho propulsante. Por último, sus usos como impactadores, penetradores, "data clippers" y su utilización para el transporte de recursos in situ (ISRU). La E-sail puede emplearse para impactar o aterrizar en cuerpos celestes con objetivos científicos o para desviar asteroides. Los "data clippers" podrían transportar información recopilada por otras sondas de regreso a la Tierra como carga útil. Además, la E-sail puede facilitar la utilización de recursos in situ al posibilitar el transporte de materiales entre distintos cuerpos celestes del sistema solar, contribuyendo a la construcción de bases o la extracción de recursos en el espacio.[3] 1.4 Contextualización del proyecto Una vez se ha presentando el concepto de E-sail y sus potenciales aplicaciones, en el resto del documento se va a realizar un estudio de esta tecnología desde un punto de vista fundamentalmente dinámico, no desde la perspectiva de sistema propulsivo ni desde el punto de vista de la mecánica orbital, que serían estudios posteriores al presente. Para ello, se ha desarrollado un modelo dinámico del cual se han hecho estudios modales y de estabilidad junto con las respectivas simulaciones numéricas para ver su comportamiento. Además, se ha dedicado un capítulo completo al estudio del control del modelo.
2 Descripción del modelo E l objetivo de este capítulo es la descripción del modelo dinámico de E-sail que se ha construido, en el que se incluyen cuerpos tanto rígidos como flexibles y el que se tendrán en cuenta tanto deformaciones axiales como a flexión, considerando un modelo plano debido a la complejidad del problema real. En este capítulo, además del modelo dinámico en forma de ecuaciones diferenciales no lineales, se presentará la linealización de este en torno a una trayectoria de referencia de mucho interés operativo que puede seguir la E-sail y en el que se centrará este documento. 2.1 Obtención de las ecuaciones del movimiento 2.1.1 Introducción teórica. Existen muchos sistemas reales que se pueden modelar con una buena aproximación mediante sistemas de un solo grado de libertad. Sin embargo, hay muchos otros en lo que esto no es posible, como será el caso del modelo que se está analizando. El comportamiento de los sistemas dinámicos con un número finitos de grados de libertad se puede representar mediante un sistema de ecuaciones diferenciales de segundo orden de tal manera que el número de ecuaciones será igual al número de grados de libertad. En el caso de vibraciones libres, se tendrá un sistema de ecuaciones homogéneas, de cuya resolución se obtendrán las frecuencias naturales y modos de vibración. El planteamiento para obtención de las ecuaciones del movimiento para los sistemas de varios grados de libertad sigue principalmente dos enfoques: • Newtoniano: En los casos más simples es fácil la obtención de las ecuaciones del movimiento a partir de las ecuaciones de equilibrio de fuerzas en los distintos cuerpos. •Lagrangiano: En estos casos se hace uso de métodos basados en la dinámica analítica. Debido a la complejidad del modelo definido se seguirá un enfoque Lagrangiano. Las ecuaciones de Lagrange para un sistema de N grados de libertad toman la forma d dt ∂L ∂˙ q−∂L ∂q+∂FR ∂˙ q=Q(2.1) donde L=T−V es el lagrangiano, siendo T la energía cinética y V la energía potencial, FR es la función de disipación de Rayleigh, q son las coordenadas generalizadas y Q las fuerzas generalizadas. 5
12 Capítulo 2. Descripción del modelo •B0=1 0 0 ... 0 0 1 0 ... 02×Nc posibilita obtener el vector Rcomo R=B0q •BA=0000100... 0 0000010... 02×Nc permite obtener el vector qA=BAq •BβA=00010... 01×Ncde igual modo permite obtener la coordenada βAdel vector de coordenadas. De este modo, la velocidad del punto Atoma la forma: rA=B0˙ q+∂Aθ ∂θ ρABθ˙ q=B0+∂Aθ ∂θ ρABθ˙ q(2.29) Por otra parte, antes de desarrollar 2.27 para calcular la velocidad de cualquier punto del tether es conveniente expresar la matriz ˙ AβAde la siguiente forma Aα=cos α−sin α sin αcos α−→ ∂Aα ∂α =−sin α−cos α cos α−sin α ˙ AβA=∂Aα ∂α α=θ+φA+βA·(˙ θ+˙ βA) = ∂Aα ∂α α=θ+φA+βA·(Bθ+BβA)˙ q con lo que ˙ rP A="B0+∂Aθ ∂θ ρABθ+∂Aα ∂α α=θ+φA+βA (rP 0+S(η)qA)(Bθ+BβA)+AβAS(η)BA#˙ q=JP A(η,q)·˙ q (2.30) donde JP A(η,q) es la matriz jacobiana que permite conseguir la proporcionalidad de la que antes se habló JP A(η,q) = B0+∂Aθ ∂θ ρABθ+∂Aα ∂α α=θ+φA+βA (rP 0+S(η)qA)(Bθ+BβA)+AβAS(η)BA(2.31) De esta forma se ha llegado a una expresión de la energía cinética del tipo TA=1 2˙ qTmA˙ q donde mA es la matriz de masa del tether Ay toma la expresión mA=ZLA 0 mtA LA (JP A(η,q))TJP A(η,q)dη+mPA(JP A(LA,q))TJP A(LA,q)(2.32) Haciendo uso del procedimiento descrito para cada uno de los tetheres, la energía cinética del sistema es T=Tvc + Nt ∑ i=1 Ti que se puede escribir de forma muy compacta en términos de las matrices de masa del sistema de manera que T=1 2˙ qTm˙ q(2.33) donde m=mvc + Nt ∑ i=1 mi(2.34)
2.1 Obtención de las ecuaciones del movimiento 13 es la matriz de masa del sistema completo. Esta última expresión 2.33 es la que se ha implimentado en MATLAB para el cálculo de la energía cinética del sistema para cualquier número de tetheres que se quiera. 2.1.4 Cálculo de la energía potencial Una vez se ha calculado la energía cinética, el siguiente paso para el desarrollo de las ecuaciones de Lagrange es el cálculo de la energía potencial. Debido al entorno de operación de la E-sail, se ha despreciado el papel de la dinámica orbital y por tanto el potencial gravitatorio, por lo que todas las contribuciones a la energía potencial del sistema serán aportadas como energía potencial de deformación de los tetheres. De esta forma, la energía potencial de un tether cualquiera i=Ase puede calcular como VA=1 2ZLA 0 EA∂uA(η,t) ∂η 2 dη+1 2ZLA 0 EI ∂2vA(η,t) ∂η22 dη(2.35) donde E en el módulo de Young del material, A es el área transversal del tether e I su inercia. Por su parte, u(η,t) = sa(η)uA(t) representa la deformación axial de cualquier punto del tether en cualquier instante de tiempo y v(η,t) = sf(η)vA(t) representa lo mismo para la deformación a flexión. La energía potencial del sistema completa será, por tanto V= Nt ∑ i=1 Vi(2.36) No obstante, se puede escribir de forma más compacta en términos de la matriz de rigidez introduciendo los vectores de funciones de forma ψa y ψf de manera que tomando como ejemplo el tether A se tiene ψA a= 0 0 0 0 sa(η) 0 0 . . . Nc×1 ψA f= 0 0 0 0 0 sf(η) 0 . . . Nc×1 de manera que, para cada tehter, son vectores de la misma dimensión que q que valen 0 en todas las componentes excepto en aquellas en las que q representa una coordenada de deformarción, en la que valdrán saosfsegún se trata de la coordenada de deformación axial o de flexión. La matriz de rigidez de un tether cualquiera será por tanto ki=ZL 0 EA∂ψi a ∂η ∂ψi a ∂η T dη+ZL 0 EI ∂2ψi f ∂η2! ∂2ψi f ∂η2!T dη(2.37) de forma que la matriz de rigidez del sistema es k= Nt ∑ i=1 ki(2.38)
14 Capítulo 2. Descripción del modelo con lo que la energía potencial del sistema completo se puede calcular como V=1 2qTkq (2.39) De nuevo, la expresión 2.39 ha sido la implementada en MATLAB para cualquier número de tetheres que se quiera. 2.1.5 Aplicación de las Ecuaciones de Lagrange Para el caso del estudio de los modos de vibración, al ser estos patrones característicos del propio sistema, se puede considerar que Q=0 en las ecuaciones 2.1. Por otra parte, el amortiguamiento que se introducirá en capítulos posteriores no cambiará los modos de vibración al tratarse de un amortiguamiento proporcional. De esta manera, las ecuaciones 2.1 quedan como d dt ∂L ∂˙ q−∂L ∂q=0(2.40) Como en el sistema de estudio la energía potencial no es función de la velocidad, V=V(˙ q) , las ecuaciones 2.40 toman la siguiente forma d dt ∂T ∂˙ q−∂T ∂q+∂V ∂q=0(2.41) El primer término de la ecuación 2.41 se puede desarrollar para llegar a una expresión más intuitiva. En primer lugar, se manipulará el término de dentro del parántesis del siguiente modo T=1 2˙ qT·(m(q)·˙ q) ∂T ∂˙ q=1 2m(q)·˙ q+1 2˙ qT·∂ ∂˙ q(m(q)·˙ q)(2.42) donde se tendrá en cuenta que la matriz de masa es simétrica, m(q) = m(q)T , debido a la simetría del sistema. A su vez, se puede simplificar el segundo término de la expresión 2.42 de la siguiente forma: 1 2˙ qT·∂ ∂˙ q(m(q)·˙ q) = 1 2vT∂ ∂˙ q(m(q)·˙ q)v=˙ q =1 2∂ ∂˙ q(vTm(q)˙ q)v=˙ q =1 2∂ ∂˙ q(˙ qTmT(q)v)v=˙ q =1 2mT(q)vv=˙ q (2.43) donde se ha hecho uso de que si f(x) = xT·u, entonces su derivada es ∂f(x) ∂x=u. Teniendo en cuenta estos últimos resultados, la ecuación 2.42 se puede desarrollar como 1 2m(q)+ mT(q)˙ q=m(q)·˙ q(2.44) De esta manera, la ecuacion 2.41 queda de la siguiente forma d dt (m(q)·˙ q)−∂T ∂q+∂V ∂q=0 m(q)¨ q+∂ ∂q(m(q)·˙ q)˙ q−∂T ∂q+∂V ∂q=0 m(q)¨ q=f(q,˙ q)(2.45)
2.2 Linealización de las ecuaciones del movimiento 15 donde f(q,˙ q) = −∂ ∂q(m(q)·˙ q)˙ q−∂T ∂q+∂V ∂q(2.46) representa las fuerzas del sistema. De este modo, las ecuaciones 2.45 denotan las ecuaciones dinámicas del sistema que han sido implementadas de manera simbólica en MATLAB de tal forma que se pueden obtener para cualquier número de tetheres. Sin embargo, no se van a adjuntar en el documento en ningún caso debido a complejidad y sobre todo tamaño de las mismas. 2.2 Linealización de las ecuaciones del movimiento El objetivo de este apartado es del linealizar el sistema de ecuaciones anteriormente obtenido alrededor de un movimiento de referencia para posteriormente estudiar la dinámica del sistema en el entorno de éste. El movimiento de referencia escogido es el de la vela girando a velocidad angular constante, sin realizar ninguna translación y con la deformación de los tetheres debida únicamente a la fuerza centrífuga, de modo que solo habrá deformación axial. De esta manera, el vector de coordenadas del movimiento de referencia se puede escribir como: q0= x=0 y=0 θ=ω·t βA=0 uA=∆l vA=0 βB=0 uB=∆l vB=0 . . . Para el cálculo del alargamiento del cable ∆l debido a la fuerza centrífuga se ha hecho uso del Principio de los Trabajos Virtuales. Si se considera por ejemplo el tether A, y se aplica el PTV se debe cumplir: ∂VA ∂uA =Qω(2.47) El trabajo virtual desarrollado por la fuerza centrífuga es δWω=ZLA 0 mtA LA ω2(ρ+η+sauA)δ(ρ+η+sauA)dη+mpAω2(ρ+LA+uA)δ(ρ+LA+uA) =ZLA 0 mtA LA ω2(ρ+η+sauA)sadη+mpAω2(ρ+LA+uA)δuA De este modo, la fuerza centrífuga desarrollada responde a la expresión Qω=ZLA 0 mtA LA ω2(ρ+η)sadη+(ρ+LA)mpAω2+ZLA 0 mpA LA ω2s2 adη+mpAω2uA Por otra parte, ∂VA ∂uA =EA LA uA(2.48)
16 Capítulo 2. Descripción del modelo con lo que despejando de la ecuación 2.47 se tiene el valor de ∆l uA=∆l=RLA 0 mtA LAω2(ρ+η)sadη+(ρ+LA)mpAω2 EA LA−hRLA 0 mtA LAω2s2 adη+mpAω2i(2.49) Sustituyendo esta solución de referencia en las ecuaciones del movimiento se obtenía una solución de ceros, con lo que queda comprobado que la trayectoria calculada es solución de las ecuaciones del movimiento. Una vez que se ha definido el movimiento de referencia alrededor del cual se van a linealizar las ecuaciones, se procede al desarrollo de esta. Las ecuaciones a linealizar son las dadas por 2.45. Para ello, en primer lugar se redefinen las coordenadas del sistema de manera que q=q0+˜ q q q ˙ q=˙ q0+˙ ˜ q q q ¨ q=¨ q0+¨ ˜ q q q donde las variables con "tilde" reperesentan pequeñas pertubaciones alrededor del movimiento de referencia. Si se procede ahora a la expansión en serie de Taylor del modelo alredor de q0 m(q0)+ Nt ∑ k=1 ∂m(q0) ∂qkq0 ˜ qk+...!(¨ q0+¨ ˜ q q q) = f(q0,˙ q0)+ ∂f ∂qq0,˙ q0 ˜ q q q+∂f ∂˙ qq0,˙ q0 ˙ ˜ q q q+... (2.50) Teniendo en cuenta que m(q0)¨ q0=f(q0,˙ q0) y que el sumatorio dará lugar a términos de segundo orden, la ecuación 2.50 se simplifica bastante m(q0)¨ ˜ q=∂f ∂qq0,˙ q0 ˜ q q q+∂f ∂˙ qq0,˙ q0 ˙ ˜ q q q m(q0)¨ ˜ q−∂f ∂qq0,˙ q0 ˜ q q q−∂f ∂˙ qq0,˙ q0 ˙ ˜ q q q=0 m0¨ ˜ q+c0˙ ˜ q+k0˜ q=0(2.51) donde: •m0=m(q0)es la matriz de masa del sistema linealizado. •c0=−∂f ∂qq0,˙ q0 es la matriz de rigidez del sistema linealizado. •k0=−∂f ∂˙ qq0,˙ q0 se corresponde con la matriz de rigidez del sistema linealizado. Sin embargo, estas matrices no son constante, si no que dependen del tiempo al tener que aproximar una trayectoria que a su vez depende el mismo. Con el objetivo de tener las matrices linealizadas del sistema invariantes con el tiempo se hará a continuación un cambio de coordenadas al sistema solidario "1" que rota según la ley θ=ω·tcon la E-sail. De este modo q=Rq ˙ q=˙ Rq +R˙ q ¨ q=¨ Rq +2˙ R˙ q+R¨ q
2.2 Linealización de las ecuaciones del movimiento 17 donde R= cos(ωt)−sin(ωt)0... 0 sin(ωt)cos(ωt)0... 0 0 0 1 ... 0 . . .. . .... .... . . 0 0 ... ... 1 Nc×Nc es la matriz de rotación para pasar al nuevo sistema de referencia y q es el vector de posición en el sistema local. De esta forma, la ecuaciones del movimiento en el nuevo sistema de referencia responden a la expresión (RTm(Rq)R)( ¨ Rq +2˙ R˙ q+R¨ q) = f(Rq,˙ Rq +R˙ q) (RTm(Rq)R)¨ q=RT[f(Rq,˙ Rq +R˙ q)−m(Rq)( ¨ Rq +2˙ R˙ q)] m(q)¨ q=f(q,˙ q)(2.52) donde se ha tenido en cuenta que la matriz de masa en el nuevo sistema de referencia es m= RTm(RT)Ry que R−1=RTal ser una matriz ortogonal. Si se linealiza ahora la ecuación 2.52 entorno a un movimiento de referencia de tal manera que q0=q0y los vectores de estado locales son q=q0+˜ q ˙ q=˙ q0+˙ ˜ q ¨ q=¨ q0+¨ ˜ q se llegará a un resultado totalmente análogo a 2.51, es decir, a una ecuación del tipo: m0¨ ˜ q+c0˙ ˜ q+k0˜ q=0(2.53) donde se tiene que •m0=m(q0)es la matriz de masa linealizada en el sistema local. •c0=−∂f ∂qq0,˙ q0 es la matriz de amortiguamiento linealizada en el sistema local. •k0=−∂f ∂˙ qq0,˙ q0 es la matriz de rigidez linealizada en el sistema local. A través de este cambio de coordenadas se ha llegado a un sistema en el que las matrices son constantes en el tiempo, por lo que se podrán calcular los modos de vibración del sistema. Al igual que para el caso no lineal, todo el desarrollo para los casos lineales en ambos sistemas de referencia ha sido implementado en MATLAB de manera simbólica de forma que para cualquier número de tetheres se podrán obtener la matrices del sistema y la ecuaciones linealizadas. De nuevo, no se mostraran en ningún caso debido al excisivo tamaño que presentan.
3 Análisis modal L os objetivos de este capítulo son por un lado la obtención de los modos de vibración del modelo dinámico de E-sail que se ha desarrollado y por otro lado el análisis de la estabilidad dinámica del sistema. Los modos de vibración representan la deformada de los sistemas cuando vibran ante cualquiera de sus frecuencias naturales, que son la frecuencias a las cuales un sistema tiende a vibrar cuando es perturbado y luego se le permite vibrar libremente, sin ninguna influencia externa adicional. Por otra parte, como se verá a lo largo de la sección, estos modos están relacionados con la respuesta del sistema y por tanto con la estabilidad del mismo. La estabilidad de un sistema dinámico se refiere a su capacidad para mantenerse bajo control y regresar a un estado deseado después de ser perturbado. En el problema que se viene desarrollando es fundamental realizar un análisis de estabilidad para que comprobar si la E-sail es capaz de mantenerse o no en la trayectoria de referencia que se ha definido por si misma. 3.1 Obtención de los modos de vibración Para el cálculo de los modos de vibración, en primer lugar se pasará el sistema 2.53 a primer orden introduciendo el siguiente cambio de variables: (˙ ˜ q=v ˙ v=¨ ˜ q=−m−1 0k0˜ q−m−1 0c0v de esta manera, se obtendrá un sistema del tipo ˙ x=Ax, típico de los sistemas de primer orden v ˙ v=0 I −m−1 0k0−m−1 0c0˜ q v ˙ x=Ax (3.1) donde x=˜ q v= ˜ q ˙ ˜ q!A=0 I −m−1 0k0−m−1 0c0(3.2) son el vector de estados del sistema y la matriz de estado del sistema respectivamente. De estaforma, se ha pasado de Nc coordenadas en un sistema de segundo orden a 2·Nc coordenadas en un sistema de primer orden. Para el cálculo de los modos de vibración, se considera la solución como una superposición de soluciones del tipo x(t) = φeλt(3.3) 19
20 Capítulo 3. Análisis modal y sustituyendo esta solución en la ecuación 3.1 se llega a que λφeλt=Aφeλt(3.4) de tal forma que Aφ=λφ (3.5) de manera que se ha llegado a que las soluciones del sistema son justamente los autovalores y los autovectores de la matriz A , pues así lo indica la expresión 3.5 que es exactamente la definición de un problema de autovalores y autovectores. Si se resuelve el problema de autovalores, se obtendrán Nc parejas de autovalores y autovectores que representan las frecuencias naturales y los modos de vibración del sistema respectivamente. En el caso de tener valores complejos, los modos aparecerán por parejas de complejos conjugados y al representarse únicamente las magnitudes de los modos, se representará exclusivamente un modo por pareja pues son exactamente idénticos. Lo mismo pasará en el caso en el que haya multiplicidades en los autovalores. 3.2 Parámetros del sistema y adimensionalización de los mismos. Para el cálculo de los modos, con el objetivo de reducir el número de parámetros de los que depende el sistema y el coste cumputacional, se van a definir una serie de parámetros adimensionales que sustituirán a varias de las magnitudes que intervienen en el problema. Los parámetros se han definido a partir de las estructuras repetitivas que aparecían en las matrices de masa, rigidez y amortiguamiento linealizadas en el sistema local y a través del análisis dimensional de las magnitudes de tal manera que se ha llegado a lo siguiente: 1. Se adimensionalizarán las masa del sistema del sistema a partir de la masa del vehículo central mvc =M 2. Adimensionalización de la masa de los theteres con la masa del vehículo central mt=a·M→a=mt M 3. Adimensionalización de las masas puntuales con la masa del vehículo central mp=b·M→b=mp M 4. Se ha adimensionalizado la inercia del vehículo central a través de su propia masa y de la longitud de los tetheres, que será la variable que adimensionalice las longitudes Ivc =c·ML2→c=Ivc ML2 5. Adimensionalización del radio del vehículo central con la longitud de los tetheres ρ=d·L→d=ρ L
3.2 Parámetros del sistema y adimensionalización de los mismos. 21 6. El parámetro e representa la adimensionalización de la rigidez axial del tether como se puede ver si se realiza un análisis dimensional A=e·Lω2M E→e=EA Lω2M 7. Análogamente al caso anterior, este parámetro representa la adimensionalización de la rigidez a flexión del tether I=z·L3ω2M E→z=EI L3ω2M De esta manera los 6 parámetros adimensionales ( a,b,c,d,e,z ) sustuirán a 7 de los del problema inicial (mt,mp,Id,ρ,A,I,E). Para poder realizar los cálculos numéricos, se van a presentar a continuación en la tabla 3.1 los valores de los magnitudes que se han tomado, que son valores realistas para el tipo de tecnología que se está estudiando según otros autores. Tabla 3.1 Valores numéricos de las magnitudes. Cuerpo Dimensión Valor Altura, hr[m] 2 Vehículo central Radio, ρ[m] 1 Densidad, ρr[kg/m3] 884 Longitud sin deformar, L0[Km] 10 Área de la sección transversal, A[mm2]4.28 ·10−3 Tetheres Inercia, I[mm4]1.47 ·10−6 Densidad, ρt[kg/m3] 7653 Módulo de Young, E[GPa] 70 Masa puntual Masa, mp[kg] 1.5 Por su parte, la velocidad de giro del vehículo central que no ha sido especificada en la tabla, tomará el valor ω=4·10−3rad/s. A la vista de estos valores, los parámetros definidos anteriormente toman los valores especificados en la tabla 3.2. Tabla 3.2 Valores numéricos de los parámetros. Parámetro Valor M[kg]5554.3358 a5.8972 ·10−5 b2.7·10−4 c5·10−9 d10−4 e0.3371 z1.1579 ·10−18 Estos valores se considerarán los valores base o nominales del problema, a expensas de que alguno de ellos pueda cambiar posteriormente.
28 Capítulo 3. Análisis modal Modo 19 Modo 20 Modo 21 Figura 3.17 Modos del 19 al 21. Modo 22 Modo 23 Modo 24 Figura 3.18 Modos del 22 al 24. Modo 25 Modo 26 Modo 27 Figura 3.19 Modos del 25 al 27.
3.3 Análisis modal 29 Modo 28 Modo 29 Modo 30 Figura 3.20 Modos del 28 al 30. Modo 31 Modo 32 Modo 33 Figura 3.21 Modos del 31 al 33. Modo 34 Modo 35 Modo 36 Figura 3.22 Modos del 34 al 36.
30 Capítulo 3. Análisis modal Modo 37 Modo 38 Modo 39 Figura 3.23 Modos del 37 al 39. Modo 40 Modo 41 Modo 42 Figura 3.24 Modos del 40 al 42. Modo 43 Modo 44 Figura 3.25 Modos del 43 al 44. Para este caso, el número de coordenadas del sistema en primer es de 2·Nc=78 con Nc=39 el número de coordenadas del sistema de segundo orden. Seleccionando solamente aquellos modos diferentes entre sí y descartando las parejas de complejo conjugados y multiplicidades, se obtienen finalmente 44 modos diferentes. Se ve como aparecen sobre todo modos como sólido rígido segun la coordenada βi y modos de deformación a flexión.
3.4 Análisis paramétrico 31 3.4 Análisis paramétrico Es en este apartado en el que se realizará el análisis de estabilidad del sistema. Este análisis es fundamental para comprobar que el modelo dinámico que se ha desarrollado es capaz de mantenerse siguiendo el movimiento de referencia por sí mismo, ya que de los contrario se haría necesario recurrir a otras soluciones. En el apartado anterior se expuso que la respuesta del sistema linealizado en primer orden era x(t) = φeλt y venía determinada por los autovalores y los autovectores de la matriz A de 3.1. Se puede ver como el único término de la respuesta que varía en el tiempo es el correspondiente a eλt . En general, se tiene que los autovalores son complejos y por tanto eλt=eat (cos (bt) +isin (bt)) (3.6) De esta forma, si todos los autovalores λk tienen parte real negativa, el sistema será estable ya que la solución exponencial será decreciente con lo que x=0 para t→∞ . Sin embargo, si un único autovalor tiene parte real positiva, el sistema será inestable y divergirá del movimiento de equilibrio en la dirección del correspondiente autovector (modo) inestable conforme pasa el tiempo. En los apartados posteriores se presentarán los autovalores del sistema para el caso base con los valores de las tablas 3.1 y 3.2. En caso de que se encuentren autovalores con parte real positiva, se realizará un estudio paramétrico en el que se variarán tanto varios de los parámetros adimensionales definidos como varios de los paramétros iniciales del sistema que no fueron sustituidos. Este estudio consistirá en variar los parámetros en un ±50 % del valor base definido en las tablas anteriormente mencionadas y calcular los autovalores para cada valor en el que se ha discretizado el intervalo, con el fin de encontrar de este modo configuraciones estables. Este análisis, además de permitir encontrar posibles configuraciones estables, permitirá observar cómo varía cada autovalor en función de los parámetros. En concordancia con el análisis modal, se realizará el estudio de estabilidad para 2,6 y 12 tetheres. 3.4.1 Caso Nt=2 Lo primero que se ha hecho ha sido calcular los autovalores del sistema en el caso al que se viene llamando base, y posteriormente se han representado su parte real e imaginaria en el entorno del lugar de las raíces como se muestra a continuación. -0.08 -0.06 -0.04 -0.02 0 0.02 0.04 0.06 0.08 -0.15 -0.1 -0.05 0 0.05 0.1 0.15 Figura 3.26 Autovalores del sistema para Nt=2.
32 Capítulo 3. Análisis modal Como se puede ver en la figura 3.26, aparecen autovalores con parte real en el semiplano de abcisas positivo, lo que quiere decir que el sistema presenta una inestabilidad inherente para los valores de los parámetros que inicialmente se han definido. Con el objetivo de buscar combinaciones de parámetros que resulten en configuraciones estables (la parte real de ningún autovalor debe aparecer en el semiplano positivo), se va a realizar un estudio de cómo varían los autovalores en función de los valores de los parámetros. En la figura 3.26 se ha asignado un color a cada autovalor, y en la leyenda de la figura aparecen aquellos autovalores con con parte real positiva, de manera que en los próximos apartados se podrá ver si la parte real de estos evoluciona a ser negativa según el barrido de los parámetros. Se va a barrer según los parámetros ω,a,b,e,L ya que se ha considerado que el resto (a excepción del parámetro z que representa la rigidez a flexión y que se comprobó que para el rango de valores en el que se hacía variar no producía ninguna variación en los autovalores debido a lo pequeña que era su magnitud) carece de importancia como para realizar su análisis. No obstante, tal y como se ha abordado el problema se podrían obtener las variaciones de los autovalores con cualquier parámetro y para cualquier número de tetheres. Por otra parte, es necesario aclarar que en los apartados posteriores simplemente se va a representar la variación de todos los autovalores a la vez en el lugar de la raíces. Sin embargo, también se ha representado la variación de cada autovalor con respecto a a cada parámetro posibilitando así una visualización más clara de la tendencia que siguen. Esto ha sido incluido en el Apéndice A. Variación con ω Como ya se ha comentado, se ha hecho variar el paramétro ω correspondiente a la velocidad de giro del vehículo central en un ±50 % de su valor inicial. Figura 3.27 Variación de los autovalores con ωpara Nt=2. Como se puede ver en la figura 3.39 siguen apareciendo autovalores con parte real positiva, con lo que el sistema seguirá comportándose de manera inestable. Variación con a El parámetro a representa la masa del tether adimensionalizada e igualmente se ha hecho variar en el mismo rango de valores que en el apartado anterior, obteniendo igualmente que el sistema no se
3.4 Análisis paramétrico 33 estabiliza, como refleja la figura 3.28. Figura 3.28 Variación de los autovalores con apara Nt=2. Variación con b En este caso, se considera la variación de la adimensionalización de las masas puntuales de manera que tampoco se consigue el comportamiento dinámico del sistema que se desea. Figura 3.29 Variación de los autovalores con bpara Nt=2. Variación con e Para la variación de la adimensionalización de la rigidez axial apenas se obtienen cambios en los autovalores como se puede ver en la figura 3.36 y además siguen apareciendo autovalores con parte real positiva.
34 Capítulo 3. Análisis modal Figura 3.30 Variación de los autovalores con epara Nt=2. Variación con L Finalmente, se ha consideración la variación de la longitud de los tetheres, obteniendo de la misma manera un comportamiento inestable del sistema. Figura 3.31 Variación de los autovalores con Lpara Nt=2. 3.4.2 Caso Nt=6 Al igual que en el caso Nt=2 , en primer lugar se han calculado los autovalores del sistema para los valores de los parámetros dados por las tablas 3.1 y 3.2 y se han representado las partes real e imaginaria para cada uno de ellos como muestra la figura 3.32.
3.4 Análisis paramétrico 35 -0.08 -0.06 -0.04 -0.02 0 0.02 0.04 0.06 0.08 -0.15 -0.1 -0.05 0 0.05 0.1 0.15 Figura 3.32 Autovalores del sistema para Nt=6. Al igual que ocurría en el caso anterior, se tiene que el sistema es inestable para esta configuración. Se pasa a continuación a la realización del análisis paramétrico para comprobar si los autovalores de la leyenda de la figura 3.32 pueden llegar a tener parte real positiva. Variación con ω Se presenta a continuación la variación de los autovalores con respecto al parámetro ω. Figura 3.33 Variación de los autovalores con ωpara Nt=6. Como se puede ver en la figura 3.33 siguen apareciendo autovalores con parte real positiva, con lo que el sistema segirá comportándose de manera inestable.
36 Capítulo 3. Análisis modal Variación con a En este caso se obtiene igualmente que el sistema no se estabiliza, como refleja la figura 3.34. Figura 3.34 Variación de los autovalores con apara Nt=6. Variación con b De nuevo, se obtiene que el sistema no se estabiliza, como plasma la figura 3.35. Figura 3.35 Variación de los autovalores con bpara Nt=6. Variación con e Igualmente, siguen apareciendo autovalores en el semiplano positivo del eje real.
3.4 Análisis paramétrico 37 Figura 3.36 Variación de los autovalores con epara Nt=6. Variación con L En este caso muy pocos autovalores varían y no se consigue la estabilidad del sistema. Figura 3.37 Variación de los autovalores con Lpara Nt=6. 3.4.3 Caso Nt=12 Como se ha venido realizando hasta el momento, se presentan en primer lugar los autovalores para el caso base.
44 Capítulo 4. Simulación numérica no lineal 0 20 40 60 80 -20 -15 -10 -5 0 5 Figura 4.4 Respuesta de la coordenada βA. Como se puede ver en la figura 4.12, la coordenada βA se mantiene en su posición de equilibrio durante los primeros minutos hasta que finalmente se aleja de manera bastante rápida. 0 20 40 60 80 -6 -4 -2 0 2 4 6 8 10 Figura 4.5 Respuesta de la coordenada uA. La coordenada uA se mantiene inicialmente en ∆l hasta que finalmente se aleja del equilibrio aparentemente buscando otro nuevo en torno a cero. 0 20 40 60 80 -2 -1.5 -1 -0.5 0105 Figura 4.6 Respuesta de la coordenada vA.
4.1 Simulaciones sin disipación 45 Tal y como plasma la figura 4.6, aunque inicialmente vA se mantiene en su equilibrio, en un cierto instante comienza a alejarse de manera exponencial. 0 20 40 60 80 -20 -15 -10 -5 0 5 Figura 4.7 Respuesta de la coordenada βB. 0 20 40 60 80 -6 -4 -2 0 2 4 6 8 10 Figura 4.8 Respuesta de la coordenada uB. 0 20 40 60 80 -2 -1.5 -1 -0.5 0105 Figura 4.9 Respuesta de la coordenada βB. Como se puede ver en las figuras 4.17,4.18 y 4.9, las respuestas son exactamente iguales a las del tether A. Se tiene por tanto que el sistema es capaz de seguir la trayectoria de referencia durante un intervalo de tiempo, concretamente durante unos 11 minutos. Tras esto, los propios errores numéricos asociados a la integración junto a la inherente inestabilidad del sistema hacen que el modelo dinámico se aleje de la trayectoria de referencia según muestran las anteriores respuestas, habiendo algunas que presentan un carácter mucho más inestable que otras. 4.1.2 Caso Nt=6 Para la situación Nt=6 se tiene un sistema de 2Nc=42 ecuaciones diferenciales no lineales. Si se resuelven se obtiene respuestas muy similares a las del caso Nt=2 , sin embargo hay algunas diferencias curiosas de comentar. Con respecto a las coordenadas del centro de gravedad del vehículo central x,y se ve en la figura 4.10 como ahora sí sufren desviaciones apreciables de su valor de referencia en la simulación, a expensas de que en los primeros minutos mantenerse en su equilibrio.
46 Capítulo 4. Simulación numérica no lineal Conforme a la coordenada θ , se ve en 4.11 como al inicio cumple con la ley de referencia pero que sin embargo, igual que en el caso de Nt=2 sufre un cambio de pendiente que los desvía de su equilibrio y en este caso además más abrupto que en el anterior. 0 20 40 60 80 -1 0 1 2 3 4 5 6 x y Figura 4.10 Respuesta de las coordenadas x,y. 0 20 40 60 80 0 5 10 15 20 25 30 35 Figura 4.11 Respuesta de la coordenada θ. Con respeto a las coordenadas βi la figura 4.12 plasma como ahora no todas las coordenadas siguen la misma evolución si no que comienzan a aparecer asimetrías en las respuestas, además de ser estas más rápidas. Esto es unos de los cambios dinámicos más significativos con respecto al caso de Nt=2 y se puede explicar a través de la mayor inestabilidad que presenta el sistema al tener un mayor número de tetheres. 0 20 40 60 80 -30 -25 -20 -15 -10 -5 0 5 Figura 4.12 Respuestas de las coordenadas βi. 0 20 40 60 80 -6 -4 -2 0 2 4 6 8 10 Figura 4.13 Respuestas de la coordenadas ui. Las coordenadas de extensión ui , cuyas respuestas se representa en la figura 4.13, presentan un comportamiento similar al caso de dos tetheres. Sin embargo ahora algunas de ellas aumentan su amplitud con el tiempo mientras que otras parecen mantenerse a lo largo de la simulación, pero todas ellas parecen oscilar en torno a cero.
4.1 Simulaciones sin disipación 47 0 20 40 60 80 -2 -1.5 -1 -0.5 0 0.5 1 1.5 2105 Figura 4.14 Respuesta de la coordenada vi. Las coordenadas vi en este caso presentan asimetrías y crecen más rápidamente que para Nt=2 . A pesar que el modelo sigue la trayectoria de referencia durante los primeros aproximadamente 11 minutos, al igual que en el caso anterior, de nuevo las perturbaciones numéricas acaban sacando a la E-sail de su equilibrio. Además, y como el análisis de estabilidad predijo, el sistema se vuelve más inestable conforme se van añadiendo más tetheres. Esto se puede ver reflejado en las figuras que muestran como las desviaciones del movimiento de referencia son mayores y además las respuestas que lo alejan son más rápidas. 4.1.3 Caso Nt=12 En esta situación, se tiene un sistema de 2Nc=78 ecuaciones diferenciales. Su resolución lleva a conclusiones exactamente análogas a las de la situación Nt=6 , al comiendo de la simulación el modelo es capaz de seguir la trayectoria de referencia pero finalmente se aleja de esta debido a las perturbaciones numéricas. Además, las respuestas (que se muestran en las figuras de abajo) son todavía más rápidas y abruptas y se tienen respuestas diferentes en las coordenadas βi y vi según los tetheres. Todo esto es lógico al tener un mayor número de tetheres pues así el análisis de estabilidad mostró que al tener un mayor números de tetheres mayor sería la inestabilidad del sistema. 0 20 40 60 80 -4 -2 0 2 4 6 8 10 x y Figura 4.15 Respuesta de las coordenadas x,y. 0 20 40 60 80 0 5 10 15 20 25 30 35 Figura 4.16 Respuesta de la coordenada θ.
48 Capítulo 4. Simulación numérica no lineal Figura 4.17 Respuesta de las coordenadas βi.Figura 4.18 Respuesta de las coordenadas ui. Figura 4.19 Respuesta de las coordenadas vi. 4.2 Simulaciones con disipación En este apartado se añadirá disipación al sistema y se repetirán las simulaciones numéricas para ver el efecto que tiene su adición en las respuestas. La disipación implica la conversión de energía mecánica en otras formas de energía, principalmente calor, debido a resistencias internas o externas en un sistema. En sistemas dinámicos y estructuras mecánicas, la disipación es crucial para absorber energía que de otro modo causaría oscilaciones indeseadas, vibraciones excesivas o daños estructurales. Introducir métodos de disipación, como amortiguadores o materiales viscoelásticos, permite controlar las respuestas dinámicas del sistema, mejorando su estabilidad y prolongando su vida útil al reducir el desgaste y la fatiga de los materiales. Recordando la expresión general de las ecuaciones de Lagrange que se presentó en el Capítulo 2: d dt ∂L ∂˙ q−∂L ∂q+∂FR ∂˙ q=Q se dijo que FR denotaba la función de disipación de Rayleigh. La función de disipación de Rayleigh representa una función potencial de la velocidad, definida de forma que su derivada parcial respecto a ˙q produce como resultado las fuerzas de amortiguamiento viscoso asociadas a la coordenada generalizada. Su expresión formal es FR=1 2c˙q2(4.4)
4.2 Simulaciones con disipación 49 donde ces el coeficiente de amortiguamiento asociado a la coordenada generalizada o coeficiente de amortiguamiento generalizado. En el caso de un sistema de n amortiguadores, cada uno de ellos con una velocidad de alargamiento ˙ δj, la función de disipación puede expresarse como FR= n ∑ j=1 1 2cj ˙ δ2 j(4.5) Aplicando este concepto al modelo dinámico de E-sail desarrollado, lo que se pretende es introducir la función de disipación de Rayleigh de tal manera que solo existirá amortiguamiento en el caso de que exista deformación, es decir, sólo existirá disipación en el sistema cuando los tetheres se deformen tanto axialmente como a flexión. De esta manera, la función de disipación de Rayleigh tendrá la siguiente forma FR= Nt ∑ i=1 1 2cfi˙v2 i+ Nt ∑ i=1 1 2cai˙u2 i(4.6) donde cfiycairepresentan el coeficiente de amortaguamiento a flexión y axial respectivamente. Para determinar estos coeficientes, se puede hacer uso de la analogía que presenta la energía disipada por deformación con la energía almacenada durante la deformación, es decir, la energía potencial elástica cuya expresión se recuerda que era Vi=1 2ZL 0 EA∂ui(η,t) ∂η 2 dη+1 2ZL 0 EI ∂2vi(η,t) ∂η22 dη(4.7) de manera que se puede puede escribir en términos de unas rigideces generalizadas de la siguiente forma Vi=1 2k∗ aiu2 i+1 2k∗ fiv2 i(4.8) donde k∗ aiyk∗ fison la rigidez generalizada axial y a flexión respectivamente. En vista a las funciones de forma que se han escogido, estas rigideces tomarán el siguiente aspecto k∗ ai=kai EA Lk∗ fi=kfi EI L3(4.9) donde kai y kfi son coeficientes adimensionales que ajustan la rigidez del sistema según las funciones de forma que se hayan escogido. A la vista de la similitud entre las expresiones de la energía potencial y la función de disipación de Rayleigh, se van a definir los coeficientes caiycficomo cai=qai EA Lcfi=qfi EI L3(4.10) de forma que en este caso qai y qfi son coeficientes que ajustan la energía disipada según la rigidez del elemento, y por análisis dimensional se medirán en segundos. Con el fin de encontrar una expresión de los coeficientes de amortiguamiento más intuitiva, se introducirá la relación de amortiguamiento ζ: ζ=c 2√km (4.11) que aplicada al caso de los tetheres quedará como ζai=cai 2pk∗ aim∗ ti ζfi=cfi 2qk∗ fim∗ ti (4.12)
50 Capítulo 4. Simulación numérica no lineal Así, combinando las expresiones 4.10 y 4.12 se llega que qai=L EA2·qk∗ aim∗ tiζaiqfi=L3 EI 2·qk∗ fim∗ tiζfi(4.13) donde m∗ ties la masa generalizada del tether. Para el cálculo de la respuesta del sistema, se va a considerar que todos los tetheres presentan la misma relación de amortiguamiento, rigideces y masas generalizadas y por tanto el mismo coeficiente de amortiguamiento. Para dar valores, en primer lugar se calculan las rigideces y masa generalizada k∗ ai=ZL 0 EA(s′ a)2dη=0.0300 N/m k∗ fi=ZL 0 EI(s′′ f)2dη=3.0870 ·10−19 N/m m∗ ti=ZL 0 mti Ls2 adη+ZL 0 mti Ls2 fdη=0.1864 kg A continuación, una vez que se han calculado los valores de k∗ ai , k∗ fi y m∗ ti y con lo valores de la tabla 3.1 se calcularán tanto los valores de los coeficientes qai y qfi así como los valores de los amortiguamientos cai,cfi escogiendo para ello unos valores de ζai y ζfi adecuados. Los resultados de estos cálculos se muestran en la tabla 4.1. Tabla 4.1 Valores del amortiguamiento. ζiqi[s] ci[kg/s] ζai=0.01 qai=0.0499 cai=0.0015 ζfi=0.01 qfi=4.6622 ·107cfi=4.7974 ·10−12 Se ha considerado ζai=ζfi=0.01 un valor adecuado según el material y la geometría. De esta manera, el sistema a integrar se convierte ahora en el siguiente: (˙ q=v ˙ v=¨ q=m−1(q)·f(q,˙ q)(4.14) donde f(q,˙ q) = −∂ ∂q(m(q)·˙ q)˙ q−∂T ∂q+∂V ∂q+∂FR ∂˙ q(4.15) de tal manera que ha añadido el término de la función de disipación de Rayleigh. Se van a mostrar a continuación las respuestas del sistema en el caso de añadir disipación. En esta ocasión, solo se analizará los caso de Nt=2 y Nt=6 ya que se ha considerado que son lo suficientemente representativos, pues para el resto de caso, añadiendo tetheres, las conclusiones obtenidas van a ser las mismas. 4.2.1 Caso Nt=2 Se muestran a continuación los resultados de la simulación para las misma condiciones de integración que en el caso de sin disipación.
4.2 Simulaciones con disipación 51 0 20 40 60 80 -4 -3.5 -3 -2.5 -2 -1.5 -1 -0.5 010-13 Sin disipación Con disipación Figura 4.20 Respuesta de la coordenada x. 0 20 40 60 80 -1 0 1 2 3 4 510-13 Sin disipación Con disipación Figura 4.21 Respuesta de la coordenada y. Con respecto a las coordenada x,y , las figuras 4.21 muestran que al igual que en el caso Nt=2 la desviaciones de su equilibrio son muy pequeñas. 0 20 40 60 80 0 5 10 15 20 25 Sin disipación Con disipación Figura 4.22 Respuesta de la coordenada θ.
52 Capítulo 4. Simulación numérica no lineal 0 20 40 60 80 -20 -15 -10 -5 0 5 Sin disipación Con disipación Figura 4.23 Respuesta de la coordenada βA. En las figuras 4.22 y 4.23 se puede ver como la adición de disipación no tiene un efecto notable en las coordenadas θ y βA pues se obtiene básicamente la misma respuesta, como es lógico ya que el tipo de disipación que se ha introducido sólo se hará presente en aquellos casos en los exista deformación por sólido elástico. 0 20 40 60 80 -6 -4 -2 0 2 4 6 8 10 Sin disipación Con disipación Figura 4.24 Respuesta de la coordenada uA. Con respecto a la coordenada uA si que se obtiene un efecto notable con la adición de de la disipación, como era de esperar. Según se dijo, los efectos de la disipación se harían notables en aquellos casos en los que exista deformación y al ser uA una coordenada directamente ligada a esta se ve como el amortiguamiento acaba haciendo nulas las amplitudes de las oscilaciones que antes aparecían y además la coordenada encuentra un nuevo equilibrio en el cero. Todo eso se puede ver en la figura 4.24.
4.2 Simulaciones con disipación 53 0 20 40 60 80 -2 -1.5 -1 -0.5 0105 Sin disipación Con disipación Figura 4.25 Respuesta de la coordenada vA. Con respecto a la coordenada vA , la figura 4.25 plasma como a pesar de ser esta una coordenada asociada a la deformación, no se obtiene un cambio notable en la respuesta ya que se sigue una tendencia exponencial y no oscilatoria. 0 20 40 60 80 -20 -15 -10 -5 0 5 Sin disipación Con disipación Figura 4.26 Respuesta de la coordenada βB. 0 20 40 60 80 -6 -4 -2 0 2 4 6 8 10 Sin disipación Con disipación Figura 4.27 Respuesta de la coordenada uB. 0 20 40 60 80 -2 -1.5 -1 -0.5 0105 Sin disipación Con disipación Figura 4.28 Respuesta de la coordenada vB.
60 Capítulo 5. Control de la E-sail A diferecencia del sistema dinámico definido por 5.1, en el definido por 5.4 la dinámica del sistema ya no viene dada exclusivamente los autovalores del la matriz A sino que ya hay que tener en cuenta los términos de la actuación. Ya que el tipo de control que se está considerando es control en bucle cerrado con retroalimentación, la ley de control que se sigue es la siguiente u=−Kx (5.5) donde K es la matriz de ganancia que permite modificar la dinámica del sistema de la siguiente forma ˙ x=Ax +Bu = (A−BK)x(5.6) Se ve como ahora el comportamiento del sistema vendrá dado los autovalores de la matriz (A−BK) , de tal manera que con una elección correcta de K es posible obtener todos los autovalores con parte real negativa y hacer estable el sistema. Para que el sistema se pueda controlar, se debe cumplir el requisito de que la matriz de controlabilidad C debe de tener tener rango máximo de n , donde n es el número de estados del sistema y C toma la expresión C=B AB A2B... An−1B(5.7) Una vez se tiene que el sistema es controlable, es posible manipular los autovalores del sistema (A−BK) mediante una gran cantidad de leyes de control u=−Kx . Sin embargo, hay que tener en cuenta que el intentar conseguir autovalores demasiado estables puede requerir un control excesivamente caro y también puede resultar en señales de actuación que exceden los valores máximos permitidos. Por tanto, es necesario encontrar un balance entre la estabilidad del sistema y la agresividad del control. La función de coste J(t) = Zt 0 xT(τ)Qx(τ)+ uT(τ)Ru(τ)dτ(5.8) equilibra el coste de la regulación efectiva de los estados con el coste del control. Las matrices Q y R ponderan el coste de las desviaciones de los estados respecto al cero y el coste de la actuación respectivamente. Estas matrices suelen ser diagonales, tal que los elementos de la diagonal pueden modificarse para cambiar la importancia relativa de los objetivos del control. La ley de control del controlador lineal-quadratic regulator (LQR) u=−Krx está diseñada para minimizar la función de coste J=l´ ımt→∞J(t) . Se le llama control LQR ya que se utiliza una ley de control lineal para controlar un sistema no lineal obtenida al minimizar una función de coste cuadrática. Como la función de coste es cuadrática, existe una solución analíitica para la ganancia del controlador óptimo Kr, dado por Kr=R−1BTX(5.9) donde Xes la solución de una ecuación algebraica de Riccati: ATX+XA −XBR−1BTX+Q=0(5.10) Resolviendo esta ecuación, se puede obtener la matriz Kr.
5.2 Aplicación del control LQR 61 5.2 Aplicación del control LQR Una vez que se ha especificado en qué consiste el control LQR, se va a procedir a la implementación de un controlador de este tipo en el modelo dinámico de E-sail que se ha desarrollado. En primer lugar se va a introducir la actuación en el sistema. Se va considerar que los actuadores de cada extremo de los tetheres desarrollan una fuerza que se mantiene perpendicular a dicho extremo en cada instante. Además, en el vehículo central actuará un par del control y en su centro de gravedad dos fuerzas según las coordenadas de su centro de gravedad x,y . En un principio, se intentó controlar el sistema sin estas dos últimas fuerzas mencionadas que pueden resultar poco realistas. Sin embargo se obtenía un rango de la matriz de controlabilidad C de Nc−2 , es decir, hacían falta otras dos fuerzas para poder así controlar el sistema. Es por esto que se añadieron estas fuerzas al sistema. Figura 5.1 Actuación del sistema. En la figura 5.1 puede verse cada una de las fuerzas anteriormente mencionadas y su actuación en el sistema. La actuación en las ecuaciones del movimiento del sistema se va a introducir a través del término de las fuerzas generalizadas de Lagrange Q . Para ello en primer se define el vector de actuación de la siguiente forma τ= Fx Fy Tvc FPA FPB . . . Na×1 donde Fx y Fy se corresponden con las fuerzas que actúan en el centro de gravedad del vehículo central, Tvc el par que actúa en el vehículo central y FPA,FPB,... las fuerzas que actúan sobre los extremos de los tetheres, desarrolladas por las unidades remotas. De esta manera, se tendrán un vector de actuación de Na=3+Ntcomponentes. A continuación, se va a calcular la aportación de cada componente del vector de actuación τ al término de las fuerzas generalizadas Q.
62 Capítulo 5. Control de la E-sail El caso de las fuerzas Fx,Fy y Tvc es trivial, de tal manera que la contribución de estas tres fuerzas se calcula de la siguiente forma δWvc =FT vcδrvc (5.11) donde rvc = x y θ Fvc = Fx Fy Tvc (5.12) son el vector de posición del vehículo y el vector de fuerzas que actúan sobre este respectivamente. Si se reescribe δrvc de la siguiente forma δrvc =∂rvc ∂qδq(5.13) y se tiene en cuenta que δWvc =QT vcδq(5.14) se llega a que la fuerza generalizada aportada por la actuación en el vehículo central es Qvc =∂rvc ∂qT Fvc (5.15) donde ∂rvc ∂q= 1 0 0 0 ... 0 0 1 0 0 ... 0 0 0 1 0 ... 0 3×Nc (5.16) De esta manera se llega a que Qvc = Fx Fy Tvc 0 . . . 0 Nc×1 (5.17)
5.2 Aplicación del control LQR 63 Figura 5.2 Fuerza aplicada en los tetheres. Para el caso de las fuerzas aplicadas en los extremos de los tetheres, se seguirá el mismo procedimiento. Si se centra el estudio en un único tether, por ejemplo el tether A, la situación sería la mostrada en la figura 5.2. A la vista de la figura, lo más sencillo es escribir la fuerza FPA como un vector unitario multiplicado por su módulo FPA=tA·FPA(5.18) donde tA es dicho vector unitario que proporciona la dirección perpendicular al extremo del tether en cualquier momento y como ya se sabe FPAes la fuerza de actuación de la unidad remota. El vector tA se puede expresar en términos de su expresión en el sistema de referencia local "2" y de la matriz de giro AβA tA=AβAtA(5.19) donde se recuerda que AβA=cos (θ+φA+βA)−sin (θ+φA+βA) sin (θ+φA+βA)cos (θ+φA+βA) es la matriz de rotación para pasar del sistema de referencia "2" al sistema inercial "0" y tA es el vector en locales que toma la expresión tA=−sin αt cos αt con αt=∂vA(η,t) ∂η η=L =s′f(η=L)·vA el ángulo que forma el extremo del tether con respecto a la dirección definida por uA. Aplicando nuevamente el TTV: δWA P= (FA P)TδrA P(5.20)
64 Capítulo 5. Control de la E-sail donde rA P es el vector de posición que localiza el extremo del tether (punto A) en el sistema inercial "0", que se recuerda que era rP=x y+AβAρ 0+η+sa(η)uA sf(η)vA rA P=rP(η=L) Si se expresa ahora δrA Pcomo δrA P=∂rA P ∂qδq(5.21) y se tiene en cuenta que δWA P= (QA P)Tδq(5.22) se llega a que la aportación de la fuerza aplicada en el extremo del tether A a las fuerzas generalizadas del sistema es QA P="∂rA P ∂qT AβA·tA#·FA P(5.23) Si se repite el procedimiento para todos los tetheres, el vector de fuerzas generalizas será por tanto Q=Qvc + Nt ∑ i=1 Qi P(5.24) Las ecuaciones del movimiento responden ahora, con la introducción de las fuerzas generalizadas, a la siguiente expresión m(q)¨ q=f(q,˙ q)+Q(q,τ)(5.25) Con el objetivo de llegar a un sistema lineal de la forma ˙ x=Ax +Bu , se va a linealizar el sistema 5.25 en torno a q0 , τ0 donde q0 es la trayectoria de referencia que se definió en el Capítulo 2 y τ0 es el valor del vector de actuación cuando el modelo sigue la trayectoria de referencia por sí mismo, es decir τ0= 0 0 . . . 0 Na×1 Antes de pasar a la linealización de las ecuaciones, se van a definir los cambios de coordenadas necesarios para obtener unas matrices A y B constantes en el tiempo. Comenzando por el vector q q=Rq ˙ q=˙ Rq +R˙ q ¨ q=¨ Rq +2˙ R˙ q+R¨ q el cambio de coordenadas para el vector de coordenadas generalizadas es exactamente análogo al cambio de coordenadas que se hizo para el sistema sin actuación, como es lógico. Con respecto al vector τ, el cambio de coordenadas necesario es τ=Ru(t)τ
5.2 Aplicación del control LQR 65 donde τ es el vector de actuación en el sistema local "1" y Ru(t) es la matriz de giro que permite pasar el vector de actuación del sistema inercial "0" al sistema "1" Ru(t) = cos(ωt)−sin(ωt)0... 0 sin(ωt)cos(ωt)0... 0 0 0 1 ... 0 . . .. . .... .... . . 0 0 ... ... 1 Na×Na De esta manera, el sistema 5.25 en el sistema global expresado en el sistema local es, recuperando la ecuación 2.52 m(q)¨ q=f(q,˙ q)+RTQ(Rq,Ruτ) m(q)¨ q=f(q,˙ q)+Q(q,τ)(5.26) Si se redefinen ahora los vectores de estado y el vector de actuación de forma que representen pequeñas variaciones alrededor de la trayectoria de referencia se llega a que q=q0+˜ q ˙ q=˙ q0+˙ ˜ q ¨ q=¨ q0+¨ ˜ q τ=τ0+˜ τ donde las variables con "tilde" son las que representan estas pertubaciones entorno al equilibrio y se tiene que q0=q0,τ0=τ0. Si se desarrolla ahora la ecuación 5.26 entorno a (q0 , τ0)≡0 de igual modo que se hizo en el Capítulo 2 m(q0)¨ ˜ q+m(q0)¨ q0=f(q0,˙ q0)+ ∂f ∂q0 ˜ q+∂f ∂˙ q0 ˙ ˜ q+Q(q0,τ0)+ ∂Q ∂q0 ˜ q+∂Q ∂τ 0 ˜ τ(5.27) Esta ecuación se simplifica teniendo en cuenta que m(q0)¨ q0=f(q0,˙ q0) y además se puede comprobar como Q(q0,τ0) = 0y∂Q ∂q0=0. Así, la ecuación 5.27 toma la expresión, recordando los resultados de 2.53: m0¨ ˜ q+c0˙ ˜ q+k0˜ q=Bx˜ τ(5.28) donde Bx=∂Q ∂τ 0tiene dimensiones Nc×Na. Se está ya en concidiones de pasar el sistema a primer orden. Si definen la matrices A y B como A=0 I −m−1 0k0−m−1 0c0B=0 m−1 0Bx yucomo u=˜ τ, se llega a la expresión ˙ x=Ax +Bu (5.29)
66 Capítulo 5. Control de la E-sail donde no hay que olvidar que el vector de estados es x= ˜ q ˙ ˜ q!(5.30) Una vez se tiene este sistema, se puede diseñar el controlador LQR sin más que utilizar el comando lqr en MATLAB de la siguiente forma: || >> K = lqr(A, B, Q, R); donde A,Bson las matrices del sistema 5.29 y Q,Rlas matrices de coste de desviación respecto a la referencia y de control respectivamente. Como ya se dijo, estas matrices suelen ser diagonales con valores en la diagonal iguales al peso del control y con el objetivo de mantener la trayectoria de referencia se tomaron Q=I2Nc×2NcR=0.01 ·INa×Na donde Ies la matriz identidad. La lógica que se siguió para llegar a estos valores fue la mostrada por el diagrama de flujo representado en 5.3. De esta manera, partiendo de unos valores iniciales muy similares a los finales, se llegó a una combinación de los pesos de Q y R que consiguen un controlador robusto y eficaz. Figura 5.3 Diagrama de flujo. Una vez se tiene la matriz de ganancia K, ya solo queda introducir la ley de control u=−Kx en el sistema 5.25. Sin embargo, hay que tener en cuenta que esta ley se ha calculado en el sistema local perturbado y la ecuación 5.25 está expresada en el sistema global, con lo que habrá que llevar dicha ley al sistema de referencia global. Para ello, se expresará la matriz Kcomo K=K1K2Na×2Nc donde K1=KNa×[1,Nc]yK2=KNa×[Nc+1,2Nc]. De esta forma, la actuación de control en el sistema global se calculará en los siguientes pasos u=−Kx u=−K1˜ q−K2˙ ˜ q(5.31) Si se tiene en cuenta ahora que ˜ q=q−q0y˙ ˜ q=˙ q−˙ q0, la expresión 5.31 queda como u=−K1(q−q0)−K2(˙ q−˙ q0)(5.32)
5.3 Simulaciones numéricas sin disipación 67 Por último, pasando el vector de estados al sistema global deshaciendo la transformación (q=R−1q ˙ q=R−1˙ q−R−1˙ RR−1q y teniendo en cuenta que q0=q0y˙ q0=˙ q0la expresión 5.33 finalmente es u=−K1(R−1q−q0)−K2(R−1˙ q−R−1˙ RR−1q−˙ q0)(5.33) Para poder introducir esta ley en el sistema 5.25, se debe deshacer el cambio u=˜ τ , de tal manera que la ley de control en el sistema global es τcont =Ru(˜ τ+˜ τ0) τcont =Ru−K1(R−1q−q0)−K2(R−1˙ q−R−1˙ RR−1q−˙ q0)) 3(5.34) donde se ha tenido en cuenta que ˜ τ0=0. Introduciendo ahora la ley de control 5.34 en el sistema 5.25, se tiene que m(q)¨ q=f(q,˙ q)+Q(q,τcont )(5.35) El sistema 5.35 es que el se integrará en los posteriores apartados para demostrar que el controlador es totalmente robusto. 5.3 Simulaciones numéricas sin disipación El objetivo de este apartado es demostrar la funcionalidad del controlador diseñado, mostrando como es capaz de controlar el sistema dinámico para cualquier número de tetheres. El sistema a integrar ahora es (˙ q=v ˙ v=¨ q=m−1(q)·f(q,˙ q)+m−1(q)·Q(q,τcont )(5.36) con f(q,˙ q) = −∂ ∂q(m(q)·˙ q)˙ q−∂T ∂q+∂V ∂q Las condiciones iniciales serán nuevamente q0(t=0) y ˙ q0(t=0) y el tiempo de integración 3·T . 5.3.1 Caso Nt=2 En primer lugar, se puede ver como el controlador consigue modificar la dinámica del sistema y estabilizarlo el pues todos los autovalores del sistema ( A−BK ) tienes parte real negativa como se puede ver en la figura 5.4.
68 Capítulo 5. Control de la E-sail -20 -10 0 10 20 -0.15 -0.1 -0.05 0 0.05 0.1 0.15 Figura 5.4 Autovalores con controlador. Se presentan a continuación las respuestas del sistema tras implementar el control LQR. Se puede ver como para todas las respuestas obtenidas, en las que la línea roja representa la respuesta correcta (trayectoria de referencia), tras integrar el sistema 5.36 el controlador consigue que el modelo siga la trayectoria de referencia sin desviaciones apreciables como bien muestran las figuras 5.5 y 5.6, que plasman el orden de las desviaciones de la referencia el cual es insignificante (las amplitudes son minúsculas, menores del orden 10−13 ), con lo que se puede considerar que es robusto y cumple la misión para la que fue diseñado. Las figuras 5.8 y 5.10 son una buena referencia para visualizar de nuevo la robustez del controlador, pues se puede ver como la coordenadas βA y vA , que son las que una mayor inestabilidad mostraban, no se desvían del equilibrio una vez se ha implementado el controlador 0 20 40 60 80 -3.5 -3 -2.5 -2 -1.5 -1 -0.5 0 0.5 10-13 Figura 5.5 Respuesta de x.
5.3 Simulaciones numéricas sin disipación 69 0 20 40 60 80 -0.5 0 0.5 1 1.5 2 2.5 3 3.5 10-13 Figura 5.6 Respuesta de y. 0 20 40 60 80 0 5 10 15 20 25 Figura 5.7 Respuesta de θ. 0 20 40 60 80 -20 -15 -10 -5 0 5 Figura 5.8 Respuesta de βA.
76 Capítulo 5. Control de la E-sail 0 5 10 15 20 25 -4 -2 0 2 4 610-13 tsim=3T tsim=25T Figura 5.23 Respuesta de x. 0 5 10 15 20 25 -6 -4 -2 0 2 410-13 tsim=3T tsim=25T Figura 5.24 Respuesta de y. 0 5 10 15 20 25 0 20 40 60 80 100 120 140 160 tsim=3T tsim=25T Figura 5.25 Respuesta de θ. 0 5 10 15 20 25 -0.05 0 0.05 0.1 0.15 0.2 tsim=3T tsim=25T Figura 5.26 Respuesta de βA. 0 5 10 15 20 25 8.52 8.54 8.56 8.58 8.6 8.62 tsim=3T tsim=25T Figura 5.27 Respuesta de uA. 0 5 10 15 20 25 -0.02 0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 tsim=3T tsim=25T Figura 5.28 Respuesta de vA. La simulación se ha realizado para Nt=2 pero las respuestas correspondientes al tether "B" no se han mostrado pues al ser la perturbaciones simétricas, las respuestas también lo son para este caso. Se ha obtenido la respuesta para dos tiempos de integración: para 3 veces el periodo, como se viene haciendo durante todo el documento, y para 25 veces el periodo. Como ya se ha dicho, se está implementando un controlador de horizonte infinito, y como se puede ver en las figuras, en la mayoría de respuestas se puede ver la diferencia entre los dos tiempos de integración. Lo normal (y lo que sucede) es que cuanto mayor sea sea el tiempo de simulación, mejor sea la respuesta que se
5.5 Robustez y coste del control 77 obtiene, puesto el objetivo del control LQR de horizonte infinito es conseguir una ley de control minimizando el coste J(t) para t→∞ . Si se quisira minimizar esta función y llegar a una ley de control en un tiempo determinado, se tendría que hacer uso del control LQR de horizonte finito. Para cuantificar cuánto mejora la respuesta al aumentar el tiempo de simulación, se ha hecho uso del error relativo que existe entre las posiciones finales calculadas y las exactas: err(%) = |q0(t=tsim)−q(t=tsim)| |q0(t=tsim)|·100 (5.38) de manera que para la integración de 3 periodos en la que tsim =3·T , se obtiene un error del 0.6484 % y para la integración de 25 periodos, en la que tsim =25 ·T, se tiene un error del 2.8470 ·10−4%. Se demuestra que el controlador es robusto para este tipo de perturbaciones a pesar de estar controlando el sistema no lineal y estar introduciendo una perturbación que se puede considerar grande. Repitiendo el proceso para unas perturbaciones en ui , de tal manera que ahora las condiones iniciales son CI ≡q0(t=0)+upert donde upert =00003·∆lm003·∆lm0...T 1×Nc las respuestas que se obtienen son las presentadas en las siguientes figuras. 0 5 10 15 20 25 -6 -5 -4 -3 -2 -1 0 1 210-11 tsim=3T tsim=25T Figura 5.29 Respuesta de x. 0 5 10 15 20 25 -1.5 -1 -0.5 0 0.5 1 1.5 10-11 tsim=3T tsim=25T Figura 5.30 Respuesta de y. 0 5 10 15 20 25 0 50 100 150 tsim=3T tsim=25T Figura 5.31 Respuesta de θ. 0 5 10 15 20 25 -1.5 -1 -0.5 0 0.5 1 1.5 tsim=3T tsim=25T Figura 5.32 Respuesta de βA.
78 Capítulo 5. Control de la E-sail 0 5 10 15 20 25 -20 -10 0 10 20 30 40 tsim=3T tsim=25T Figura 5.33 Respuesta de uA. 0 5 10 15 20 25 -30 -20 -10 0 10 20 30 tsim=3T tsim=25T Figura 5.34 Respuesta de vA. En este caso, se obtienen unos errores algo más elevados, pues tiene un error del 5.7225 % para el caso de tsim =3Ty del 0.0034 % para tsim =25T. Considerando un tiempo de integración mayor se llegaría un menor valor del error. A pesar de esto, de puede considerar que el controlador puede seguir siendo eficaz para conrregir este tipo de perturbaciones. Por último, con respecto a las perturbaciones en vi, las condiciones iniciales son CI ≡q0(t=0)+vpert donde vpert =0 0 0 0 0 300m0 0 300m...T 1×Nc Los resultando se muestran en las siguientes figuras. 0 5 10 15 20 25 -10 -8 -6 -4 -2 0 2 410-12 tsim=3T tsim=25T Figura 5.35 Respuesta de x. 0 5 10 15 20 25 -1.5 -1 -0.5 0 0.5 1 1.5 210-12 tsim=3T tsim=25T Figura 5.36 Respuesta de y.
5.5 Robustez y coste del control 79 0 5 10 15 20 25 0 50 100 150 tsim=3T tsim=25T Figura 5.37 Respuesta de θ. 0 5 10 15 20 25 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 tsim=3T tsim=25T Figura 5.38 Respuesta de βA. 0 5 10 15 20 25 2 4 6 8 10 12 14 16 tsim=3T tsim=25T Figura 5.39 Respuesta de uA. 0 5 10 15 20 25 -50 0 50 100 150 200 250 300 tsim=3T tsim=25T Figura 5.40 Respuesta de vA. Aunque las perturbaciones son bastantes grandes, de 300 m, el controlador es capaz de corregirlas de manera que se obtiene un error del 0.2366 % para la simulación de 3 periodos y del 1.3614 ·10−4% para la de 25. 5.5.2 Coste del control El objetivo de apartado es cuantificar el coste asociado al control de la E-sail. En concreto se ha calculado el trabajo, potencia y fuerzas implicadas en controlar el modelo cuando se aplican perturbaciones de la misma manera que se aplicaron en el apartado anterior. Para calcular la potencia instantánea asociado al control, simplemente hay que multiplicar la ecuación 5.35 por ˙ qT. De esta forma ˙ qT·m(q)¨ q=˙ qT·f(q,˙ q)+ ˙ qT·Q(q,τcont )(5.39) donde Pcont(t) = ˙ qT·Q(q,τcont )(5.40) es la potencia asociada al control. Una vez se tiene la potencia, el trabajo desarrollado por los actuadores es sencillo de calcular y se hace de la siguiente forma Wcont =Zt 0 Pcont(ˆ t)dˆ t=Zt 0 ˙ qT·Q(q,τcont )dˆ t(5.41)
80 Capítulo 5. Control de la E-sail Como el cálculo se realiza computacionalmente, t se corresponderá con el tiempo de la simulación. En este caso se ha efectuado utilizando la orden trapz de MATLAB. Por último, las fuerzas implicadas en el control, o fuerzas desarrolladas por los actuadores, se pueden calcular sin más que evaluar el vector de actuacion τcont =τcont(q,˙ q,t) , de manera que se obtiene la fuerza que se está desarrollando en cada instante de tiempo conociendo la ley de control. Como se puede ver, para calcular todas esta magnitudes es preciso en primer lugar integrar las ecuaciones del movimiento para obtener los vectores q y ˙ q . El procedimiento que se ha seguido ha sido discretizar un intervalo de un cierto tipo de perturbaciones (perturbaciones en βi , ui y vi ) y para cada uno de los valores del intervalo se han resuelto las ecuaciones 5.35 para obtener el vector de posición y velocidad. Con estos valores se ha calculado el trabajo desarrollado para corregir la perturbación ( W ), la potencia máxima instantánea ( max(Pcont(t)) = Pmax ) y la fuerza máxima desarrollada por las unidades remotas ( max(τcont (q,˙ q,t)) = |Fmax| ). Este método se ha seguido para los casos de 2,6 y 12 tetheres para poder ver así las diferencias asociadas al coste de control. Comenzando por las perturbaciones βi , se han tenido en cuenta valores de estas comprendidas entre ±10◦. Los resultados obtenido se muestran en las figuras 5.41,5.42 y 5.43. -10 -5 0 5 10 0 50 100 150 200 Nt=2 Nt=6 Nt=12 Figura 5.41 Trabajo desarrollado según perturbaciones en βi. -10 -5 0 5 10 0 20 40 60 80 100 120 140 Nt=2 Nt=6 Nt=12 Figura 5.42 Potencia máxima desarrollada según perturbaciones en βi.
5.5 Robustez y coste del control 81 -10 -5 0 5 10 0 10 20 30 40 50 Nt=2 Nt=6 Nt=12 Figura 5.43 Fuerza máxima desarrollada según perturbaciones en βi. Se ve que el comportamiento es el esperado: las magnitudes son máximas cuando la pertubación es máxima y mínima cuando son nulas, pues en este caso apenas se requiere de control. Se muestran en la tabla 5.1 valores concretos de las magnitudes calculadas. En ella se puede ver para cada número de tetheres los intervalos entre los que se encuentran la energía, potencia y fuerzas desarrolladas por los actuadores. Tabla 5.1 Valores de las magnitudes según perturbaciones en βi. Perturbaciones en βi Ntmax(W)[J] min(W)[J] max(Pmax)[W] min(Pmax)[W] max(Fmax)[N] min(Fmax)[N] 2 73.4442 1.8651 89.3509 0.0642 17.2828 0.4431 6 129.6165 3.2426 131.3686 0.1540 29.9163 0.7671 12 175.9262 4.2202 122.7631 0.2843 42.1816 1.0816 Con respecto a las perturbaciones según ui , al igual que se hizo en el apartado anterior, se ha añadido a la condición inicial de referencia un vector de perturbaciones tal que a las coordenadas ui se sumarán ±2·∆lm. Los resultados se muestran en la figuras 5.44, 5.45 y 5.46. -20 -10 0 10 20 0 1 2 3 4 5 6 7 8104 Nt=2 Nt=6 Nt=12 Figura 5.44 Trabajo desarrollado según perturbaciones en ui.
82 Capítulo 5. Control de la E-sail -20 -10 0 10 20 0 1 2 3 4 5 6 7 8104 Nt=2 Nt=6 Nt=12 Figura 5.45 Potencia máxima desarrollada según perturbaciones en ui. -20 -10 0 10 20 0 50 100 150 200 Nt=2 Nt=6 Nt=12 Figura 5.46 Fuerza máxima desarrollada según perturbaciones en ui. En cuanto a valores concretos, se muestran en la tabla 5.2. Tabla 5.2 Valores de las magnitudes según perturbaciones en ui. Perturbaciones en ui Ntmax(W)[J] min(W)[J] max(Pmax)[W] min(Pmax)[W] max(Fmax)[N] min(Fmax)[N] 2 13149.362 333.45 12858.0828 7.3177 160.8501 4.1244 6 39462.2778 999.9813 38283.9415 21.6801 159.6394 4.0933 12 78956.4159 2001.819 76563.591 43.2525 164.8681 4.2274 Por último, para las perturbaciones en vi se ha realizado un barrido en ±100 m de pertubación, de manera que los resultados obtenidos se muestran en las gráficas 5.47,5.48 y 5.49.
5.5 Robustez y coste del control 83 -100 -50 0 50 100 0 0.5 1 1.5 2 2.5 3105 Nt=2 Nt=6 Nt=12 Figura 5.47 Trabajo desarrollado según perturbaciones en ui. -100 -50 0 50 100 0 1 2 3 4 5105 Nt=2 Nt=6 Nt=12 Figura 5.48 Potencia máxima desarrollada según perturbaciones en vi. -100 -50 0 50 100 0 200 400 600 800 1000 1200 Nt=2 Nt=6 Nt=12 Figura 5.49 Fuerza máxima desarrollada según perturbaciones en vi.
84 Capítulo 5. Control de la E-sail En cuanto a los valores límite de las magnitudes, se muestran en la tabla 5.3. Tabla 5.3 Valores de las magnitudes según perturbaciones en vi. Perturbaciones en vi Ntmax(W)[J] min(W)[J] max(Pmax)[W] min(Pmax)[W] max(Fmax)[N] min(Fmax)[N] 2 44224.3666 626.9467 80209.267 22.3243 1001.616 14.5162 6 132602.8084 1881.1746 240603.5564 67.16607 1001.5152 14.5147 12 265271.8641 3763.5373 481142.5708 135.6768 1001.3808 14.5127 Se puede ver en las gráficas y tablas de resultados anteriores como existen asimetrías en la evolución de las magnitudes. Esto es lógico si se tiene en cuenta los cambios que aparecen en la dinámica del sistema cuando se introducen perturbaciones, y esto se acentúa aún más cuando se aumenta el número de tetheres. Para terminar, es interesante realizar un estudio del trabajo máximo desarrollado por tether (es decir max(W)/Nt ), teníendose en cuenta que esta es una magnitud que se ha calculado a nivel de sistema completo. Tabla 5.4 Trabajo unitario desarrollado por cada tether. NtPerturbación Coste unitario por tether[J] 2 36.7221 6βi21.6028 12 14.6605 2 6574.696 6ui6577.0463 12 66579.7012 2 22112.1833 6vi22100.4681 12 22105.9887 Se puede ver en la tabla 5.4 como el coste unitario en términos de trabajo se mantiene prácticamnete constante o incluso disminuye al aumentar el número de tetheres. Esto es interesante porque para un mismo tipo de perturbaciones, el coste del control no aumentaría si aumenta el número de tetheres, como en primera instancia podría resultar lógico.
6 Conclusiones finales y líneas futuras E n este último capítulo del documento, se van a reunir las conclusiones a las que se ha llegado en este proyecto y además se propondrán nuevas líneas de continuación o mejora que podrían surgir del trabajo que se ha realizado. 6.1 Conclusiones finales Recordando lo que se ha desarrollado a lo largo del documento, lo que primero que se hizo fue definir el modelo dinámico de E-sail que posteriormente sería el objeto de todo el estudio que se ha hecho. El modelo incluía un cuerpo rígido central y un número variable de tetheres en el que se tenía en cuenta deformación axial y a flexión. Posteriormente, el modelo se linealizó en torno a un equilibrio correspondiente a un giro con velocidad angular constante ω con la única deformación axial de los tetheres producida por la propia fuerza centrífuga. Una vez se tenía el modelo linealizado, se realizó un estudio modal para un número variable de tetheres de manera que los modos de representaron en el Capítulo 3. También en este capítulo, se obtuvieron los autovalores del sistema y se vió como el sistema era inestable, por lo que se decidió realizar un barrido de diferentes parámetros del sistema para ver si de esta manera se podía llegar a configuraciones estables, observando de igual modo que la inestabilidad seguía presente al haber siempre autovalores con parte real positiva. Tras esto, en el Capítulo 4, lo siguiente fue integrar el sistema de ecuaciones diferenciales no lineales que definían el modelo de dos formas: sin disipación y anádiendo disipación. En la integración lo que se hizo fue dar unas condiciones iniciales correspondientes al instante inicial de la trayectoria de referencia, de manera que al ser esta solución del sistema, el modelo debería seguirla durante todo el tiempo de ingración. Sin embargo, esto no currió pues como desveló el Capítulo 3, el sistema era inherentemente inestable y las propias perturbaciones numéricas asociadas a la integración sacaban al sistema de la trayectoria de referencia. A pesar de que las integraciones con disipación presentaban algo más de estabilidad, igualmente el sistema se desviaba de su equilibrio. Por último, con el objetivo de conseguir un modelo estable, se aplicaron los conceptos de control LQR al sistema. Se definió la actuación en el modelo y se introdujo en las ecuaciones de movimiento a través del término de las fuerzas generalizadas. Este nuevo sistema con actuación se volvió a linealizar en torno al movimiento de referencia y se obtuvo la ley de control correspondiente a este tipo de controladores. Esta ley se introdujo en el sistema no lineal y se repitieron las simulaciones numéricas obteniendo que el sistema ahora sí que seguía la trayectoria de referencia durante todo el periodo de integración. Además, se comprobó la robustez del control y se calcularon magnitudes asociadas al coste del control tales como energía, potencia y fuerzas implicadas. 85
92 Apéndice A. Variación de los autovalores según los parámetros 4.5 5 5.5 6 6.5 7 10-5 -0.08 -0.06 -0.04 -0.02 0 Parte real Parte imaginaria Figura A.29 Variación de λ11 con a. 4.5 5 5.5 6 6.5 7 10-5 0 0.02 0.04 0.06 0.08 Parte real Parte imaginaria Figura A.30 Variación de λ12 con a. 4.5 5 5.5 6 6.5 7 10-5 -0.08 -0.06 -0.04 -0.02 0 Parte real Parte imaginaria Figura A.31 Variación de λ13 con a. 4.5 5 5.5 6 6.5 7 10-5 0 0.02 0.04 0.06 0.08 Parte real Parte imaginaria Figura A.32 Variación de λ14 con a. 4.5 5 5.5 6 6.5 7 10-5 0 0.02 0.04 0.06 0.08 0.1 0.12 Parte real Parte imaginaria Figura A.33 Variación de λ15 con a. 4.5 5 5.5 6 6.5 7 10-5 -0.12 -0.1 -0.08 -0.06 -0.04 -0.02 0 Parte real Parte imaginaria Figura A.34 Variación de λ16 con a.
93 4.5 5 5.5 6 6.5 7 10-5 0 0.02 0.04 0.06 0.08 0.1 0.12 Parte real Parte imaginaria Figura A.35 Variación de λ17 con a. 4.5 5 5.5 6 6.5 7 10-5 -0.12 -0.1 -0.08 -0.06 -0.04 -0.02 0 Parte real Parte imaginaria Figura A.36 Variación de λ18 con a. Variación con b 2.2 2.4 2.6 2.8 3 3.2 10-4 -1 -0.5 0 0.5 1 Parte real Parte imaginaria Figura A.37 Variación de λ1con b. 2.2 2.4 2.6 2.8 3 3.2 10-4 -1 -0.5 0 0.5 1 Parte real Parte imaginaria Figura A.38 Variación de λ2con b. 2.2 2.4 2.6 2.8 3 3.2 10-4 0 0.2 0.4 0.6 0.8 1 10-4 Parte real Parte imaginaria Figura A.39 Variación de λ3con b. 2.2 2.4 2.6 2.8 3 3.2 10-4 -1 -0.8 -0.6 -0.4 -0.2 010-4 Parte real Parte imaginaria Figura A.40 Variación de λ4con b.
94 Apéndice A. Variación de los autovalores según los parámetros 2.2 2.4 2.6 2.8 3 3.2 10-4 0 0.5 1 1.5 2 2.5 3 3.5 10-3 Parte real Parte imaginaria Figura A.41 Variación de λ5con b. 2.2 2.4 2.6 2.8 3 3.2 10-4 -3.5 -3 -2.5 -2 -1.5 -1 -0.5 010-3 Parte real Parte imaginaria Figura A.42 Variación de λ6con b. 2.2 2.4 2.6 2.8 3 3.2 10-4 0 1 2 3 410-3 Parte real Parte imaginaria Figura A.43 Variación de λ7con b. 2.2 2.4 2.6 2.8 3 3.2 10-4 -4 -3 -2 -1 010-3 Parte real Parte imaginaria Figura A.44 Variación de λ8con b. 2.2 2.4 2.6 2.8 3 3.2 10-4 0 0.005 0.01 0.015 Parte real Parte imaginaria Figura A.45 Variación de λ9con b. 2.2 2.4 2.6 2.8 3 3.2 10-4 -0.015 -0.01 -0.005 0 Parte real Parte imaginaria Figura A.46 Variación de λ10 con b.
95 2.2 2.4 2.6 2.8 3 3.2 10-4 0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 Parte real Parte imaginaria Figura A.47 Variación de λ11 con b. 2.2 2.4 2.6 2.8 3 3.2 10-4 -0.07 -0.06 -0.05 -0.04 -0.03 -0.02 -0.01 0 Parte real Parte imaginaria Figura A.48 Variación de λ12 con b. 2.2 2.4 2.6 2.8 3 3.2 10-4 0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 Parte real Parte imaginaria Figura A.49 Variación de λ13 con b. 2.2 2.4 2.6 2.8 3 3.2 10-4 -0.07 -0.06 -0.05 -0.04 -0.03 -0.02 -0.01 0 Parte real Parte imaginaria Figura A.50 Variación de λ14 con b. 2.2 2.4 2.6 2.8 3 3.2 10-4 0 0.05 0.1 0.15 Parte real Parte imaginaria Figura A.51 Variación de λ15 con b. 2.2 2.4 2.6 2.8 3 3.2 10-4 -0.15 -0.1 -0.05 0 Parte real Parte imaginaria Figura A.52 Variación de λ16 con b.
96 Apéndice A. Variación de los autovalores según los parámetros 2.2 2.4 2.6 2.8 3 3.2 10-4 0 0.05 0.1 0.15 Parte real Parte imaginaria Figura A.53 Variación de λ17 con b. 2.2 2.4 2.6 2.8 3 3.2 10-4 -0.15 -0.1 -0.05 0 Parte real Parte imaginaria Figura A.54 Variación de λ18 con b. Variación con b 2.2 2.4 2.6 2.8 3 3.2 10-4 -1 -0.5 0 0.5 1 Parte real Parte imaginaria Figura A.55 Variación de λ1con b. 2.2 2.4 2.6 2.8 3 3.2 10-4 -1 -0.5 0 0.5 1 Parte real Parte imaginaria Figura A.56 Variación de λ2con b. 2.2 2.4 2.6 2.8 3 3.2 10-4 0 0.2 0.4 0.6 0.8 1 10-4 Parte real Parte imaginaria Figura A.57 Variación de λ3con b. 2.2 2.4 2.6 2.8 3 3.2 10-4 -1 -0.8 -0.6 -0.4 -0.2 010-4 Parte real Parte imaginaria Figura A.58 Variación de λ4con b.
97 2.2 2.4 2.6 2.8 3 3.2 10-4 0 0.5 1 1.5 2 2.5 3 3.5 10-3 Parte real Parte imaginaria Figura A.59 Variación de λ5con b. 2.2 2.4 2.6 2.8 3 3.2 10-4 -3.5 -3 -2.5 -2 -1.5 -1 -0.5 010-3 Parte real Parte imaginaria Figura A.60 Variación de λ6con b. 2.2 2.4 2.6 2.8 3 3.2 10-4 0 1 2 3 410-3 Parte real Parte imaginaria Figura A.61 Variación de λ7con b. 2.2 2.4 2.6 2.8 3 3.2 10-4 -4 -3 -2 -1 010-3 Parte real Parte imaginaria Figura A.62 Variación de λ8con b. 2.2 2.4 2.6 2.8 3 3.2 10-4 0 0.005 0.01 0.015 Parte real Parte imaginaria Figura A.63 Variación de λ9con b. 2.2 2.4 2.6 2.8 3 3.2 10-4 -0.015 -0.01 -0.005 0 Parte real Parte imaginaria Figura A.64 Variación de λ10 con b.
98 Apéndice A. Variación de los autovalores según los parámetros 2.2 2.4 2.6 2.8 3 3.2 10-4 0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 Parte real Parte imaginaria Figura A.65 Variación de λ11 con b. 2.2 2.4 2.6 2.8 3 3.2 10-4 -0.07 -0.06 -0.05 -0.04 -0.03 -0.02 -0.01 0 Parte real Parte imaginaria Figura A.66 Variación de λ12 con b. 2.2 2.4 2.6 2.8 3 3.2 10-4 0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 Parte real Parte imaginaria Figura A.67 Variación de λ13 con b. 2.2 2.4 2.6 2.8 3 3.2 10-4 -0.07 -0.06 -0.05 -0.04 -0.03 -0.02 -0.01 0 Parte real Parte imaginaria Figura A.68 Variación de λ14 con b. 2.2 2.4 2.6 2.8 3 3.2 10-4 0 0.05 0.1 0.15 Parte real Parte imaginaria Figura A.69 Variación de λ15 con b. 2.2 2.4 2.6 2.8 3 3.2 10-4 -0.15 -0.1 -0.05 0 Parte real Parte imaginaria Figura A.70 Variación de λ16 con b.
99 2.2 2.4 2.6 2.8 3 3.2 10-4 0 0.05 0.1 0.15 Parte real Parte imaginaria Figura A.71 Variación de λ17 con b. 2.2 2.4 2.6 2.8 3 3.2 10-4 -0.15 -0.1 -0.05 0 Parte real Parte imaginaria Figura A.72 Variación de λ18 con b. Variación con e 0.3 0.35 0.4 -1 -0.5 0 0.5 1 Parte real Parte imaginaria Figura A.73 Variación de λ1con e. 0.3 0.35 0.4 -1 -0.5 0 0.5 1 Parte real Parte imaginaria Figura A.74 Variación de λ2con e. 0.3 0.35 0.4 0 0.2 0.4 0.6 0.8 1 10-4 Parte real Parte imaginaria Figura A.75 Variación de λ3con e. 0.3 0.35 0.4 -1 -0.8 -0.6 -0.4 -0.2 010-4 Parte real Parte imaginaria Figura A.76 Variación de λ4con e.
100 Apéndice A. Variación de los autovalores según los parámetros 0.3 0.35 0.4 0 0.5 1 1.5 2 2.5 3 3.5 10-3 Parte real Parte imaginaria Figura A.77 Variación de λ5con e. 0.3 0.35 0.4 -3.5 -3 -2.5 -2 -1.5 -1 -0.5 010-3 Parte real Parte imaginaria Figura A.78 Variación de λ6con e. 0.3 0.35 0.4 0 1 2 3 410-3 Parte real Parte imaginaria Figura A.79 Variación de λ7con e. 0.3 0.35 0.4 -4 -3 -2 -1 010-3 Parte real Parte imaginaria Figura A.80 Variación de λ8con e. 0.3 0.35 0.4 0 0.002 0.004 0.006 0.008 0.01 0.012 Parte real Parte imaginaria Figura A.81 Variación de λ9con e. 0.3 0.35 0.4 -0.012 -0.01 -0.008 -0.006 -0.004 -0.002 0 Parte real Parte imaginaria Figura A.82 Variación de λ10 con e.
101 0.3 0.35 0.4 0 0.01 0.02 0.03 0.04 0.05 0.06 Parte real Parte imaginaria Figura A.83 Variación de λ11 con e. 0.3 0.35 0.4 -0.06 -0.05 -0.04 -0.03 -0.02 -0.01 0 Parte real Parte imaginaria Figura A.84 Variación de λ12 con e. 0.3 0.35 0.4 -0.06 -0.05 -0.04 -0.03 -0.02 -0.01 0 Parte real Parte imaginaria Figura A.85 Variación de λ13 con e. 0.3 0.35 0.4 0 0.01 0.02 0.03 0.04 0.05 0.06 Parte real Parte imaginaria Figura A.86 Variación de λ14 con e. 0.3 0.35 0.4 0 0.05 0.1 0.15 Parte real Parte imaginaria Figura A.87 Variación de λ15 con e. 0.3 0.35 0.4 -0.15 -0.1 -0.05 0 Parte real Parte imaginaria Figura A.88 Variación de λ16 con e.
108 Apéndice B. Respuestas con control 0 20 40 60 80 -30 -25 -20 -15 -10 -5 0 0 20 40 60 80 -10 -5 0 5 10 0 20 40 60 80 0 0.5 1 1.5 2105 0 20 40 60 80 -30 -25 -20 -15 -10 -5 0 0 20 40 60 80 -5 0 5 10 0 20 40 60 80 0 0.5 1 1.5 2105 0 20 40 60 80 -30 -25 -20 -15 -10 -5 0 5 0 20 40 60 80 -10 -5 0 5 10 0 20 40 60 80 -2 -1.5 -1 -0.5 0105 0 20 40 60 80 -30 -25 -20 -15 -10 -5 0 0 20 40 60 80 -5 0 5 10 0 20 40 60 80 0 0.5 1 1.5 2105 0 20 40 60 80 -30 -25 -20 -15 -10 -5 0 0 20 40 60 80 -5 0 5 10 0 20 40 60 80 0 0.5 1 1.5 2105 0 20 40 60 80 -30 -25 -20 -15 -10 -5 0 0 20 40 60 80 -5 0 5 10 0 20 40 60 80 0 0.5 1 1.5 2105
109 0 20 40 60 80 -30 -25 -20 -15 -10 -5 0 0 20 40 60 80 -10 -5 0 5 10 0 20 40 60 80 0 0.5 1 1.5 2105 0 20 40 60 80 -30 -25 -20 -15 -10 -5 0 5 0 20 40 60 80 -5 0 5 10 0 20 40 60 80 -2 -1.5 -1 -0.5 0105 0 20 40 60 80 -30 -25 -20 -15 -10 -5 0 0 20 40 60 80 -5 0 5 10 0 20 40 60 80 0 0.5 1 1.5 2105 Figura B.1 Respuestas para Nt=12 con control.
Índice de Figuras 1.1 Arquitectura de la E-sail 2 2.1 Modelo dinámico de la E-sail 8 2.2 Sistemas de referencia 10 3.1 Modos del 1 al 3 22 3.2 Modos del 4 al 6 22 3.3 Modos del 7 al 9 23 3.4 Modos del 1 al 3 23 3.5 Modos del 4 al 6 23 3.6 Modos del 7 al 9 24 3.7 Modos del 10 al 12 24 3.8 Modos del 13 al 15 24 3.9 Modos del 16 al 18 25 3.10 Modos del 19 al 21 25 3.11 Modos del 1 al 3 26 3.12 Modos del 4 al 6 26 3.13 Modos del 7 al 9 26 3.14 Modos del 10 al 12 27 3.15 Modos del 13 al 15 27 3.16 Modos del 16 al 18 27 3.17 Modos del 19 al 21 28 3.18 Modos del 22 al 24 28 3.19 Modos del 25 al 27 28 3.20 Modos del 28 al 30 29 3.21 Modos del 31 al 33 29 3.22 Modos del 34 al 36 29 3.23 Modos del 37 al 39 30 3.24 Modos del 40 al 42 30 3.25 Modos del 43 al 44 30 3.26 Autovalores del sistema para Nt=231 3.27 Variación de los autovalores con ωpara Nt=232 3.28 Variación de los autovalores con apara Nt=233 3.29 Variación de los autovalores con bpara Nt=233 3.30 Variación de los autovalores con epara Nt=234 3.31 Variación de los autovalores con Lpara Nt=234 111
112 Índice de Figuras 3.32 Autovalores del sistema para Nt=635 3.33 Variación de los autovalores con ωpara Nt=635 3.34 Variación de los autovalores con apara Nt=636 3.35 Variación de los autovalores con bpara Nt=636 3.36 Variación de los autovalores con epara Nt=637 3.37 Variación de los autovalores con Lpara Nt=637 3.38 Autovalores del sistema para Nt=12 38 3.39 Variación de los autovalores con ωpara Nt=12 38 3.40 Variación de los autovalores con apara Nt=12 39 3.41 Variación de los autovalores con bpara Nt=12 39 3.42 Variación de los autovalores con epara Nt=12 40 3.43 Variación de los autovalores con Lpara Nt=640 4.1 Respuesta de la coordenada x42 4.2 Respuesta de la coordenada y43 4.3 Respuesta de la coordenada θ43 4.4 Respuesta de la coordenada βA44 4.5 Respuesta de la coordenada uA44 4.6 Respuesta de la coordenada vA44 4.7 Respuesta de la coordenada βB45 4.8 Respuesta de la coordenada uB45 4.9 Respuesta de la coordenada βB45 4.10 Respuesta de las coordenadas x,y46 4.11 Respuesta de la coordenada θ46 4.12 Respuestas de las coordenadas βi46 4.13 Respuestas de la coordenadas ui46 4.14 Respuesta de la coordenada vi47 4.15 Respuesta de las coordenadas x,y47 4.16 Respuesta de la coordenada θ47 4.17 Respuesta de las coordenadas βi48 4.18 Respuesta de las coordenadas ui48 4.19 Respuesta de las coordenadas vi48 4.20 Respuesta de la coordenada x51 4.21 Respuesta de la coordenada y51 4.22 Respuesta de la coordenada θ51 4.23 Respuesta de la coordenada βA52 4.24 Respuesta de la coordenada uA52 4.25 Respuesta de la coordenada vA53 4.26 Respuesta de la coordenada βB53 4.27 Respuesta de la coordenada uB53 4.28 Respuesta de la coordenada vB53 4.29 Respuesta de la coordenada vB55 4.30 Respuestas para Nt=657 5.1 Actuación del sistema 61 5.2 Fuerza aplicada en los tetheres 63 5.3 Diagrama de flujo 66 5.4 Autovalores con controlador 68 5.5 Respuesta de x68 5.6 Respuesta de y69
Índice de Figuras 113 5.7 Respuesta de θ69 5.8 Respuesta de βA69 5.9 Respuesta de uA70 5.10 Respuesta de vA70 5.11 Respuesta de βB70 5.12 Respuesta de uB70 5.13 Respuesta de vB71 5.14 Autovalores con controlador 71 5.15 Respuestas para Nt=6con control 73 5.16 Autovalores con control 73 5.17 Respuesta de x74 5.18 Respuesta de y74 5.19 Respuesta de θ74 5.20 Respuesta de βA74 5.21 Respuesta de uA74 5.22 Respuesta de vA74 5.23 Respuesta de x76 5.24 Respuesta de y76 5.25 Respuesta de θ76 5.26 Respuesta de βA76 5.27 Respuesta de uA76 5.28 Respuesta de vA76 5.29 Respuesta de x77 5.30 Respuesta de y77 5.31 Respuesta de θ77 5.32 Respuesta de βA77 5.33 Respuesta de uA78 5.34 Respuesta de vA78 5.35 Respuesta de x78 5.36 Respuesta de y78 5.37 Respuesta de θ79 5.38 Respuesta de βA79 5.39 Respuesta de uA79 5.40 Respuesta de vA79 5.41 Trabajo desarrollado según perturbaciones en βi80 5.42 Potencia máxima desarrollada según perturbaciones en βi80 5.43 Fuerza máxima desarrollada según perturbaciones en βi81 5.44 Trabajo desarrollado según perturbaciones en ui81 5.45 Potencia máxima desarrollada según perturbaciones en ui82 5.46 Fuerza máxima desarrollada según perturbaciones en ui82 5.47 Trabajo desarrollado según perturbaciones en ui83 5.48 Potencia máxima desarrollada según perturbaciones en vi83 5.49 Fuerza máxima desarrollada según perturbaciones en vi83 A.1 Variación de λ1con ω87 A.2 Variación de λ2con ω87 A.3 Variación de λ3con ω87 A.4 Variación de λ4con ω87 A.5 Variación de λ5con ω88 A.6 Variación de λ6con ω88
114 Índice de Figuras A.7 Variación de λ7con ω88 A.8 Variación de λ8con ω88 A.9 Variación de λ9con ω88 A.10 Variación de λ10 con ω88 A.11 Variación de λ11 con ω89 A.12 Variación de λ12 con ω89 A.13 Variación de λ13 con ω89 A.14 Variación de λ14 con ω89 A.15 Variación de λ15 con ω89 A.16 Variación de λ16 con ω89 A.17 Variación de λ17 con ω90 A.18 Variación de λ18 con ω90 A.19 Variación de λ1con a90 A.20 Variación de λ2con a90 A.21 Variación de λ3con a90 A.22 Variación de λ4con a90 A.23 Variación de λ5con a91 A.24 Variación de λ6con a91 A.25 Variación de λ7con a91 A.26 Variación de λ8con a91 A.27 Variación de λ9con a91 A.28 Variación de λ10 con a91 A.29 Variación de λ11 con a92 A.30 Variación de λ12 con a92 A.31 Variación de λ13 con a92 A.32 Variación de λ14 con a92 A.33 Variación de λ15 con a92 A.34 Variación de λ16 con a92 A.35 Variación de λ17 con a93 A.36 Variación de λ18 con a93 A.37 Variación de λ1con b93 A.38 Variación de λ2con b93 A.39 Variación de λ3con b93 A.40 Variación de λ4con b93 A.41 Variación de λ5con b94 A.42 Variación de λ6con b94 A.43 Variación de λ7con b94 A.44 Variación de λ8con b94 A.45 Variación de λ9con b94 A.46 Variación de λ10 con b94 A.47 Variación de λ11 con b95 A.48 Variación de λ12 con b95 A.49 Variación de λ13 con b95 A.50 Variación de λ14 con b95 A.51 Variación de λ15 con b95 A.52 Variación de λ16 con b95 A.53 Variación de λ17 con b96 A.54 Variación de λ18 con b96 A.55 Variación de λ1con b96 A.56 Variación de λ2con b96
Índice de Figuras 115 A.57 Variación de λ3con b96 A.58 Variación de λ4con b96 A.59 Variación de λ5con b97 A.60 Variación de λ6con b97 A.61 Variación de λ7con b97 A.62 Variación de λ8con b97 A.63 Variación de λ9con b97 A.64 Variación de λ10 con b97 A.65 Variación de λ11 con b98 A.66 Variación de λ12 con b98 A.67 Variación de λ13 con b98 A.68 Variación de λ14 con b98 A.69 Variación de λ15 con b98 A.70 Variación de λ16 con b98 A.71 Variación de λ17 con b99 A.72 Variación de λ18 con b99 A.73 Variación de λ1con e99 A.74 Variación de λ2con e99 A.75 Variación de λ3con e99 A.76 Variación de λ4con e99 A.77 Variación de λ5con e100 A.78 Variación de λ6con e100 A.79 Variación de λ7con e100 A.80 Variación de λ8con e100 A.81 Variación de λ9con e100 A.82 Variación de λ10 con e100 A.83 Variación de λ11 con e101 A.84 Variación de λ12 con e101 A.85 Variación de λ13 con e101 A.86 Variación de λ14 con e101 A.87 Variación de λ15 con e101 A.88 Variación de λ16 con e101 A.89 Variación de λ17 con e102 A.90 Variación de λ18 con e102 A.91 Variación de λ1con L102 A.92 Variación de λ2con L102 A.93 Variación de λ3con L102 A.94 Variación de λ4con L102 A.95 Variación de λ5con L103 A.96 Variación de λ6con L103 A.97 Variación de λ7con L103 A.98 Variación de λ8con L103 A.99 Variación de λ9con L103 A.100 Variación de λ10 con L103 A.101 Variación de λ11 con L104 A.102 Variación de λ12 con L104 A.103 Variación de λ13 con L104 A.104 Variación de λ14 con L104 A.105 Variación de λ15 con L104 A.106 Variación de λ16 con L104
116 Índice de Figuras A.107 Variación de λ17 con L105 A.108 Variación de λ18 con L105 B.1 Respuestas para Nt=12 con control 109
Índice de Tablas 3.1 Valores numéricos de las magnitudes 21 3.2 Valores numéricos de los parámetros 21 4.1 Valores del amortiguamiento 50 5.1 Valores de las magnitudes según perturbaciones en βi81 5.2 Valores de las magnitudes según perturbaciones en ui82 5.3 Valores de las magnitudes según perturbaciones en vi84 5.4 Trabajo unitario desarrollado por cada tether 84 117