Modelado y simulación de un vehículo auscultador para una vía ferroviaria a escala

i Proyecto Fin de Grado Grado en Ingeniería de las Tecnologías Industriales Intensificación Mecánica-Máquinas Modelado y simulación de un vehículo auscultador para una vía ferroviaria a escala. Autor: Miguel Rodríguez Gómez Tutor: José Luis Escalona Franco Dpto. de Ingeniería Mecánica y Fabricación Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 20 2 1 iii Proyecto Fin de Grado Grado en Ingeniería de las Tecnologías Industriales Modelado y simulación de un vehículo auscultador para una vía ferroviaria a escala. Autor: Miguel Rodríguez Gómez Tutor: José Luis Escalona Franco Catedrático de Universidad Dpto. de Ingeniería Mecánica y Fabricación Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2021 v Proyecto Fin de Grado: Modelado y simulación de un vehículo auscultador para una vía ferroviaria a escala. Autor: Miguel Rodríguez Gómez Tutor: José Luis Escalona Franco El tribunal nombrado para juzgar el Proyecto arriba indicado, compuesto por los siguientes miembros: Presidente: Vocales: Secretario: Acuerdan otorgarle la calificación de: Sevilla, 2021 El Secretario del Tribunal vii A mi familia A mis amigos A mis profesores ix Agradecimientos uiero agradecer a todas las personas que me han acompañado hasta ahora y me han apoyado para llegar aquí. A mi familia, en especial, a mis padres, que me han educado y apoyado, abriéndome oportunidades que me han hecho llegar a donde estoy ahora. A mis compañeros y amigos, con los que he disfrutado los buenos momentos, y los que me han ayudado en los malos. A mis profesores, quienes me han enseñado conocimientos teóricos y también a cómo ser crítico y tomar decisiones en la vida. A mi tutor José Luis, que ha estado disponible en cualquier momento y me ha brindado la oportunidad de aprender realizando este proyecto. A Pedro Urda, que se ha puesto a completa disposición para lo que necesite. Miguel Rodríguez Gómez Sevilla, 2021 Q 5.3. Desarrollo del programa 29 5.4. Resultados 30 5.5. Conclusiones 32 Referencias 35 Anexo A Códigos de Matlab 37 A.1 Archivo del cálculo simbólico “auscultador_simbolico.m” 37 A.2 Archivo de Parámetros “Parametros.m” 42 A.3 Archivo a ejecutar para simular “simulacion.m” 43 A.4 Archivo con las ecuaciones de movimiento “EcsMovimiento.m” 49 A.5 Archivo para el cálculo de la fuerza generalizada de rozamiento “F_roz.m” 50 A.6 Archivo para el post-procesado de los resultados “PostProceso.m” 51 xvii ÍNDICE DE TABLAS Tabla 2.1 Vectores de Posición Relativa de los Centros de Gravedad 6 Tabla 2.2 Vectores posición relativa entre centros de gravedad y centros de ruedas 6 Tabla 2.3 Radios de las Ruedas 7 Tabla 2.4 Vectores de Posición Relativa de los Extremos del Muelle 7 Tabla 2.5 Valores de Masa e Inercia de los Sólidos 9 xix ÍNDICE DE FIGURAS Figura 1.1 Ancho de vía 1 Figura 1.2 Nivelación 1 Figura 1.3 Alineamiento 2 Figura 1.4 Peralte 2 Figura 2.1 Identificación y Numeración de Sólidos 4 Figura 2.2 Coordenadas del Sistema 5 Figura 2.3 Vectores para definir la Posición de los Centros de Gravedad y Ruedas. 7 Figura 2.4 Vista frontal Muelle 8 Figura 2.5 Vista Inferior Tren Delantero 8 Figura 2.6 Vista lateral 8 Figura 3.1 Sistema de referencia de la línea media 11 Figura 3.2 Sistemas de Referencia de la Vía 12 Figura 3.3 Vectores para localizar los centros de los carriles 12 Figura 3.4 Puntos de contacto del vehículo con los carriles 14 Figura 4.1 Ley de Coulomb 21 Figura 4.2 Ley de Creep-Coulomb 21 Figura 5.1 Geometría de la línea media de la vía 26 Figura 5.2 Ángulo de yaw del sistema de referencia de la vía 26 Figura 5.3 Ángulo de pitch del sistema de referencia de la vía 26 Figura 5.4 Ángulo de roll del sistema de referencia de la vía 26 Figura 5.5 Nivelación de la vía 27 Figura 5.6 Peralte de la vía 27 Figura 5.7 Ancho de la vía 27 Figura 5.8 Alineamiento de la vía 27 Figura 5.9 Trayectoria del centro de gravedad del chasis del vehículo. 30 Figura 5.10 Ángulo de yaw del sólido 2 30 Figura 5.11 Ángulo de pitch del sólido 2 30 Figura 5.12 Ángulo de roll del sólido 2. 30 Figura 5.13 Ángulo de yaw del sólido 3 31 Figura 5.14 Ángulo de pitch del sólido 3 31 Figura 5.15 Ángulo de roll del sólido 3. 31 Figura 5.16 Distancia transversal entre los sólidos 4 y 5 31 Figura 5.18 Desviaciones del ángulo de pitch del tren delantero 32 Figura 5.17 Desviaciones del ángulo de yaw del tren delantero 32 xxi Notación  Ángulo de giro respecto a un eje X (Ángulo de roll)  Ángulo de giro respecto a un eje Y (Ángulo de pitch)  Ángulo de giro respecto a un eje Z (Ángulo de yaw) r Vector de posición en coordenadas globales r Vector de posición en coordenadas locales i j G G Vector que une los centros de gravedad de los sólidos i y j en coordenadas locales i j G O Vector que une el centro de gravedad del sólido i y el centro de la rueda j en coordenadas locales i Índice para indicar el número del sólido k Rigidez del muelle 0 l Longitud natural del muelle m Masa de un sólido G I  Tensor de inercia de un sólido expresado en el sistema de referencia colocado en su centro de gravedad y solidario al sólido q Vector de coordenadas generalizadas q  Primera derivada respecto al tiempo q  Segunda derivada respecto al tiempo s Vector de coordenadas no generalizadas t Tiempo C Ecuaciones de restricción q C Jacobiano de las restricciones respecto a las coordenadas generalizadas s C Jacobiano de las restricciones respecto a las coordenadas no generalizadas t C Derivada parcial de las restricciones respecto al tiempo cos Función coseno sen Función seno , i j A Matriz de rotación para girar un vector expresado según el sistema de referencia j al sistema de referencia i i A Matriz de rotación para girar un vector expresado según el sistema de referencia i al sistema de referencia fijo global H Jacobiano de las traslaciones de un sólido en su centro de gravedad en coordenadas globales h Derivada temporal del jacobiano de las traslaciones de un sólido en su centro de gravedad en coordenadas globales G Jacobiano de las rotaciones de un sólido expresado en las coordenadas locales del sólido g Derivada temporal del jacobiano de las rotaciones de un sólido en coordenadas locales del sólido ω Velocidad angular de un sólido en coordenadas locales ω  Matriz antisimétrica asociada al vector de velocidad angular de un sólido en coordenadas locales L Distancia entre la línea media de la vía y el centro de la sección de un carril en su posición de diseño n m  0 Matriz de ceros de tamaño n m  n n  I Matriz identidad de tamaño n n  P a Aceleración de un punto P en coordenadas globales α Aceleración angular de un sólido en coordenadas locales F Resultante de las fuerzas aplicadas a un sólido en coordenadas globales M Resultante de los momentos aplicados a un sólido en coordenadas locales ˆ M Matriz de masa no generalizada M Matriz de masa generalizada ˆ Q Vector de fuerza no generalizada Q Vector de fuerza generalizada U Energía potencial del muelle  Multiplicador de Lagrange R Radio de una rueda x Coordenada X y Coordenada Y z Coordenada Z s Longitud de arco F   q Derivada parcial de F respecto de q dF dt Derivada absoluta de F respecto al tiempo f  Derivada primera espacial de f f  Derivada segunda espacial de f T A Traspuesta de una matriz  a b Producto vectorial de los vectores a y b  Para todo   , a b  Pertenece al intervalo definido por los extremos a y b , incluyéndolos < Menor que > Mayor que 1 1 INTRODUCCIÓN l ferrocarril es un medio de transporte muy importante en la actualidad, y se ha ido desarrollando durante varios siglos. Con los avances tecnológicos, las velocidades de los trenes han ido aumentando, pudiendo llegar fácilmente a los 300 km/h. Con ello, también han aumentado las medidas de seguridad, entre las que se encuentra la calidad de la vía. Así, para poder garantizar que la vía sea segura de circular, es necesario poder determinar las irregularidades que aparecen en ella por el uso, siendo necesario la auscultación geométrica de la misma. En este sentido, se está desarrollando un proyecto de investigación [3] titulado Desarrollo de observadores en sistemas embarcados para la estimación de la respuesta dinámica de vehículos ferroviarios, en el que se pretende analizar mediante sensores instalados en el vehículo su dinámica, y a partir de ella poder estimar las irregularidades de las vías. Existen 4 irregularidades, dos relativas y dos absolutas: - Ancho de Vía: Determina la distancia entre los carriles, siendo entonces una medida relativa. Conociendo las coordenadas de posición de los carriles izquierdo y derecho, se puede determinar según la expresión: izq der Ancho y y   (1.1) Figura 1.1 Ancho de vía - Nivelación: Diferencia de cota entre los carriles izquierdo y derecho. Es una medida relativa. Se puede calcular como: izq der Nivelación z z   (1.2) Figura 1.2 Nivelación - Alineamiento: Posición horizontal de la línea media de la vía. Es una medida absoluta. Se calcula de la forma: 2 izq der y y Alineamiento   (1.3) E Introducción 2 Figura 1.3 Alineamiento - Peralte: Cota vertical de la línea media de la vía. Es una medida absoluta. Se calcula como: 2 izq der z z Peralte   (1.4) Figura 1.4 Peralte Para llevar a cabo ensayos, se dispone de una vía a escala de geometría conocida, en la que se pueden inducir de forma manual irregularidades, para ser medidas. También, se ha desarrollado y fabricado un vehículo autónomo, cuyo diseño está contenido en el Trabajo Fin de Grado [5], con sensores incorporados, de manera que es capaz de medir directamente las irregularidades relativas, mediante un LVDT para el ancho de vía y un inclinómetro para la Nivelación. En cambio, para las irregularidades absolutas, emplea una estación total, con toda la logística que esta conlleva para preparar. Se quiere desarrollar una forma alternativa de medir las irregularidades absolutas, mediante el uso de una IMU instalada en el vehículo, de manera que sea totalmente autónomo. Para ello, se necesita realizar un modelo del vehículo para poder relacionar las aceleraciones y velocidades angulares medidas por el sensor con las irregularidades que las producen, y es de lo que se encarga este Trabajo Fin de Grado. Para ello, se va a emplear la teoría de dinámica de sistemas multicuerpos, en la que se definen unas coordenadas, las restricciones a las que estén sometidas, y a partir de las cuales se pueden escribir las ecuaciones de movimiento del vehículo. Finalmente, se realizará una simulación usando el programa Matlab®, de la que se podrá extraer la evolución de las coordenadas con el tiempo, y compararlas con la geometría de la vía. En un futuro, se pretende usar este modelo para que, dada la evolución de las coordenadas, se puedan cuantificar las irregularidades absolutas de la vía. 3 2 DEFINICIÓN DEL MODELO CINEMÁTICO DEL VEHÍCULO AUSCULTADOR n este capítulo se va a explicar el funcionamiento del vehículo, y, a partir de él, se definirán los sólidos móviles que lo conforman, y un conjunto mínimo de coordenadas para poder definir la posición de cualquier punto respecto a un sistema de referencia global fijo. Además, se definirán los parámetros geométricos necesarios para ello, así como los parámetros inerciales necesarios para estudiar la dinámica del vehículo. 2.1. Descripción del vehículo El vehículo consta de una estructura de aluminio, que se denominará chasis del vehículo, sobre un conjunto de ruedas dispuestas en un tren trasero y en un tren delantero. El tren trasero es el tren motriz del vehículo. Tiene 2 ruedas de goma, cuyos ejes están conectados a un motor eléctrico a través de una pequeña transmisión por correas y un diferencial, que ajusta los giros de las ruedas convenientemente para poder tomar curvas. Tanto las ruedas como la transmisión y el motor están sujetados en una posición fija al chasis. El tren delantero está unido al chasis a través de un par esférico, pudiendo así girar independientemente al chasis, para adaptarse a las elevaciones y peraltes de la vía. El tren está compuesto por una guía, a la que la suspensión derecha está fijada, y a lo largo de la cual la suspensión izquierda puede deslizar de forma independiente, adaptándose así al ancho de la vía. La suspensión izquierda descansa sobre 2 ruedas verticales de plástico. La más adelantada es la encargada de seguir la trayectoria del carril, y la trasera evita que el tren delantero vuelque al avanzar el vehículo. La suspensión derecha descansa sobre una única rueda de plástico, que sigue la trayectoria del carril. Además, cada suspensión delantera tiene unas ruedas horizontales para que las suspensiones se adapten al ancho de vía y deslicen sobre la guía. Finalmente, las suspensiones delanteras están unidas entre sí a través de un muelle, cuya función es tirar hacia dentro las suspensiones para que las ruedas laterales estén siempre en contacto con el carril a lo largo del movimiento del mismo. En definitiva, el vehículo tiene 7 ruedas, dos motrices montadas en el tren trasero, y 5 montadas en el tren delantero, diseñado de tal forma que, al avanzar el vehículo, las ruedas se puedan adaptar a la geometría variable de los carriles y a sus irregularidades. 2.2. Definición de sólidos Atendiendo al funcionamiento descrito en el apartado anterior, se definen los siguientes sólidos móviles, cada uno con su sistema de referencia local situado en el centro de gravedad de cada sólido: - Sólido ②: Chasis del vehículo. Puede tener cualquier posición y orientación respecto al sistema de referencia global fijo. Su sistema de referencia servirá de punto de partida para encontrar al resto de sólidos. A él pertenecen las dos ruedas motrices. - Sólido ③: Guía en el tren delantero. Conectada al chasis a través de un par esférico situado en el centro de gravedad de la guía, liberando exclusivamente la orientación relativa entre ambos. - Sólido ④: Suspensión izquierda. El único movimiento permitido de este sólido es el deslizamiento o traslación a lo largo de la dirección de la guía. Incluye todas las ruedas montadas en esta suspensión. - Sólido ⑤: Suspensión derecha. Está fijado a la guía. Incluye todas las ruedas montadas en esta suspensión. E 11 3 RESTRICCIONES l vehículo debe rodar por los carriles de la vía, lo que se traduce en restricciones al movimiento. Para definirlas, primero es necesario conocer la geometría de la las líneas que pasan a través de los centros de las secciones de los carriles izquierdo y derecho, y a continuación, se impondrá que los puntos de las ruedas del vehículo pertenezcan a estas líneas, teniendo restricciones de punto a línea. Además, de cara a la simulación del movimiento, es necesario imponer el movimiento de una coordenada a través de una restricción de movilidad. 3.1. Geometría de la vía La geometría la vía está compuesta por la geometría de diseño o ideal, identificada a través de la línea media de la vía y el ancho de vía de diseño, y las irregularidades, superpuestas a la primera. 3.1.1 Sistemas de referencia Para definir la geometría ideal, se definen los sistemas de referencia: - Sistema de referencia de la línea media de la vía , , t t t x y z : Sistema que recorre la línea media de la vía de tal forma que la dirección del eje X es tangente a la vía y apunta en la dirección de avance del vehículo sobre la vía. Su orientación respecto al sistema de referencia global viene determinada por los ángulos de Euler , , t t t    , que representan los giros respecto a los ejes X, Y y Z respectivamente, en orden Z-X-Y: Figura 3.1 Sistema de referencia de la línea media - Sistema de referencia del carril en posición ideal izquierdo , , iizq iizq iizq x y z : Sistema de referencia situado en el centro de la sección del carril izquierdo en su posición ideal. Es paralelo al sistema de referencia de la línea media de la vía. - Sistema de referencia del carril en posición ideal derecho , , ider ider ider x y z : Sistema de referencia situado en el centro de la sección del carril derecho en su posición ideal. Es paralelo al sistema de referencia de la línea media de la vía. La siguiente Figura representa todos los sistemas vistos en el plano que contiene las secciones de los carriles: E Restricciones 12 Figura 3.2 Sistemas de Referencia de la Vía 3.1.2 Vectores de posición Para localizar los centros de las secciones de los carriles, se definen los vectores de la Figura, cada uno expresado respecto al sistema de referencia correspondiente. Cada sección se identifica mediante el parámetro arco s, que indica la longitud recorrida desde el punto de referencia inicial hasta la sección determinada. Figura 3.3 Vectores para localizar los centros de los carriles 13 Modelado y simulación de un vehículo auscultador para una vía ferroviaria a escala. - Vector t r  : Vector de posición del sistema de referencia de la línea media de la vía, expresado respecto al sistema global:       T t x y z R s R s R s      r (3.1) - Vector iizq r  : Vector de posición del carril izquierdo en su posición ideal o de diseño. Expresado respecto al sistema de referencia de la línea media de la vía:   0 0 T iizq Lr (3.2) Donde 67.5 L  mm. - Vector ider r  : Vector de posición del carril derecho en su posición ideal o de diseño. Expresado respecto al sistema de referencia de la línea media de la vía:   0 0 T ider L r (3.3) - Vector rizq r  : Vector de posición del carril izquierdo en su posición real. Expresado respecto al sistema de referencia del carril izquierdo en posición ideal:     0 T rizq izq izq y s z s      r (3.4) - Vector rder r  : Vector de posición del carril derecho en su posición real. Expresado respecto al sistema de referencia del carril derecho en posición ideal:     0 T rder der der y s z s      r (3.5) Para poder sumar los vectores para llegar a la posición de los puntos P y Q, es necesario expresarlos todos en coordenadas globales, usando las matrices de rotación de cada sistema de referencia. 3.1.3 Matrices de rotación - Sistema de referencia de la línea media de la vía: Su orientación está determinada por los ángulos de Euler , , t t t    , definidos en la Figura 3.1. Así, cada ángulo tiene una rotación simple asociada, cuyas matrices de rotación son:                         1 0 0 cos 0 sen cos sen 0 0 cos sen 0 1 0 sen cos 0 0 sen cos sen 0 cos 0 0 1 t t t t t t t t t t t t t t t s s s s s s s s s s s s                                                    A A A (3.6) Por tanto, la matriz de rotación asociada a este sistema es: 1, t t t t t      A A A A A (3.7) - Sistemas de referencia de los carriles izquierdo y derecho en posición ideal: Estos sistemas son paralelos al sistema de la línea media de la vía, por lo que su matriz de rotación es idéntica: ideal t  A A (3.8) Con todo, las posiciones de los centros de las secciones de los carriles en coordenadas globales son:     izq t t iizq ideal rizq der t t ider ideal rder s s       r r A r A r r r A r A r (3.9) Restricciones 14 3.2. Restricciones por contacto El vehículo tiene 5 ruedas verticales que deben seguir a las curvas tridimensionales definidas anteriormente, que representan las curvas que siguen los centros de los carriles izquierdo y derecho respectivamente. Es decir, hay 5 puntos, denominados P, Q, S, T, U, uno de cada rueda, que están obligados a pertenecer a dichas curvas según se corresponda con una rueda en la parte izquierda o en la parte derecha del vehículo, mirando en la dirección de avance del vehículo (sentido positivo del eje 2 x ): Figura 3.4 Puntos de contacto del vehículo con los carriles Ahora bien, el tren delantero tiene 2 ruedas horizontales que impiden el deslizamiento transversal de las ruedas verticales asociadas sobre la curva, mientras que el resto (las ruedas motrices y la rueda soporte), no cuentan con dicha restricción, pudiendo deslizar de forma transversal sobre las curvas. Esto da lugar a dos pares cinemáticos diferentes. 3.2.1 Pares cinemáticos de los puntos P, Q y U Los puntos P y Q pertenecen a los sólidos ④ y ⑤ respectivamente, y representan los puntos de contacto de las ruedas verticales de centros 3 O y 4 O . Como estas ruedas tienen asociadas a las ruedas horizontales, no pueden deslizar transversalmente sobre las curvas de los carriles, por lo que dichos puntos deben pertenecer a la curva directamente:     P P izq P Q Q der Q s s       r r 0 C C r r 0 (3.10) Por otra parte, las ruedas motrices traseras en principio pueden deslizar transversalmente sobre las curvas de los carriles. Para que este movimiento sea estable, es necesario que se controle de alguna forma, introduciendo, por ejemplo, la fuerza de rozamiento al deslizamiento transversal de la rueda, o también impidiendo el deslizamiento de una de las ruedas. Para elegir entre las dos opciones presentadas, hay que tener en cuenta que la única función del tren motriz es la de hacer que el vehículo avance, y como las orientaciones del chasis y del tren delantero auscultador están liberados a través del par esférico, los deslizamientos transversales del tren motriz no afectan a la dinámica del tren delantero. Entonces, se opta por la segunda opción, que es más sencilla pero menos exacta, aunque esta pérdida de exactitud no interfiere con el tren delantero. Así, se elige la rueda motriz derecha, cuyo punto de contacto es el punto U, perteneciente al sólido ②, por lo que la restricción asociada es: 15 Modelado y simulación de un vehículo auscultador para una vía ferroviaria a escala.   U U der U s    C r r 0 (3.11) Los vectores de posición de los puntos P, Q y U se pueden obtener sabiendo los vectores de posición de los centros de gravedad de los sólidos y las matrices de rotación asociadas a cada uno, que se calculan en el capítulo siguiente, de la forma:             4 5 2 4 3 5 5 4 4 2 2 2 2 0 0 0 0 0 0 T P G T Q G T U G R R R             4 3 r r A G O r r A G O r r A G O (3.12) Por tanto, se tienen en total 9 restricciones, en las cuales aparecen 3 coordenadas no generalizadas adicionales a las coordenadas generalizadas previamente definidas: P s , Q s y U s que corresponden a las longitudes de arco de los puntos de contacto pertenecientes a las curvas de los carriles. 3.2.2 Pares cinemáticos de los puntos S y T El punto S se corresponde con el punto de contacto de la rueda soporte, y pertenece al sólido ④. El punto T se corresponde con el punto de contacto de la rueda motriz izquierda, y pertenece al sólido ②. Estas ruedas, de forma cilíndrica, son capaces de deslizar transversalmente sobre las curvas de los raíles: Las restricciones asociadas a estos puntos tienen la misma forma que las anteriores, ya que un punto de la rueda debe coincidir con un punto de la curva del raíl:     S S izq S T T izq T s s       C r r 0 C r r 0 (3.13) Pero ahora, la posición de los puntos de contacto de las ruedas es diferente:             4 2 4 5 2 1 0 0 0 0 0 0 0 0 T T S G S T T T G T R y R y           4 5 2 1 r r A G O r r A G O (3.14) Entonces, se tienen 6 restricciones adicionales, donde aparecen 4 nuevas coordenadas no generalizadas: S s y T s que representan las longitudes de arco de los puntos de contacto pertenecientes a las curvas de los carriles, e , S T y y , que representan la coordenada Y de los puntos de contacto medidas desde el centro de cada rueda, de forma que cuantifican el deslizamiento transversal sobre las curvas. 3.3. Restricción de movilidad En total, se han definido hasta ahora 15 restricciones, y a las 10 coordenadas generalizadas definidas en (2.1) se le añaden 7 coordenadas no generalizadas, que se agrupan en el vector: S T P Q S T U y y s s s s s      s (3.15) Por tanto, el vehículo tiene un total de 10 7 15 2    grados de libertad, que se corresponden al avance del vehículo sobre la vía, mediante la coordenada 2 G x y al ángulo de yaw del tren delantero   3  , gracias al deslizamiento transversal de la rueda soporte. Para poder realizar la simulación, es necesario imponer uno de estos grados de libertad, para lo que se elige el avance del vehículo, imponiendo la restricción de movilidad: Restricciones 16   exp 0 mov T T C s s t    (3.16) Donde la función   exp T s t es una función obtenida en uno de los ensayos del vehículo. En definitiva, el modelo a simular está descrito con 17 coordenadas, y tiene 1 grado de libertad, por estar sometido a las 16 restricciones:   , , mov P Q S T U C t                       C C C q s 0 C C C (3.17) 17 4 CÁLCULO SIMBÓLICO DE LAS ECUACIONES DE MOVIMIENTO DEL VEHÍCULO na vez definido el modelo cinemático del vehículo, hay que calcular de manera simbólica sus ecuaciones de movimiento. Para ello, hay que definir las posiciones de todos los sistemas de referencia y las matrices de rotación asociadas a cada uno, para determinar los jacobianos que transforman las ecuaciones de Newton-Euler en las ecuaciones de movimiento en coordenadas generalizadas, y permiten también calcular las fuerzas generalizadas. 4.1. Cinemática 4.1.1 Matrices de rotación Para empezar, es necesario calcular las matrices de rotación de cada sólido, que transformen un vector expresado según el sistema de referencia local al sistema de referencia global, denotado como sistema 1. Para ello, se usan rotaciones sucesivas empleando matrices de giro intermedias cuando sean necesarias. - Sólido ②: La orientación de este sólido viene definida por las coordenadas 2 2 2 , ,    , es decir, por sus ángulos de Euler. Así, cada ángulo representa una rotación simple alrededor de los ejes X, Y, y Z respectivamente, en el orden Z- Y-X. Las matrices de las rotaciones simples asociadas a cada ángulo son: 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 1 0 0 cos 0 sen cos sen 0 0 cos sen 0 1 0 sen cos 0 0 sen cos sen 0 cos 0 0 1                                                    A A A (4.1) Por lo tanto, la matriz de rotación del sólido es: 2 2 2 2 1,2      A A A A A (4.2) - Sólido ③: La orientación de este sólido viene definida por las coordenadas 3 3 3 , ,    , que representan los ángulos de Euler medidos respecto al sistema de referencia global. Análogamente al sólido ②, las rotaciones simples son: 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 1 0 0 cos 0 sen cos sen 0 0 cos sen 0 1 0 sen cos 0 0 sen cos sen 0 cos 0 0 1                                                    A A A (4.3) Por tanto, la matriz de rotación del sólido es: 3 3 3 3 1,3      A A A A A (4.4) - Sólidos ④ y ⑤: Estos sólidos no tienen coordenadas asociadas a su orientación, pues son paralelos al sistema 3, ya que el único movimiento que tiene permitido el sólido ④ es la traslación a lo largo del eje Y3, y el sólido ⑤ está fijado al sólido ③, por lo que ambos son paralelos al sistema de referencia 3 Entonces: U Cálculo Simbólico de las Ecuaciones de Movimiento del Vehículo 18 4 3 5 3   A A A A (4.5) 4.1.2 Posición de los centros de gravedad de los sólidos A continuación, se calcula la posición de los centros de gravedad de todos los sólidos a partir del sistema de referencia global. En primer lugar, la posición del centro de gravedad del sólido ② viene directamente determinada por las coordenadas 2 2 2 , , G G G x y z , de forma que: 2 2 2 T G G G x y z      2 G r (4.6) A partir de esta, se pueden determinar el resto haciendo uso de los vectores definidos en el apartado 2.4.1, girándolos utilizando la matriz de rotación adecuada, de tal forma que se obtienen: 2 5 3 4 5 2 3 5          3 G G 2 3 G G 3 5 G G 5 4 r r A G G r r A G G r r A G G (4.7) 4.1.3 Jacobianos de traslación y rotación. Matriz de transformación Una vez conocidas la posición y matriz de rotación de todos los sólidos, se puede proceder al cálculo de los jacobianos. El jacobiano de las traslaciones del sólido i , donde i  2, 3, 4, 5 viene dado por la expresión: i i    r H q (4.8) La derivada temporal del jacobiano de las traslaciones del sólido i es:   i i i      H q H h q   (4.9) El jacobiano de las rotaciones del sólido i en coordenadas locales viene dado por la expresión: i i    ω G q  (4.10) La derivada temporal del jacobiano de las rotaciones del sólido i es:   i i i i         G q ω G g q q   (4.11) Para los dos últimos es necesario calcular previamente la velocidad angular de cada sólido, expresada en coordenadas locales. Utilizando las matrices de rotación, se tiene: T i i i  ω A A   (4.12) Donde i ω  es la matriz antisimétrica asociada al vector i ω . Entonces, se puede construir este vector de la forma: 19 Modelado y simulación de un vehículo auscultador para una vía ferroviaria a escala.       3,2 1,3 2,1 T i i i i      ω ω ω ω    (4.13) Finalmente, para calcular la derivada temporal de la matriz de rotación del sólido i , se utiliza la expresión: 10 1 i i j jj d q dq  A A   (4.14) Una vez calculados los jacobianos y sus derivadas, se define la matriz de transformación y su derivada temporal como: 2 2 2 2 3 3 3 3 5 5 5 5                                                H h G g H h L G L l g H h G g    (4.15) 4.2. Dinámica La dinámica del sólido i viene determinada por sus ecuaciones de Newton-Euler: 3 3 3 3 3 1 3 3 ˆ ˆ ˆ i i i i i G i i i G i G i i m                                   apl V a I 0 0 FM a Q Q 0 I ω I ω M α   (4.16) Donde se utiliza la aceleración angular en coordenadas locales para que el tensor de inercia, en dichas coordenadas, sea constante y no dependa del tiempo. En contraposición, aparecen las fuerzas de inercia cuadráticas en velocidad   ˆ i V Q , que se trata como una fuerza aplicada adicional. Aplicando la transformación definida en (4.15) a las matrices resultantes de agrupar las ecuaciones todos los sólidos, se llega a las ecuaciones de movimiento en coordenadas generalizadas. 4.2.1 Matriz de masa Partiendo de la matriz de masa ˆ i M de cada sólido, se puede montar la matriz del sistema completo: 2 6 6 6 6 6 6 6 6 3 6 6 6 6 6 6 6 6 4 6 6 6 6 6 6 6 6 5 ˆ ˆ ˆˆ ˆ                            M 0 0 0 0 M 0 0 M0 0 M 0 0 0 0 M (4.17) Finalmente, aplicando la transformación, se obtiene la matriz de masa en coordenadas generalizadas: ˆ T  M L ML (4.18) 4.2.2 Fuerzas aplicadas generalizadas En el modelo aparecen 3 fuerzas externas: la gravedad, la fuerza del muelle y las fuerzas de inercia cuadráticas en velocidad. Además, al tener restricciones, aparecen fuerzas de reacción. La Simulación 26 - Cuarta columna: Ancho de vía. - Quinta columna: Peralte.  Archivo “velProfile.txt” Este archivo contiene información sobre el movimiento del vehículo, en unidades del Sistema Internacional. Contiene 4 columnas: - Primera columna: Coordenada   s t . Longitud de arco de la línea media de la vía, como en el caso anterior. - Segunda columna:   v t . Velocidad de avance del vehículo. - Tercera columna:   a t . Aceleración de avance del vehículo. - Cuarta columna: t . Tiempo transcurrido desde el inicio del viaje. Con estos datos, se pueden representar: Figura 5.1 Geometría de la línea media de la vía Figura 5.4 Ángulo de roll del sistema de referencia de la vía Figura 5.3 Ángulo de pitch del sistema de referencia de la vía Figura 5.2 Ángulo de yaw del sistema de referencia de la vía t (rad) t (rad) t (rad) 27 Modelado y simulación de un vehículo auscultador para una vía ferroviaria a escala. 5.2. Cálculos previos Los datos proporcionados no son suficientes para poder realizar la simulación, pero a partir de ellos, se pueden calcular el resto. 5.2.1 Irregularidades Los vectores rizq r y rder r contienen la información sobre las irregularidades de las vías, pero no de la forma en la que se nos proporciona como dato. Aun así, según las expresiones del cálculo de las irregularidades expuestas en el capítulo 1, es fácil pasar de una a otra resolviendo el sistema:                 1 1 0 0 Alineamiento 2 2 Ancho de vía 1 1 0 0 Peralte 1 1 0 0 2 2 Nivelación 0 0 1 0 izq der izq der y s s y s s z s s z s s                                                 (5.1) -10 0 10 20 30 40 50 60 70 80 90 Longitud de arco (m) -6 -5 -4 -3 -2 -1 0 1 2 310-3 Peralte -10 0 10 20 30 40 50 60 70 80 90 Longitud de arco (m) -10 -8 -6 -4 -2 0 2 4 610-3 Nivelación -10 0 10 20 30 40 50 60 70 80 90 Longitud de arco (m) -2 -1.5 -1 -0.5 0 0.5 1 1.5 10-3 Ancho de vía -10 0 10 20 30 40 50 60 70 80 90 Longitud de arco (m) -0.02 -0.015 -0.01 -0.005 0 0.005 0.01 0.015 0.02 Alineamiento Figura 5.5 Nivelación de la vía Figura 5.6 Peralte de la vía Figura 5.7 Ancho de la vía Figura 5.8 Alineamiento de la vía La Simulación 28 5.2.2 Derivadas espaciales Como en las ecuaciones de movimiento aparecen la derivada segunda de las restricciones, es necesario conocer las derivadas espaciales de los vectores izq r y der r . Como dichas funciones se calculan de forma numérica, el cálculo de las derivadas se hará usando diferencias finitas centradas, de tal forma que:             1 1 1 1 1 1 1 1 i i i i i i i i i i s s ss s s s ss s                r r r r r r (5.2) 5.2.3 Interpolación A la hora de simular, el integrador numérico tiene que ser capaz de evaluar todas las funciones que necesite en cualquier punto. Como todas las funciones asociadas a la geometría de la vía, así como las funciones de la restricción de movilidad están definidas por puntos, es necesario interpolar. Para ello, se van a emplear splines, utilizando el comando “spline” de Matlab. 5.2.4 Condiciones iniciales Para cerrar el problema definido por las ecuaciones diferenciales-algebraicas de movimiento del vehículo, es necesario definir las condiciones iniciales de posición y velocidad. Para ello, se supone que el vehículo parte del reposo desde una posición tal que el centro de gravedad del sólido ② esté sobre el punto cero de la vía, es decir, donde el parámetro arco es nulo. Como parte del reposo, las condiciones iniciales en velocidad son: 0 0   q 0 s 0   (5.3) En cuanto a las condiciones iniciales de posición, las coordenadas no generalizadas son fáciles de determinar. Como el tramo en el que empieza a circular el vehículo es una recta, cuya dirección está alineada con el eje X del sistema de referencia global, los parámetros arco de cada punto corresponden a las coordenadas X de dichos puntos respecto del sistema de referencia del sólido ②, que es esta posición es paralelo al global. Además, se supone que el vehículo se coloca de tal forma que el deslizamiento transversal de las ruedas es nulo. Así, en metros se tiene que: 0 0 0 0 0 0 0 2 3 3 5 5 4 4 3 2 3 3 5 5 4 2 3 3 5 5 4 4 5 2 1 2 2 0 0 0.169 0.169 0.120 0.100 0.100 S T P x x x x Q x x x S x x x x T x U x y y s G G G G G G G O s G G G G G O s G G G G G G G O s G O s G O                       (5.4) Ahora bien, para las coordenadas generalizadas, no es fácil en primera instancia calcularlas, debido a que el vehículo se debe adaptar a las irregularidades desde un principio. Para facilitar el cálculo, se utiliza una función auxiliar que sea nula hasta aproximadamente lo que ocupe el vehículo, es decir, hasta 0.2 s  m, y pase a valer 1 de forma suave en un intervalo de s pequeño, como pueden ser 0.3 m. Para hacer la transición de forma suave, se emplea un polinomio de tercer grado, en el que se impone que la derivada en s  0.2 y 0.5 sea nula, y que 29 Modelado y simulación de un vehículo auscultador para una vía ferroviaria a escala. para 0.2 s  la función valga 0, y para 0.5 s  la función valga 1. Imponiendo estas condiciones, se obtiene el polinomio:     3 2 2000 700 200 52 0.2,0.5 27 9 9 27 f s s s s s        (5.5) Por tanto, multiplicando esta función por las funciones de las irregularidades         , , , izq ider izq ider y s y s z s z s hace que dichas funciones sean nulas en la posición inicial del vehículo, y que vuelvan a ser iguales a partir de 0.5 s  m. Entonces, como el vehículo se encuentra en una recta alineada con el eje X global, todos los sistemas de los sólidos son paralelos entre sí, y paralelos al sistema global, siendo nulos todos los ángulos de giro, y el resto de coordenadas son sencillas de calcular, siendo: 2 2 2 0 0 0 0 0 0 0 0 2 1 0 1 2 1 2 2 2 3 3 3 0 4 3 5 4 0m 0.0045m 0.06m 0 rad 0 rad 0 rad 0 rad 0 rad 0 rad 2 0.137m G G y G z y y x y L G O z R G O w L G O G O                         (5.6) 5.3. Desarrollo del programa El programa principal de la simulación es el llamado “simulacion.m”, y es el que se debe ejecutar. En él, primero se declaran todos los parámetros necesarios a través de la función “Parametros.m”. A continuación, se cargan los archivos con los datos de partida, y luego se realizan todos los cálculos previos definidos en el apartado anterior. Después, integra las ecuaciones de movimiento, definidas en la función “EcsMovimiento.m”, y finalmente realiza el post-procesado de los resultados la función “PostProceso.m”, y se hace una animación del movimiento del vehículo mediante la función “Animacion.m”. La función “EcsMovimiento.m” resuelve el sistema algebraico que resulta en cada iteración de la integración numérica de las ecuaciones diferenciales. Para ello, hace uso de varias funciones: - “mat_L.m”: Calcula la matriz de transformación L . - “mat_dL.m”: Calcula la derivada de la matriz de transformación l . - “mat_MM”: Calcula la matriz de masa de las ecuaciones de Newton-Euler ˆ M . - “F_inerc_NO_General.m”: Calcula las fuerzas de inercia cuadráticas en velocidad ˆ V Q . - “F_gravedadGeneral.m”: Calcula las fuerzas generalizadas debidas a la gravedad g Q . - “F_muelleGeneral.m”: Calcula las fuerzas generalizadas debidas al muelle k Q . - “F_roz.m”: Calcula las fuerzas generalizadas debidas al rozamiento transversal roz Q . La Simulación 30 - “Cq.m”: Calcula el jacobiano de las restricciones q C . - “dCq.m”: Calcula la derivada temporal del jacobiano de las restricciones q C  . - “Cs.m”: Calcula el jacobiano de las restricciones s C . - “dCs.m”: Calcula la derivada temporal del jacobiano de las restricciones s C  . - “dCt.m”: Calcula la derivada temporal de la parcial de las restricciones respecto al tiempo t C  . Todas estas funciones han sido generadas a partir del archivo de cálculo simbólico “auscultador_simbolico.m”, menos la función “F_roz.m”. Esta función realiza todos los cálculos expuestos en la sección 4.2.2.3, y utiliza la función “jacobS.m” para calcular el jacobiano de las traslaciones en el punto S de contacto, que también ha sido generada a partir del cálculo simbólico. 5.4. Resultados Realizando la simulación durante el trayecto de 50 segundos, se obtienen los resultados: Figura 5.9 Trayectoria del centro de gravedad del chasis del vehículo. Figura 5.12 Ángulo de roll del sólido 2. Figura 5.11 Ángulo de pitch del sólido 2 Figura 5.10 Ángulo de yaw del sólido 2 2 (rad) 2 (rad) 0 5 10 15 20 25 30 35 40 45 50 Tiempo (s) -1.5 -1 -0.5 0 0.5 1 1.5 2 (rad) 2(t) 31 Modelado y simulación de un vehículo auscultador para una vía ferroviaria a escala. Figura 5.16 Distancia transversal entre los sólidos 4 y 5 Se puede observar comparando estas gráficas con las gráficas de la sección 5.1, que el vehículo sigue correctamente a la vía, validando que la simulación es correcta. Por otra parte, la coordenada w está directamente relacionada con el ancho de vía, por su definición. Además, al ser el ángulo t  nulo en todo el trayecto, el ángulo de roll del tren delantero   3  está directamente relacionado con la nivelación. Además, si a los dos ángulos restantes del tren delantero se le restan los ángulos debidos a seguir la geometría de diseño de la vía, se obtienen las desviaciones que deben entonces estar relacionadas con las irregularidades absolutas: 0 5 10 15 20 25 30 35 40 45 50 Tiempo (s) 0.135 0.1355 0.136 0.1365 0.137 0.1375 0.138 0.1385 0.139 w (m) w(t) Figura 5.15 Ángulo de roll del sólido 3. Figura 5.14 Ángulo de pitch del sólido 3 Figura 5.13 Ángulo de yaw del sólido 3 3 (rad) 0 5 10 15 20 25 30 35 40 45 50 Tiempo (s) -0.06 -0.05 -0.04 -0.03 -0.02 -0.01 0 0.01 0.02 3 (rad) 3(t) 3 (rad) La Simulación 32 Se observa que las desviaciones de 3  guardan relación con el peralte de la vía, pero las desviaciones de 3  no tienen relación aparente con el alineamiento. 5.5. Conclusiones Una vez realizada la simulación, se puede afirmar que el modelo empleado representa con bastante precisión al vehículo real. El modelo utilizado es una simplificación del problema de contacto entre superficies estudiado mediante el método de las restricciones por contacto, donde se tienen en cuenta la geometría de las superficies del carril y las ruedas. En el modelo no se han utilizado las superficies de contacto, sino que se han simplificado por líneas, pasando de restricciones de contacto entre dos superficies a restricciones de punto a línea. Adicionalmente, el problema del contacto puede ser analizado usando el método elástico, donde las fuerzas de contacto se calculan previamente y se aplican al vehículo, lo que requiere del cálculo de estas fuerzas, tanto normales como tangenciales, teniendo en cuenta la geometría de las superficies de contacto. Este problema, más complicado, representa más verazmente la realidad, pero no ha sido necesario utilizarlo. Según los resultados obtenidos, se demuestra que las irregularidades en la vía tienen un efecto visible en la dinámica del vehículo, por lo que el objetivo siguiente de averiguar las irregularidades a partir de la dinámica del vehículo es posible. 0 5 10 15 20 25 30 35 40 45 50 Tiempo (s) -0.1 -0.05 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 Desviaciones de 3 0 5 10 15 20 25 30 35 40 45 50 Tiempo (s) -0.06 -0.05 -0.04 -0.03 -0.02 -0.01 0 0.01 0.02 0.03 Desviaciones de 3 Figura 5.18 Desviaciones del ángulo de yaw del tren delantero Figura 5.17 Desviaciones del ángulo de pitch del tren delantero 33 Modelado y simulación de un vehículo auscultador para una vía ferroviaria a escala. 35 REFERENCIAS [1] J. L. Escalona, Dinámica de Sistemas Multicuerpo, Sevilla. [2] J. L. Escalona, A Multibody Formulation for the Dynamic Analysis of Railway Vehicles, Sevilla. [3] R. Chamorro, J. L. Escalona, «Desarrollo de observadores en sistemas embarcados para la estimación de la respuesta dinámica de vehículos ferroviarios.» Memoria Cinetífico- Técnica del Proyecto de Investigación TRA2017, 2017. [4] H. Sugiyama, K. E. Zaazaa, A. A. Shabana, Railroad Vehicle Dynamics. A Computational Approach, New York: Taylor & Francis Group, 2008, p. 12. [5] R. Álvarez. , «Diseño, fabrica ción, montaje y puesta a punto de un vehículo autónomo para la auscultación geométrica de una vía ferroviaria a escala.,» Trabajo Fin de Grado, 2019. Anexo A Códigos de Matlab 42 O, O, O, O, O, O, dsU*ddsU_rder]; dC_s = dC_s_veh - dC_Via; %Rozamiento HS = jacobian(rS,q); %Creación de funciones matlabFunction(MM,'file','mat_MM'); matlabFunction(C_q,'file','Cq'); matlabFunction(dC_q,'file','dCq'); matlabFunction(C_s,'file','Cs'); matlabFunction(dC_s,'file','dCs'); matlabFunction(dC_t,'file','dCt'); matlabFunction(Qg,'file','F_gravGeneral'); matlabFunction(Qk,'file','F_muelleGeneral'); matlabFunction(QQv,'file','F_inerc_NO_General'); matlabFunction(L,'file','mat_L'); matlabFunction(l,'file','mat_dL','Optimize',false); matlabFunction(HS,'file','jacobS'); A.2 Archivo de Parámetros “Parametros.m” function Param = Parametros %Radios ruedas Param.R1 = 40e-3; Param.R2 = 40e-3; Param.R3 = 12e-3; Param.R4 = 12e-3; Param.R5 = 10e-3; %Masas Param.g = 9.81; Param.m2 = 6; Param.m3 = 0.3; Param.m4 = 0.7; Param.m5 = 0.4; %Inercias Param.IG2xx = 2.09e-2; Param.IG2xy = 7.64e-5; Param.IG2xz = 9.09e-4; Param.IG2yy = 5.80e-2; Param.IG2yz = 2.44e-4; Param.IG2zz = 7.22e-2; Param.IG3xx = 1.90e-3; Param.IG3xy = 0; Param.IG3xz = 0; Param.IG3yy = 2.75e-5; Param.IG3yz = 0; Param.IG3zz = 1.89e-3; Param.IG4xx = 3.11e-4; Param.IG4xy = 1.98e-5; Param.IG4xz = -1.56e-5; Param.IG4yy = 1.17e-3; Param.IG4yz = -8.25e-6; Param.IG4zz = 1.46e-3; Param.IG5xx = 1.52e-4; Param.IG5xy = -2.99e-7; Param.IG5xz = -2.12e-8; Param.IG5yy = 1.12e-4; Param.IG5yz = 7.44e-6; Param.IG5zz = 2.47e-4; 43 Modelado y simulación de un vehículo auscultador para una vía ferroviaria a escala. %Vectores de posición Param.G2G3 = [0.165, 0.005, -0.022]'; Param.G3G5 = [0.003 ,-0.069, -0.028]'; Param.G5G4 = [-0.044, 0.137, 0.001]'; Param.G2O1 = [-0.1, 0.072, -0.048]'; Param.G2O2 = [-0.1, -0.063, -0.048]'; Param.G4O5 = [-0.004, -0.001, -0.029]'; Param.G4O3 = [0.045, -0.001, -0.027]'; Param.G5O4 = [0.001, 0.001, -0.026]'; Param.uC = [0.044, -0.046, 0.001]'; Param.uD = [0, 0.046, 0.002]'; %Muelle Param.k = 60; Param.L0 = 19.05e-3; %Via Param.L = 0.0675; Param.mu = 0.4; Param.vlim = 0.05; A.3 Archivo a ejecutar para simular “simulacion.m” global Param Param = Parametros; T_In = load('trackInput.txt'); T_Irr = load('trackIrr.txt'); V_In = load('velProfile.txt'); sV = V_In(:,1); tV = V_In(:,4); aV = V_In(:,3); Param.sV_spline = spline(tV,sV); Param.aV_spline = spline(tV,aV); sIrr = T_Irr(:,1); Alineamiento = T_Irr(:,2); Nivelacion = T_Irr(:,3); Ancho = T_Irr(:,4); Peralte = T_Irr(:,5); sIn = T_In(:,1); phit = T_In(:,2); tett = T_In(:,3); psit = T_In(:,4); Rx = T_In(:,9); Ry = T_In(:,10); Rz = T_In(:,11); Anexo A Códigos de Matlab 44 Param.psit_spline = spline(sIn,psit); Param.tett_spline = spline(sIn,tett); AUX = [0.5, 0.5, 0, 0; 1, -1, 0, 0; 0, 0, 0.5, 0.5; 0, 0, 1, -1]; incog = zeros(length(T_Irr(:,1)),4); for i = 1:length(sIrr) Irr = [Alineamiento(i), Ancho(i), Peralte(i), Nivelacion(i)]'; incog(i,:) = (AUX\Irr)'; end yizq = incog(:,1); yder = incog(:,2); zizq = incog(:,3); zder = incog(:,4); %Suavizar irregularidades al principio j = 1; f = ones(length(sIrr),1); while sIrr(j) <= 0.2 f(j) = 0; j = j+1; end while sIrr(j) > 0.2 && sIrr(j) <= 0.5 f(j) = -(2000/27)*sIrr(j)^3 + (700/9)*sIrr(j)^2 - (200/9)*sIrr(j) + (52/27); j = j+1; end yizq = yizq.*f; zizq = zizq.*f; yder = yder.*f; zder = zder.*f; Param.yizq_spline = spline(sIrr, yizq); Param.yder_spline = spline(sIrr, yder); Param.zizq_spline = spline(sIrr, zizq); Param.zder_spline = spline(sIrr, zder); %Interpolación líneas de los carriles L = Param.L; rizq = zeros(length(sIn),3); rder = zeros(length(sIn),3); 45 Modelado y simulación de un vehículo auscultador para una vía ferroviaria a escala. for i = 1:length(sIn) rt = [Rx(i), Ry(i), Rz(i)]'; riizq = [0, L, 0]'; rider = [0, -L, 0]'; rrizq = [0, ppval(Param.yizq_spline,sIn(i)), ppval(Param.zizq_spline,sIn(i))]'; rrder = [0, ppval(Param.yder_spline,sIn(i)), ppval(Param.zder_spline,sIn(i))]'; Aphit = [1, 0, 0; 0, cos(phit(i)), -sin(phit(i)); 0, sin(phit(i)), cos(phit(i))]; Atett = [cos(tett(i)), 0, sin(tett(i)); 0, 1, 0; -sin(tett(i)), 0, cos(tett(i))]; Apsit = [cos(psit(i)), -sin(psit(i)), 0; sin(psit(i)), cos(psit(i)), 0; 0, 0, 1]; At = Apsit*Atett*Aphit; Aideal = At; rizq(i,:) = (rt + At*riizq + Aideal*rrizq)'; rder(i,:) = (rt + At*rider + Aideal*rrder)'; end %Derivadas por diferencias finitas ds_rizq = zeros(length(sIn),3); ds_rder = zeros(length(sIn),3); for i = 2:length(sIn)-1 ds_rizq(i,:) = (rizq(i+1,:)-rizq(i-1,:))/(sIn(i+1)-sIn(i-1)); ds_rder(i,:) = (rder(i+1,:)-rder(i-1,:))/(sIn(i+1)-sIn(i-1)); end ds_rizq(i+1,:) = ds_rizq(i,:); ds_rizq(1,:) = ds_rizq(2,:); dds_rizq = zeros(length(sIn),3); dds_rder = zeros(length(sIn),3); for i = 2:length(sIn)-1 dds_rizq(i,:) = (ds_rizq(i+1,:)-ds_rizq(i-1,:))/(sIn(i+1)-sIn(i-1)); dds_rder(i,:) = (ds_rder(i+1,:)-ds_rder(i-1,:))/(sIn(i+1)-sIn(i-1)); end dds_rizq(i+1,:) = dds_rizq(i,:); dds_rizq(1,:) = dds_rizq(2,:); Param.rizqx_spline = spline(sIn, rizq(:,1)); Param.rderx_spline = spline(sIn, rder(:,1)); Param.rizqy_spline = spline(sIn, rizq(:,2)); Param.rdery_spline = spline(sIn, rder(:,2)); Param.rizqz_spline = spline(sIn, rizq(:,3)); Param.rderz_spline = spline(sIn, rder(:,3)); Anexo A Códigos de Matlab 46 Param.ds_rizqx_spline = spline(sIn, ds_rizq(:,1)); Param.ds_rderx_spline = spline(sIn, ds_rder(:,1)); Param.ds_rizqy_spline = spline(sIn, ds_rizq(:,2)); Param.ds_rdery_spline = spline(sIn, ds_rder(:,2)); Param.ds_rizqz_spline = spline(sIn, ds_rizq(:,3)); Param.ds_rderz_spline = spline(sIn, ds_rder(:,3)); Param.dds_rizqx_spline = spline(sIn, dds_rizq(:,1)); Param.dds_rderx_spline = spline(sIn, dds_rder(:,1)); Param.dds_rizqy_spline = spline(sIn, dds_rizq(:,2)); Param.dds_rdery_spline = spline(sIn, dds_rder(:,2)); Param.dds_rizqz_spline = spline(sIn, dds_rizq(:,3)); Param.dds_rderz_spline = spline(sIn, dds_rder(:,3)); %Condiciones iniciales xG20 = 2; yG20 = Param.L-Param.G2O1(2); zG20 = Param.R1-Param.G2O1(3); phi20 = 0; tet20 = 0; psi20 = 0; phi30 = 0; tet30 = 0; psi30 = 0; w0 = Param.G5G4(2); yS0 = 0; yT0 = 0; sP0 = Param.G2G3(1)+Param.G3G5(1)+Param.G5G4(1)+Param.G4O3(1); sQ0 = Param.G2G3(1)+Param.G3G5(1)+Param.G5O4(1); sS0 = Param.G2G3(1)+Param.G3G5(1)+Param.G5G4(1)+Param.G4O5(1); sT0 = Param.G2O1(1); sU0 = Param.G2O2(1); q0 = [xG20, yG20, zG20, phi20, tet20, psi20, phi30, tet30, psi30, w0]'; s0 = [yS0, yT0, sP0, sQ0, sS0, sT0, sU0]'; dq0 = zeros(10,1); ds0 = zeros(7,1); y0 = [q0; s0; dq0; ds0]; %Integración numérica t0 = 0; tend = 50; tspan = t0:0.01:tend; options = odeset('RelTol',1e-3,'AbsTol',1e-3); [t, y] = ode15s(@EcsMovimiento, tspan, y0); q = y(:,1:10); s = y(:,11:17); dq = y(:,18:27); ds = y(:,28:34); %Post-Proceso 47 Modelado y simulación de un vehículo auscultador para una vía ferroviaria a escala. [ddq, dds, lambda] = PostProceso(t,q,dq,s,ds); Animacion(t,q); %Representaciones gráficas %Coordenadas figure(2) plot3(q(:,1),q(:,2),q(:,3)) axis('equal') title('Trayectoria') xlabel('x_{G_{2}} (m)') ylabel('y_{G_{2}} (m)') zlabel('z_{G_{2}} (m)') figure(3) plot(t,q(:,4)) title('\phi_{2}(t)') xlabel('Tiempo (s)') ylabel('\phi_{2} (rad)') figure(4) plot(t,q(:,5)) title('\theta_{2}(t)') xlabel('Tiempo (s)') ylabel('\theta_{2} (rad)') figure(5) plot(t,q(:,6)) title('\psi_{2}(t)') xlabel('Tiempo (s)') ylabel('\psi_{2} (rad)') figure(6) plot(t,q(:,7)) title('\phi_{3}(t)') xlabel('Tiempo (s)') ylabel('\phi_{3} (rad)') figure(7) plot(t,q(:,8)) title('\theta_{3}(t)') xlabel('Tiempo (s)') ylabel('\theta_{3} (rad)') figure(8) plot(t,q(:,9)) title('\psi_{3}(t)') xlabel('Tiempo (s)') ylabel('\psi_{3} (rad)') figure(9) plot(t,q(:,10)) title('w(t)') xlabel('Tiempo (s)') ylabel('w (m)') figure(10) plot(t,s(:,1)) Anexo A Códigos de Matlab 48 title('yS') figure(11) plot(t,s(:,2)) title('yT') figure(12) plot(t,s(:,3)) title('sP') figure(13) plot(t,s(:,4)) title('sQ') figure(14) plot(t,s(:,5)) title('sS') figure(15) plot(t,s(:,6)) title('sT') figure(16) plot(t,s(:,7)) title('sU') %Geometría vía figure(20) plot3(Rx,Ry,Rz) axis('equal') title('Línea media Vía') xlabel('R_{x} (m)') ylabel('R_{y} (m)') zlabel('R_{z} (m)') figure(21) plot(sIn,phit) title('\phi_{t}(s)') xlabel('Longitud de arco (m)') ylabel('\phi_{t} (rad)') figure(22) plot(sIn,tett) title('\theta_{t}(s)') xlabel('Longitud de arco (m)') ylabel('\theta_{t} (rad)') figure(23) plot(sIn,psit) title('\psi_{t}(s)') xlabel('Longitud de arco (m)') ylabel('\psi_{t} (rad)') figure(24) plot(sIrr,Alineamiento) title('Alineamiento') xlabel('Longitud de arco (m)') ylabel('Alineamiento (m)') 49 Modelado y simulación de un vehículo auscultador para una vía ferroviaria a escala. figure(25) plot(sIrr,Ancho) title('Ancho de vía') xlabel('Longitud de arco (m)') ylabel('Ancho (m)') figure(26) plot(sIrr,Peralte) title('Peralte') xlabel('Longitud de arco (m)') ylabel('Peralte (m)') figure(27) plot(sIrr,Nivelacion) title('Nivelación') xlabel('Longitud de arco (m)') ylabel('Nivelación (m)') for i =1:length(t) psi_des(i) = q(i,9) - ppval(Param.psit_spline,ppval(Param.sV_spline,t(i))); tet_des(i) = q(i,8) - ppval(Param.tett_spline,ppval(Param.sV_spline,t(i))); end figure(28) plot(t,psi_des) title('Desviaciones de \psi_{3}') xlabel('Tiempo (s)') ylabel('\psi_{3} (rad)') figure(29) plot(t,tet_des) title('Desviaciones de \theta_{3}') xlabel('Tiempo (s)') ylabel('\theta_{3} (rad)') A.4 Archivo con las ecuaciones de movimiento “EcsMovimiento.m” function dy = EcsMovimiento(t,y) t q = y(1:10); s = y(11:17); dq = y(18:27); ds = y(28:34); N = 15; MM = mat_MM; Qg = F_gravGeneral(q); Qk = F_muelleGeneral(q); Anexo A Códigos de Matlab 50 Qroz = F_roz(q,dq,s,N); QQv = F_inerc_NO_General(q,dq); L = mat_L(q); l = mat_dL(q,dq); M = L'*MM*L; Qv = L'*(QQv - MM*l*dq); C_q = Cq(q,s); dC_q = dCq(q,dq,s,ds); dC_t = dCt(t); C_s = Cs(q,s); dC_s = dCs(q,dq,s,ds); A = [M, zeros(10,7), C_q'; zeros(7,10), zeros(7,7), C_s'; C_q, C_s, zeros(16,16)]; b = [Qg + Qk + Qv + Qroz; zeros(7,1); -dC_q*dq-dC_s*ds-dC_t]; x = A\b; ddq = x(1:10); dds = x(11:17); lambda = -x(18:33); dy = [dq; ds; ddq; dds]; A.5 Archivo para el cálculo de la fuerza generalizada de rozamiento “F_roz.m” function QS = F_roz(q,dq,s,N) global Param mu = Param.mu; vlim = Param.vlim; phi3 = q(7); tet3 = q(8); psi3 = q(9); Aphi3 = [1, 0, 0; 0, cos(phi3), -sin(phi3); 0, sin(phi3), cos(phi3)]; Atet3 = [cos(tet3), 0, sin(tet3); 0, 1, 0; -sin(tet3), 0, cos(tet3)]; Apsi3 = [cos(psi3), -sin(psi3), 0; sin(psi3), cos(psi3), 0; 0, 0, 1]; A4 = Apsi3*Atet3*Aphi3; H = jacobS(q,s); vS = H*dq; dir3 = A4*[0,1,0]'/norm(A4*[0,1,0]'); 51 Modelado y simulación de un vehículo auscultador para una vía ferroviaria a escala. vdesS = (dir3'*vS)*dir3; if norm(vdesS) < vlim FtS = -mu*N*vdesS/vlim; else FtS = -mu*N*vdesS/norm(vdesS); end QS = H'*FtS; end A.6 Archivo para el post-procesado de los resultados “PostProceso.m” function [ddq, dds, lambda] = PostProceso(t,q,dq,s,ds) ddq = zeros(length(t),10); dds = zeros(length(t),7); lambda = zeros(length(t),16); N = 15; for i = 1:length(t) MM = mat_MM; Qg = F_gravGeneral(q(i,:)); Qk = F_muelleGeneral(q(i,:)); Qroz = F_roz(q(i,:),dq(i,:)',s(i,:),N); QQv = F_inerc_NO_General(q(i,:),dq(i,:)); L = mat_L(q(i,:)); l = mat_dL(q(i,:),dq(i,:)); M = L'*MM*L; Qv = L'*(QQv - MM*l*dq(i,:)'); C_q = Cq(q(i,:),s(i,:)); dC_q = dCq(q(i,:),dq(i,:),s(i,:),ds(i,:)); dC_t = dCt(t(i)); C_s = Cs(q(i,:),s(i,:)); dC_s = dCs(q(i,:),dq(i,:),s(i,:),ds(i,:)); A = [M, zeros(10,7), C_q'; zeros(7,10), zeros(7,7), C_s'; C_q, C_s, zeros(16,16)]; b = [Qg + Qk + Qv + Qroz; zeros(7,1); -dC_q*dq(i,:)'-dC_s*ds(i,:)'-dC_t]; x = A\b; ddq(i,:) = x(1:10)'; dds(i,:) = x(11:17)'; lambda(i,:) = -x(18:33)'; end