VEHÍCULOS ESPACIALES Y MISILES Práctica 4 de STK: Representación de la actitud mediante STK. Estudio de la dinámica de actitud mediante MATLAB/STK. NOMBRE DEL ALUMNO:

En esta práctica se emplearán los conceptos de representación de actitud aprendidos en la asignatura y se utilizará el interfaz Matlab-STK para demostrar algunas aplicaciones avanzadas. Se recomienda consultar la práctica 2 (disponible en la página web) si no se recuerda como interconectar Matlab y STK. Evaluación: Se pide rellenar el boletín de prácticas con ciertos resultados analíticos u obtenidos de STK, realizando los pasos indicados en la práctica.

1. Representación de actitud en STK A continuación vamos a aprender las opciones básicas de representación de actitud de STK. En primer lugar creamos un escenario nuevo con las opciones por defecto y en él creamos un satélite llamado “Sat” en órbita de semieje mayor 8800 km con excentricidad 0.2, y el resto de sus elementos orbitales (que no son muy relevantes a efectos de esta demostración) a elección del alumno. Pinchamos en el icono “View From/To” de la ventana 3D, y pinchamos nuestro satélite en “View from” y en “View to”. Cambiar la vista para ver desde el sistema de referencia inercial, pinchando en “Reference frame” la opción “Earth Inertial Axes”. Pulsamos OK. Hacer zoom hacia el satélite y encuadrarlo en la pantalla. Observemos que por defecto tiene una actitud que permanece prácticamente constante relativa a la órbita (girando con ella) y sus paneles solares siempre apuntan al sol. Esto se debe al perfil de actitud por defecto, que estudiamos a continuación. Para estudiar adecuadamente la actitud, primero añadimos una representación de los ejes cuerpo, ejes de la órbita (VVLH, Vehicle Velocity Local Horizontal) y del vector del Sol: para ello hacemos doble click en el satélite para obtener sus propiedades y pinchamos en “3D->Vector”. Pinchamos “Show” en Body Axes , en VVLH y en “Sun Vector”. Pinchamos OK y simulamos. Animar la simulación. Observar en primer lugar el comportamiento de los ejes VVLH: el eje Z apunta a la tierrra, el eje X está contenido en el plano vertical de la velocidad del vehículo (coincidirá con ella en una órbita circular, pero no en una órbita elíptica) y el eje Y es normal a la órbita. Observemos que los ejes cuerpos permanecen con actitud casi constante respecto a estos ejes (sería totalmente constante si no hubiera excentricidad). Observamos también que los paneles solares siguen constantemente al vector del Sol. Para cambiar el comportamiento por defecto, pinchamos “Attitude” en las propiedades del satélite.

El desplegable de la parte superior de la ventana nos permite especificar la actitud de tres maneras, Standard, Real-Time y Multi Segment. 1. Pinchando en la opción Standard, aparecen 3 tipos de perfil de actitud: • Basic: Perfil predefinido por STK (no obedece necesariamente las ecuaciones de la dinámica de actitud). El tipo (type) que viene por defecto, Nadir alignment with ECF velocity constraint, alinea el satélite con la órbita, manteniendo el eje Z cuerpo apuntando hacia el nadir (eje Z de VVLH). De ahí viene “... nadir constraint”. La otra dirección la determina la velocidad, que en una órbita excéntrica no es constante respecto a la órbita, de ahí la lenta variación del sistema de referencia. Probemos a cambiarlo seleccionado en “type”, el perfil “Sun alignment with nadir constraint”. Pinchamos Ok para cerrar la ventana de propiedades y observemos como cambia la actitud del vehículo. NOTA: El hecho de que los paneles solares apunten al Sol viene dado por la propiedad “Model pointing” en 3D Graphics. Otros modelos permiten apuntar antenas, etc.. a diferentes blancos. Volvamos a la pantalla de propiedades de actitud. • Target Pointing: Durante los intervalos de tiempo seleccionados se sustituye al perfil básico con un perfil de apuntamiento al objetivo seleccionado. Para probarlo, creemos una base en Baikonur. Pinchamos “Attitude” en las propiedades del satélite, marcamos “Override Basic attitude for selected targets” y pinchamos “Select targets” en el cuadro “Target Pointing”. Pinchamos la base y la desplazamos a “Assigned Targets” con la flecha. Pulsamos OK. Podemos insertar un vector que vaya desde el satélite hacia la base: pinchamos en “3D>Vector” en las propiedades del satélite. Pinchamos “Add”, seleccionamos Satellite/To Vectors/Tyuratam_Baikonur_CIS_LC en el árbol “Available”, y lo pasamos a “Selected” pinchando la flecha. Pulsamos OK. Seleccionamos el vector Tyuratam_Baikonur_CIS_LC y le cambiamos el color, a blanco (por ejemplo). Pulsamos OK para cerrar la ventana de propiedades y simulamos. Observar que el perfil de actitud de apuntamiento a la base sólo prevalece sobre el perfil básico de apuntamiento al Sol cuando Baikonur es accesible desde la posición del satélite; entre ambos perfiles hay un cambio casi instantáneo. • Precomputado: Proporcionamos la información de actitud en el tiempo a través de un archivo de actitud. Este archivo contiene un vector de tiempos y una matriz de cuaterniones que da la actitud para cada instante de tiempo del vector de tiempos. No se utilizará esta opción directamente en esta práctica, sino que se impondrá a través del interfaz MatlabSTK. 2. La opción Real-Time permite hacer cálculos de actitud en tiempo real intercambiando información con Matlab a través del interfaz Matlab-STK. 3. La opción Multi-Segment permite ordenar en el tiempo diferentes perfiles de actitud. Cuando pinchamos más de un tipo de perfil de actitud bajo la opción Standard, STK genera automáticamente un perfil de actitud multi segmento. Pinchamos “Attitude” en las propiedades del satélite y elegimos “Multi Segment”. Observamos cómo STK ha generado un perfil multisegmento seleccionando Baikonur sólo cuando Baikonur es accesible para el satélite. Como se ha podido comprobar, las transiciones entre perfiles de actitud son instantáneas. Pueden sin embargo crearse perfiles de actitud que simulen una transición más realista teniendo en cuenta la dinámica del vehículo. Estos perfiles de transición se insertan entre los perfiles de actitud mediante perfiles multisegmento, y se generan mediante el Attitude Simulator, cuyo uso escapa del alcance de esta práctica.

2. Añadir el telescopio espacial Hubble. Fijar como objetivo la estrella Sirio (sirius) que habrá que añadir y en 3D graphics-model pointing, fijar que las antenas (HGA) apunten a Canarias (crear la base, buscando por Maspalomas, el observatorio canario). Añadir los vectores adecuados para verificar que esto es cierto (para apuntar a Sirio crear un vector de type “to star”). Mostrar la simulación al profesor.

2. Ejemplos: Representación de actitud a través de Matlab/STK Vamos a representar en STK distintos perfiles de actitud, generados en Matlab. Recordemos que para establecer la interconexión Matlab-STK hay que ejecutar los archivos de inicialización tal como se describía en la anterior práctica. Lo haremos a través del comando StkSetAttitudeCBI: stkSetAttitudeCBI(objPath, cb, time, quats, aFilePath): Manda a STK información de actitud en el tiempo, respecto al sistema de referencia inercial J2000. • ObjPath: Cadena de texto con la ruta del satélite al que se aplica. • Cb: Cadena de texto con la el nombre del cuerpo central. Normalmente será ‘Earth’ • Time: Vector 1xN de tiempos en los que se aplica cada dato de actitud. Medido en segundos desde la época. • quats: Matriz de 4xN, la columna i-ésima incluye el cuaternión asociado al instante dado por el elemento i-ésimo del vector Time. OJO: Tal como STK define cuaterniones, la cuarta componente es la parte escalar y las tres primeras la parte vectorial! • aFilePath: Nombre del archivo de actitud generado. Si no se especifica el argumento, le asignará un nombre automático. Cuando ejecutamos la orden stkSetAttitudeCBI, ésta genera un archivo de actitud y lo asigna al satélite como Perfil de Actitud Precomputado para que lo siga durante los tiempos determinados en el vector de tiempos. Ejemplo: stkSetAttitudeCBI('*/Satellite/Sat', 'Earth', T, quat) Desde STK, borramos las bases y el satélite “Sat” y creamos un nuevo satélite llamado una vez más “Sat”. Los elementos orbitales del satélite pueden escogerse como quiera el alumno, pues no serán importantes en el resto de la práctica. Ahora ajustamos las propiedades gráficas para ver la actitud del satélite. Pinchamos en el icono “View From/To” de la ventana 3D, y pinchamos nuestro satélite en “View from” y en “View to”. Cambiar la vista para ver desde el sistema de referencia inercial, pinchando en “Reference frame” la opción “Earth Inertial Axes”. Pulsamos OK. Para ver la traza que siguen los ejes cuerpo respecto a los ejes inerciales J2000, pinchamos en “3D->Vector” en las propiedades del satélite. Pinchamos “Show” en los sistemas “J2000” y “Body Axes”. Seleccionamos “Body Axes” y en el cuadro “Persistence”, pinchamos “Show” y seleccionamos “Trace” en el desplegable “Connect”. Pulsamos OK para cerrar la ventana de propiedades. También es útil poder vigilar los cuaterniones y otros parámetros en tiempo real. Para ello en 3D graphics/Data Display se puede seleccionar “Attitude Quaternions” (también Euler Angles).

Ejemplo 1. En primer lugar, vamos a generar un perfil de actitud constante en el tiempo, respecto al sistema inercial J2000. Para ello, dada la actitud en ángulos de Euler, la transformamos a cuaternión con la función “AtbEulerToQuat”, generamos el vector de tiempos y la matriz de actitud en función del tiempo, incluyendo en todas las columnas el cuaternión asociado a la actitud constante. Escribimos el siguiente script (en la pág. web se puede encontrar como Ej1Prac3.m): euler_const=[10 45 120]'*pi/180; quat_const=atbEulerToQuat(euler_const, 313); %secuencia 3-1-3 tiempos=[0:1:20000]; quat=zeros(4,length(tiempos)); for i=1:length(tiempos) quat(:,i)=quat_const; end stkSetAttitudeCBI('*/Satellite/Sat', 'Earth', tiempos, quat);

Ejecutar el script de Matlab, y pulsar reset y play en STK. Visualizar el resultado. Ojo: Verificar que la animación se visualiza en el periodo de tiempo correcto (desde el inicio de la práctica hasta pasados 20000 segundos). Una forma de ver que se ha cargado el perfil de actitud es en las propiedades del satélite, Basic/Attitude, aparecerá marcado precomputed por una serie de tiempos, con un archivo .a que debe haberse creado conteniendo el perfil de actitud para cada instante. Ejemplo 2. Ahora creamos un perfil de actitud con velocidad angular constante. Para ello, a partir de una actitud inicial, integramos con ode45 las ecuaciones diferenciales cinemáticas en cuaterniones. Escribimos el siguiente script (en la pág. web se puede encontrar como Ej2Prac3.m): function Ej2Prac4(); euler_ini=[10 90 -45]'*pi/180; q_ini=atbEulerToQuat(euler_ini, 313); tiempos=[0:1:20000]; [tiempos,quat]=ode45(@fun_cinematica, tiempos, q_ini); stkSetAttitudeCBI('*/Satellite/Sat', 'Earth', tiempos, quat') %-------------------------------------------------------------------------function q_dot=fun_cinematica(t,q); omega1=0.002; omega2=0.003; omega3=0.001; MatrizCine=0.5*[

0 -omega3 omega2 -omega1

q_dot=MatrizCine*q;

omega3 0 -omega1 -omega2

-omega2 omega1 0 -omega3

omega1 omega2 omega3 0

; ... ; ... ; ... ];

Observar que utilizando el argumento “t” en la función de las ecuaciones cinemáticas pueden escribirse el vector velocidad angular como función del tiempo. Ejecutar el programa de Matlab, y pulsar reset y play en STK. Ejemplo 3. Creamos ahora un perfil de actitud de maniobra. La maniobra partirá de unas actitudes inicial y final dadas, en un tiempo de maniobra dado. Para ello, usamos la siguiente propiedad: si partimos de un cuaternión inicial q0 a otro final, qF , existe un cuaternión que representa la rotación de velocidad angular mínima qrot que cumple: q F = qrot q0 Por tanto: * q0 1 * qrot = q F = qF = q F q0 2 q0 q0 Estas multiplicaciones de cuaterniones pueden hacerse fácilmente en Matlab mediante el comando AtbQuatXQuat del toolbox “AeroToolbox”, de manera que qrot = q F q0 * se escribe en Matlab *

qrot = atbQuatXquat (qF , q0 ) . Una vez obtenemos el cuaternión de rotación, extraemos sus eje y ángulo de Euler, y construimos el cuaternión de rotación instantánea: ⎡→ ⎛ ωt ⎞⎤ ⎢ e sin ⎜ 2 ⎟⎥ ⎝ ⎠⎥ q rot (t ) = ⎢ ⎛ ⎢ cos⎜ ωt ⎞⎟ ⎥ ⎢⎣ ⎝ 2 ⎠ ⎥⎦ donde la velocidad angular viene dada por:

θ

, tMAN siendo θ el ángulo de Euler y tMAN el tiempo de maniobra dado.

ω=

De esta manera, la actitud en función del tiempo es: q F (t ) = qrot (t )q0 En este ejemplo programamos una maniobra de una actitud inicial a otra final, dada en ángulos de Euler, que se realiza en un tiempo de maniobra tMAN, y permanece en la actitud final otro segmento de tiempo tMAN. Las actitudes inicial y final se pasan a cuaterniones mediante el comando “atbEulerToQuat” del toolbox AeroToolBox, al cual hay que dar como argumento la secuencia elegida para los ángulos de Euler. Escribimos el siguiente script (en la pág. web se puede encontrar como Ej3Prac3.m): dtr = pi/180; % conversion de grados a rad % Datos de partida euler_0=[-45 0 180]'*dtr; euler_F=[30 30 45]'*dtr; tMAN=3000;

% Vector de tiempos de duracion 2*tMAN, con 2 puntos por segundo tiempos=linspace(0,2*tMAN,tMAN*4); % Traducimos actitudes inicial y final a cuaterniones q0=atbEulerToQuat(euler_0, 313); qF=atbEulerToQuat(euler_F, 313); conj_q0=[-q0(1) -q0(2) -q0(3) q0(4)]'; qrot_F=atbQuatXquat(qF, conj_q0); % Obtenemos eje y angulo de Euler theta=2*acos(qrot_F(4)); e=[qrot_F(1) qrot_F(2) qrot_F(3)]'/sin(theta/2); omega=theta/tMAN; quat=zeros(4,length(tiempos)); % Escribimos informacion de actitud en la matriz "quat" for i=1:length(tiempos) t=tiempos(i); if tclose scenario) y ejecutar el script “iniciaSTK” para inicializar las simulaciones. Se genera el satélite que va a ser usado en todas las simulaciones. Como antes visualizar el satélite directamente en la ventana 3-D, en el sistema de referencia inercial, y mostrar los ejes J2000 y ejes cuerpo. Ejecutar el script “cuerpoLibre1.m”. Este script genera simulaciones para una rotación en torno al eje menor, intermedio y mayor (las tres primeras gráficas), así como una rotación en torno al eje mayor con precesión (las otras gráficas). Se utiliza un modelo sin disipación de energía. Las gráficas no dan mucha información, así que para visualizarlas usar el comando stkSetAttitudeCBI('*/Satellite/Pruebas_Actitud', 'Earth', times, quat); sustituyendo “quat” por qMay', qMen', y qInt' respectivamente (hay que trasponer

los vectores para que tengan la dimensión correcta), para ver los tres casos. La simulación sólo dura 200 segundos por lo que hay que ralentizar considerablemente la velocidad de simulación. NOTA: Si diera un error diciendo que no puede escribir en cierto archivo, avisar al profesor. 4. ¿Qué rotaciones son estables, qué rotaciones son inestables y por qué?

Ver ahora el caso de rotación con precesión, representado por qPrec stkSetAttitudeCBI('*/Satellite/Pruebas_Actitud', 'Earth', times, qPrec');

y para ver los “conos” que representan el movimiento, hacer lo siguiente: -Añadir la visualización del vector velocidad angular (en 3-D Graphics, Vector, Add y buscar AngVelocity Vector), seleccionar persistencia, y connect:trace. -Volver a Vecotr, Add, y elegir “Launch Vector Geometry Tool”. Seleccionar AngVelocity y pinchar el icono de “copiar” (el primero). Cerrar la “Vector Geometry Tool”. Seleccionar el nuevo vector creado AngVelocity1. Como antes seleccionar persistencia y connect:trace, pero ahora seleccionar también Axes: J2000. Si se han llevado a cabo estos pasos correctamente se visualizará un circulo deslizando sobre otro, son los vistos en teoría. 5. ¿Es este el caso prolato u oblato? Modificar el tensor de inercia en el script para observar el otro caso.

Cerrar las gráficas de MATLAB y borrar los vectores de velocidad angular. Ejecutar ahora el script “cuerpoLibre2.m”. Este script ejecuta un modelo con disipación de energía, tal como se vio en teoría (con una burbuja de combustible en el depósito que disipa energía por rozamiento). La figura 1 muestra la estabilidad (asintótica) del eje mayor (x), frente a una pequeña perturbación.

La figura 2 muestra la inestabilidad del eje menor (z). La figura 3 es del mismo caso y muestra como se disipa lentamente la energía (mientras que el emomento cinético permanece constante). La figura 4 muestra como las curvas polodia se vuelven convergentes al eje mayor, dibujando la elipse del momento cinético en el proceso. La figura 5 la inestabilidad del eje intermedio (y). En todos los casos el vehículo termina rotando en torno al eje mayor, aunque es imprevisible, si se comienza con una rotación intermedia o menor, si el resultado final será una rotación “positiva” o “negativa”. Esto se ejemplifica con la figura 6, con unas condiciones muy muy parecidas a las de la figura 2 pero que termina en una rotación positiva (la figura 2 termina en una rotación negativa). Ésto es un ejemplo de caos. Para ver estas simulaciones en STK se usa de nuevo el comando stkSetAttitudeCBI('*/Satellite/Pruebas_Actitud', 'Earth', times, quat);

cambiando quat por qMay', qMen', qInt', qMen2' respectivamente. Se puede visualizar la curva polodia (para los casos de ejes menor e intermedio) dibujando el vector velocidad angular y añadiendo persistencia en la traza, como antes. Aunque dichas curvas estarán dibujadas en ejes cuerpo y por tanto “se moverán”. Correr ahora el script cuerpoForz.m. Simula el satélite en rotación sometido a un par de perturbación. Se muestran tres figuras, la primera son las velocidades angulares, la segunda los ángulos evolucionando respecto al tiempo, y la tercera el ángulo uno frente al dos. Si la velocidad angular de rotación es suficientemente alta, el efecto del par perturbador es pequeño y la curva de los ángulos es un epiciclo de pequeña amplitud (si no es suficientemente alta, hay efectos no lineales y la curva es más compleja y de gran amplitud). La forma exacta del epiciclo depende de la relación entre los momentos de inercia. Visualizarlo en STK: stkSetAttitudeCBI('*/Satellite/Pruebas_Actitud', 'Earth', times, q’);

Para ver bien en STK la curva que describen los ángulos y que se ve en la figura 3 de MATLAB, se puede hacer de la siguiente forma: en el satélite, 3D-Graphics, Vector, eliminar todos los vectores y sistemas de referencia presentes, y añadir el eje z de los ejes cuerpo. Marcar show en persistence, connect:trace, y como ejes los J2000. 6. Probar a modificar la velocidad angular del satélite, haciéndola menor (linea 20 del script cuerpoForz) o a cambiar las inercias (linea 12) para ver como se modifican los resultados. Si se baja lo suficiente la velocidad angular inicial (por ejemplo hasta 0.1), se ve que se pierde rigidez giroscópica. Finalmente ejecutar el script cuerpoGG. Se simula el efecto del gradiente gravitatorio en tres casos: sin otros momentos perturbadores (qg1, primera gráfica), sometido a un momento perturbador constante durante diez segundos iniciales (qg2, segunda gráfica), y sometido al mismo momento pero con un volante de inercia rotando con eje de giro perpendicular al plano orbital (qg3, tercera gráfica). Comparar los tres casos y visualizarlo en STK, teniendo en pantalla los ejes cuerpo y los ejes VVLH. Para poder simular en STK es necesario escribir el cuaternión de actitud (qg1,qg2 o qg3) en relación a la órbita para lo cual hay que ejecutar qg1_orb=qLVLH(times,a,raan,inc,omegaperi,theta0,qg1);

y ya se puede visualizar en STK: stkSetAttitudeCBI('*/Satellite/Pruebas_Actitud', 'Earth', times, qg1_orb);

Igualmente para qg2 y qg3, pero el vector de tiempo correcto para estos vectores de actitud es T. OJO: La duración de la simulación es muy corta. 7. ¿Cuál es la posición de equilibrio (qué ejes están apuntando en qué dirección)? ¿A qué se debe la diferencia entre el caso 2 y 3?