Análisis de la influencia del ciclista en el comportamiento dinámico de la bicicleta

Equation Chapter 1 Section 1 Trabajo Fin de Grado Grado en Ingeniería en Tecnologías Industriales Análisis de la influencia del ciclista en el comportamiento dinámico de la bicicleta Autor: Francisco José Begines Gavira Tutor: Joaquín Ojeda Granja Dpto. Ingeniería Mecánica y Fabricación Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2019 iii Trabajo Fin de Grado Grado en Ingeniería en Tecnologías Industriales Análisis de la influencia del ciclista en el comportamiento dinámico de la bicicleta Autor: Francisco José Begines Gavira Tutor: Joaquín Ojeda Granja Profesor Contratado Doctor Interino Dpto. de Ingeniería Mecánica y Fabricación Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2019 v Trabajo Fin de Grado: Análisis de la influencia del ciclista en el comportamiento dinámico de la bicicleta El tribunal nombrado para juzgar el Proyecto arriba indicado, compuesto por los siguientes miembros: Presidente: Vocales: Secretario: Acuerdan otorgarle la calificación de: Sevilla, 2019 El Secretario del Tribunal Autor: Francisco José Begines Gavira Tutor: Joaquín Ojeda Granja vii A mi familia y maestros ix RESUMEN Este trabajo tiene objetivo principal el análisis de la influencia del comportamiento del ciclista en la dinámica de la bicicleta. Para ello se hará uso de la dinámica de sistemas multicuerpo. La dinámica directa permite resolver mediante ecuaciones diferenciales obtener la respuesta de un sistema mecánico en función de las acciones que éste sufra como puede ser un giro en el manillar o el pedaleo del ciclista. Motivo por el cuál la dinámica directa es de gran utilidad ya que permite preveer cualquier hipotético caso. Para plantear el sistema de ecuaciones diferenciales, se utilizarán las ecuaciones de Newton-Euler, para las cuales será necesario el uso de un programa de cálculo numérico, en este caso MATLAB, en el cual se trabajará de forma simbólica. Una vez se tenga resuelta la dinámica, se podrá evaluar el comportamiento de la bicicleta en función de la velocidad inicial de la bicicleta o de las acciones del ciclista (pedaleo, vuelco y giro del manillar) Este estudio se centrará en el comportamiento de la bicicleta cuando está girando. Para ver la influencia del ciclista se analizarán diferentes casos como son, variar la masa del ciclista, variar la inercia y por último variar la posición del centro de gravedad. Los cuales no serán determinantes en el comportamiento de la bicicleta, tal y como se muestra en este trabajo. Figura 15d Evolución temporal de 𝒚𝒄 31 Figura 15e Evolución temporal de 𝝋 32 Figura 15f Evolución temporal de 𝜺 32 Figura 15g Evolución temporal de 32 Figura 15h Evolución temporal de 32 Figura 16a Evolución temporal de 𝜸 33 Figura 16b Evolución temporal de 𝝍 33 Figura 16c Evolución temporal de 𝒙𝒄 33 Figura 16d Evolución temporal de 𝒚𝒄 33 Figura 16e Evolución temporal de 𝝋 33 Figura 16f Evolución temporal de 𝜺 33 Figura 16g Evolución temporal de 34 Figura 16h Evolución temporal de 34 𝜸󰇗 𝜸󰇗 󰇗 𝝍 󰇗 𝝍 xvii NOTACIÓN Ai Matriz i sen Función seno cos Función coseno q Coordenadas generalizadas 𝑞󰇗 Velocidades generalizadas 𝑞󰇘 Aceleraciones generalizadas 𝐴𝑇 Matriz traspuesta 𝛼 Aceleración angular 𝑖 Matriz/vector en coordenadas locales 𝜕𝑦/𝜕𝑥 Derivada de y con respecto a x 𝐶(𝑝,𝑝󰇗) Matriz/vector C dependiente de 𝑝 𝑦 𝑝󰇗 Notación 18 1 1 INTRODUCCIÓN a bicicleta es uno de los vehículos más estudiados a lo largo de la historia, pero la mayoría de los estudios tratan sobre el comportamiento de esta en una recta o incluso de su comportamiento a bajas velocidades, pero es difícil encontrar estudios sobre la influencia del ciclista mientras la bicicleta se encuentra en una curva. Motivo por el cual se ha decidido realizar un análisis sobre la influencia de la influencia del cicilista, ya sea por su complexión física (masa y posición del centro de gravedad), como por sus acciones (pedaleo, vuelco y giro del manillar). Para el análisis se hará uso de la dinámica de sistemas multicuerpos. Para resolver este tipo de problemas hay varios métodos, pero existen dos metodologías que destacan entre todas, la primera de ellas consiste en eliminar las fuerzas de reacción de las ecuaciones de la dinámica, expresando estas ecuaciones en términos de los grados de libertad, este método utiliza las variables asociadas a las articulaciones como grados de libertad para poder expresar el sistema de ecuaciones en función de estos grados de libertad. El uso de las variables asociadas a las articulaciones tiene la gran ventaja de reducir el número de coordenadas, y por lo tanto ecuaciones para poder resolver el sistema. Sin embargo, tiene la desventaja de tener unas ecuaciones no lineales y complejas que dificultarán su integración, a esta técnica se le conoce como embedding techniques. La otra metodología es la denominada como formulación aumentada, en esta formulación, las ecuaciones están expresadas en coordenadas en las que no todas son independientes. Este planteamiento tiene como desventajas la aparición de más coordenadas y ecuaciones además de que la resolución de estas ecuaciones requiere un algoritmo de resolución complejo. Sin embargo, tiene la ventaja de que el cálculo matricial tiene poca complejidad, además los programas basados en esta metodología suelen ser más fáciles de manejar. Una vez se ha expuesto las dos metodologías más extendidas en la resolución de problemas de dinámica de sistemas multicuerpo, se ha decidido desarrollar la segunda opción, por los motivos explicados anteriormente, aunque se va a utilizar un número mínimo de coordenadas. Dentro de esta metodología hay también una gran variedad de formas de resolver el problema dinámico, como son las ecuaciones de Newton-Euler, ya que es uno de los métodos más utilizados y fiables que existen. Este método tiene la gran desventaja de que aparecen los multiplicadores de Lagrange y sus ecuaciones tienen un alto grado de no-linealidad, pero como ventaja tiene que es un método con una formulación simple, además de utilizar un número mínimo de coordenadas. Hay métodos alternativos que se basan en, a través de las restricciones, intentar separar los multiplicadores de Lagrange, evitando que así aparezcan en la ecuación final de la dinámica, pero es necesario un trabajo matricial extenso, por estos motivos es por lo que se ha decidido utilizar las ecuaciones de Newton-Euler. Para resolver el problema será necesario un modelo dinámico que se comporte de la manera más aproximada posible a una bicicleta y la identificación de las restricciones necesarias para definir el comportamiento. La formulación de las restricciones y modelado de la bicicleta están basados en documentos sobre el comportamiento de la bicicleta del profesor Jose Luis Escalona, concretamente en los documentos [1] y [2]. El desarrollo de la cinemática se encuentra en el capítulo 2. Una vez se ha identificado el problema cinemático, se resolverá la dinámica con las ecuaciones de Newton- Euler, cuyas ecuaciones se desarrollarán en el capítulo 3, para ello será necesario un programa de cálculo numérico, en este caso MATLAB, y definir sus restricciones de forma simbólica para poder hacer posteriormente la integración de las ecuaciones diferenciales. Todo lo neceario para la resolución del problema mediante el programa de cálculo se explicará con más detenimiento en el capítulo 4. Cuando ya se pueda ejecutar el programa, será la hora de simular en diferentes casos el comportamiento dinámico de la bicicleta cuando se varían las propiedades del ciclista, todos los resultados serán analizados en el capítulo 5, realizando una conclusión final de estos en el capítulo 6. L Introducción 2 3 2. CINEMÁTICA DE LA BICICLETA ara el análisis de la influencia del ciclista en la dinámica de la bicicleta se hará uso de la dinámica de sistemas multicuerpos. El principal objetivo de la dinámica de sistemas multicuerpos es resolver la cinemática y dinámica de la bicicleta. Para será necesario la determinación de un modelo que simula el movimiento de una bicicleta. A partir del cual identificar las restricciones que definen el problema. 2.1 Modelo de la bicicleta Se definen 4 sólidos como partes integrantes del modelo: - Ruedas: sólidos 2 y 5, siendo el 2 la rueda trasera y el 5 la delantera. - Cuadro de la bicicleta: sólido 4. - Manillar: sólido 3 En cuanto a los sistemas de coordenadas, se han considerado los siguientes: - <X Y Z> será el sistema de referencia global. - < xi1 yi1 zi1> y <𝑥𝑖2 𝑦𝑖2 𝑧𝑖2> serán sistemas auxiliares - Los sistemas <xi yi zi> para i=2,3,4,5 serán los sistemas de referencia local situados en el centro de gravedad de cada sólido Figura 1. Modelo de la bicicleta [2] J.L. Escalona . En la figura 1, se puede apreciar tanto los sólidos definidos, como los sistemas de referencia. P Cinemática de la bicicleta 4 Para el estudio se definen las siguientes coordenadas: -xc, yc definen la posición del punto de contacto (C) de la rueda trasera en el plano X Y -𝜑 (phi) es el ángulo de guiñada, es decir, el ángulo que forma el eje xi1 del sistema intermedio de referencia con el eje X -𝜃 (theta) es el ángulo de vuelco, que se forma entre el eje 𝑧𝑖2 del segundo sistema intermedio con el eje 𝑧𝑖1 del primer sistema intermedio Con estos dos ángulos se puede determinar la orientación del plano que contiene el cuadro de la bicicleta. -𝜓 (psi) es el ángulo de rodadura, que forma el eje z2 del sistema asociado a la rueda trasera con el eje 𝑧𝑖2 del segundo sistema intermedio -𝛽 (beta) es el ángulo que forma el eje z3 del sistema del marco con el eje 𝑧𝑖2 del segundo sistema intermedio. -𝛾 (gamma) es el ángulo de dirección que forma el eje x4 del sistema de referencia del manillar con el eje x3 del sistema de referencia del marco. - 𝜀 (épsilon), ángulo de rodadura, es el ángulo que forma el eje z5 del sistema de referencia asociado a la rueda delantera con el eje z4 del manillar. Una vez definido estos parámetros se define el siguiente vector de coordenadas, al que se va a denominar como q: Conociendo todas las componentes de este vector se podrá calcular la posición y orientación de cualquier sólido de la bicicleta. 2.2 Orientación de sólidos Para conocer la orientación de cada sólido es necesario conocer las matrices de giro de éste, que se definen a continuación 𝐴𝑝ℎ𝑖=(cos(𝜑)−𝑠𝑒𝑛(𝜑)0 𝑠𝑒𝑛(𝜑)cos(𝜑)0 0 0 1) (2.2) 𝐴𝑡ℎ𝑒𝑡𝑎=(1 0 0 0 cos(𝜃)−sin(𝜃) 0sin(𝜃)cos(𝜃)) (2.3) 𝐴𝑝𝑠𝑖=(cos(𝜓)0 𝑠𝑒𝑛(𝜓) 0 1 0 −𝑠𝑒𝑛(𝜓)0 cos(𝜓)) (2.4) 𝐴𝑏𝑒𝑡𝑎=(cos(𝛽)0 𝑠𝑒𝑛(𝛽) 0 1 0 −𝑠𝑒𝑛(𝛽)0 cos(𝛽)) (2.5) 𝑞=[𝑥𝑐 𝑦𝑐 𝜑 𝜃 𝜓 𝛽 𝛾 𝜀]𝑇 (2.1) 5 5 Análisis de la influencia del ciclista en el comportamiento dinámico de la bicicleta 𝐴𝑔𝑎𝑚𝑚𝑎=(cos(𝛾)−𝑠𝑒𝑛(𝛾)0 𝑠𝑒𝑛(𝛾)cos(𝛾)0 0 0 1) (2.6) 𝐴𝑒𝑝𝑠𝑖𝑙𝑜𝑛=(cos(𝜀)0 𝑠𝑒𝑛(𝜀) 0 1 0 −𝑠𝑒𝑛(𝜀)0 cos(𝜀)) (2.7) Una vez ya han sido obtenidas estas matrices de giro de los ángulos, el siguiente paso será multiplicar las matrices necesarias para llegar a definir la orientación del sólido, en el caso de la bicicleta serían los siguientes: 𝐴𝑖1=𝐴𝑝ℎ𝑖 (2.8) 𝐴𝑖2=𝐴𝑝ℎ𝑖∗𝐴𝑡ℎ𝑒𝑡𝑎 (2.9) 𝐴2=𝐴𝑝ℎ𝑖∗𝐴𝑡ℎ𝑒𝑡𝑎∗𝐴𝑝𝑠𝑖 (2.10) 𝐴3=𝐴𝑝ℎ𝑖∗𝐴𝑡ℎ𝑒𝑡𝑎∗𝐴𝑏𝑒𝑡𝑎 (2.11) 𝐴4=𝐴𝑝ℎ𝑖∗𝐴𝑡ℎ𝑒𝑡𝑎∗𝐴𝑏𝑒𝑡𝑎∗A𝑔𝑎𝑚𝑚𝑎 (2.12) 𝐴5=𝐴𝑝ℎ𝑖∗𝐴𝑡ℎ𝑒𝑡𝑎∗𝐴𝑏𝑒𝑡𝑎∗A𝑔𝑎𝑚𝑚𝑎∗𝐴𝑒𝑝𝑠𝑖𝑙𝑜𝑛 (2.13) Debido a que se va a utilizar un programa de cálculo numérico, las derivadas de las matrices de los sólidos son necesarias para poder obtener posteriormente las velocidades angulares. Para el cálculo de estas derivadas se utiliza la siguiente expresión: 𝐴𝑖󰇗=𝑑𝐴𝑖 𝑑𝑡=𝜕𝐴𝑖 𝜕𝑞∙𝜕𝑞 𝜕𝑡=𝜕𝐴𝑖 𝜕𝑞∙𝑞󰇗 󰇗 (2.14) Para obtener la expresión de las velocidades angulares es necesario conocer las ecuaciones fundamentales de la cinemática de los sólidos rígidos. 𝑉0 󰇍 󰇍 󰇍 =𝑉𝑃 󰇍 󰇍 󰇍 󰇍 +𝑤 󰇍 󰇍 ∧𝑂𝑃 󰇍 󰇍 󰇍 󰇍 󰇍 (2.15) 𝑅𝑖𝑃=𝑅𝑖+𝐴𝑖𝑢𝑖𝑃 (2.16) Cinemática de la bicicleta 6 Si deriva con respecto al tiempo la ec 2.16: 𝑅𝑖𝑃 󰇗=𝑅𝑖󰇗+𝐴𝑖󰇗𝑢𝑖𝑃 (2.17) Comparando esta expresión con la velocidad relativa del punto P con respecto al O dado por la ec 2.15 se deduce que: 𝑤 󰇍 󰇍 ∧𝑢𝑖𝑃=𝐴𝑖󰇗𝑢𝑖𝑃 2.18) Donde a 𝑂𝑃 󰇍 󰇍 󰇍 󰇍 󰇍 se le ha denominado como 𝑢𝑖𝑃. Para el primer término se puede agrupar quitando el producto vectorial, de manera que el vector 𝑤 󰇍 󰇍 se convierta en la matriz antisimétrica 𝑤, además si 𝑢𝑖𝑃 se pasa el vector se obtiene a partir del vector en locales se obtiene la siguiente expresión 𝑤𝐴𝑖𝑇𝑢𝑖𝑃=𝐴𝑖󰇗𝑢𝑖𝑃 (2.19) Si se despeja de la ec 2.19, se obtiene la expresión de la velocidad angular antisimétrica Una vez obtenida la velocidad angular antisimétrica, se tendrá una matriz de la cual se podrá obtener los términos de la velocidad angular, ya que 𝑤 tiene la siguiente forma. 𝑤 =( 0 −𝑤𝑧𝑤𝑦 𝑤𝑧0 −𝑤𝑥 −𝑤𝑦𝑤𝑥0) (2.21) Así el vector de velocidad angular queda definido como 𝑤=[𝑤𝑥 𝑤𝑦 𝑤𝑧] La velocidad angular también se puede expresar como 𝑤=𝐺∙𝑞󰇗, donde: 𝐺𝑖=𝜕𝑤 𝜕𝑞󰇗 Esta matriz G,que es el jacobiano de w con respecto a q, será necesaria en cálculos posteriores. (2.22) 2.3 Posición, velocidad y aceleración del centro de gravedad de los sólidos Las posiciones del centro de gravedad de los sólidos de la bicicleta se obtienen de la siguiente forma: 𝑟𝑔2=𝑟𝑐+𝐴𝑖2∙𝑟𝑔2,𝑖2 (2.23) 𝑟𝑔3=𝑟𝑔2+𝐴3∙𝑟𝑔3,𝑖2 (2.24) 𝑤 =𝐴𝑖󰇗∙𝐴𝑇 (2.20) 7 7 Análisis de la influencia del ciclista en el comportamiento dinámico de la bicicleta 𝑟𝑔4=𝑟𝑔2+𝐴3∙𝑟𝑔4,𝑖2 (2.25) 𝑟𝑔5=𝑟𝑔4+𝐴4∙𝑟𝑔5,4 (2.26) Siendo 𝑟𝑐 el vector [xc yc 0]𝑇, 𝐴𝑖2 la matriz definida previamente y 𝑟𝑔2𝑖2=[0 0 𝑅𝑡]𝑇, es el vector en coordenadas locales del centro de gravedad de la rueda trasera con radio 𝑅𝑡, 𝑟𝑔3𝑖2=[𝑥3 0 𝑧3] y 𝑟𝑔4𝑖2= [𝑥4 0 𝑧4] representan las componentes de los vectores que van desde G2 hasta G3 y G4 con respecto al sistema de referencia del cuadro, respectivamente, y 𝑟𝑔5,4=[𝑥5 0 𝑧5] representa las componentes del vector que va desde G4 hasta G5 en el sistema de referencia 4 del manillar. Todos estos vectores se pueden ver en la figura 1. Se observa que se han seleccionado vectores cuyas componentes son constantes durante el movimiento de la bicicleta para que sus derivadas resulten nulas en el cálculo de las velocidades. Una vez que se tiene la posición de estos puntos, se procede a calcular la velocidad: 𝑣𝑔𝑖=𝑑𝑣𝑔𝑖 𝑑𝑡 =𝜕𝑟𝑔𝑖 𝜕𝑝 𝑑𝑝 𝑑𝑡 (2.27) 𝑣𝑔𝑖=𝐻𝑖∗𝑝󰇗 𝐻𝑖=𝜕𝑟𝑔𝑖 𝜕𝑝 i=2,3,4 y 5 (2.28) Siendo H el jacobiano de los vectores de posición con respecto a p, el vector de coordenadas. Para las acelaraciones la forma de proceder es la misma, siendo la acelaración la derivada con respecto al tiempo de la velocidad 𝑎𝑔=𝜕𝐻 𝜕𝑡∗𝑝󰇗+𝐻∗𝜕𝑝󰇗 𝜕𝑡=𝐻𝑝󰇗 +ℎ𝑝󰇘 (2.29) donde h es la derivada con respecto al tiempo de H. Para las aceleraciones angulares: 𝛼=𝑑𝐺𝑖𝑝 𝑑𝑡 =𝐺𝑖𝑑𝑝󰇗 𝑑𝑡+𝜕(𝐺𝑖𝑝󰇗) 𝜕𝑝 𝜕𝑝 𝜕𝑡=𝐺𝑖𝑝󰇘+𝑔𝑖𝑝󰇗 donde: (2.30) 𝑔𝑖=𝜕(𝐺𝑖𝑝󰇗) 𝜕𝑝 (2.31) 2.4 Restricciones cinemáticas Como ya se ha introducido antes, el problema tendrá un determinado número de coordenadas independientes, a partir de las cuales se podrán obtener las demás, denominadas dependientes. Estas deben cumplir ciertas condiciones que se denominan restricciones cinemáticas. Las restricciones cinemáticas son ecuaciones algebraicas, en general no lineales, que deben satisfacer las coordenadas en todo instante. En el modelo de la bicicleta tenemos dos tipos de restricciones: 1. Restricciones de contacto. Garantizan que la rueda delantera tenga en todo momento un punto en contacto con el suelo. Para el caso de la bicicleta, por las coordenadas seleccionadas (pues el punto de contacto C no tiene componente en Z), está garantizado que la rueda trasera mantiene un punto de contacto con el suelo y por tanto no requiere restricciones. Dinámica de la bicicleta 14 Donde: 𝑀=𝐿𝑀 𝐿𝑇 𝑄=𝐿𝑄 𝑄𝑣=𝐿(𝑄𝑣−𝑀 𝑙𝑞󰇗) (3.6) 3.1 Fuerzas de reacción El vector de fuerzas generalizas Q se puede descomponer como la suma de las fuerzas aplicadas generalizadas 𝑄𝑎𝑝 y las fuerzas de reacción generalizadas 𝑄𝑟𝑒𝑎𝑐, de tal forma que Q se puede definir como: 𝑄=𝑄𝑟𝑒𝑎𝑐+𝑄𝑎𝑝 (3.7) En un sistema mecánico sujeto sólo a restricciones holónomas, las fuerzas de reacción se pueden obtener como: 𝑄𝑟𝑒𝑎𝑐=−𝐶𝑝𝑇𝜆 (3.8) Donde Cp es la matriz jacobiana de las ecuaciones de restricción con respecto al vector de coordenadas y 𝜆 es el vector de multiplicadores de Lagrange. En el caso de un sistema sujeto a restricciones holónomas y no holónomas las fuerzas de reacción se calculan como: 𝑄𝑟𝑒𝑎𝑐=−𝐷𝑇𝜆 (3.9) donde la matriz D agrupa la matriz jacobiana de las restricciones holónomas con respecto al vector de coordenadas y la matriz jacobiana de las restricciones no holónomas con respecto al vector de velocidades generalizadas. En el caso de la bicicleta se trata de la matriz de coeficientes de las ecuaciones de restricción a nivel de velocidad, como sigue: 𝐷=[𝐶𝑃𝑐𝑜𝑛 𝐵] (3.10) 3.2 Fuerzas aplicadas El vector de fuerzas aplicadas generalizadas Qap se obtiene como: 𝑄𝑎𝑝=𝑄𝑔𝑟𝑎𝑣+𝑄𝑎𝑐𝑐+𝑄𝑎𝑒𝑟𝑜 (3.11) donde 𝑄𝑔𝑟𝑎𝑣 es la fuerza generalizada debida a la gravedad, 𝑄𝑎𝑒𝑟𝑜 es la fuerza generalizada debida a la resistencia aerodinámica y 𝑄𝑎𝑐𝑐 es la fuerza generalizada debida a la acción del ciclista en la bicicleta durante su conducción. Para calcular la expresión de los vectores de fuerza generalizada se va a utilizar el principio de las potencias virtuales. Según este principio, el vector de fuerza generalizada es el vector que multiplicado escalarmente por el vector de velocidades virtuales 𝑝󰇗∗ da lugar a la potencia virtual que desarrolla la fuerza o par. A continuación, se desarrolla el calculo de los vectores. 15 Análisis de la influencia del ciclista en el comportamiento dinámico de la bicicleta 3.2.1 Fuerza generalizada de gravedad Los pesos de los sólidos móviles que forman la bicicleta son vectores aplicados en los centros de gravedad y que llevan dirección –Z, como sigue: 𝐹𝑔𝑟𝑎𝑣 𝑖=[00 −𝑚𝑖𝑔] ,𝑖=2,3,4 𝑦 5 (3.12) Por tanto, la potencia virtual que desarrollan los pesos se obtiene como: 𝑊󰇗𝑔𝑟𝑎𝑣 ∗=∑(𝑣𝑖𝑔 ∗)𝑇𝐹𝑔𝑟𝑎𝑣 𝑖 5 𝑖=2 (3.13) Donde se ha utilizado la ec 4.3 para el cálculo de la velocidad virtual del centro de gravedad de los sólidos. El vector Qgrav queda entonces: 𝑄𝑔𝑟𝑎𝑣=∑(𝐻𝑖)𝑇𝐹𝑔𝑟𝑎𝑣 𝑖 5 𝑖=2 (3.14) 3.2.2 Fuerza generalizada de acción del ciclista Para el manejo de la bicicleta se supone que el ciclista ejerce tres pares de control tal y como se puede ver en la figura 2. El par de dirección Mdir aplicado en dirección z4, el par de pedaleo Mped aplicado en dirección y2 y el par de vuelco Mvue aplicado en dirección x3. El valor de estos pares es en principio desconocido. Figura 3. Acciones del ciclista [3] J.L. Escalona La potencia virtual de estos pares se calcula de la forma siguiente: 𝑊󰇗𝑎𝑐𝑐=(𝑤4∗)𝑇[00 𝑀𝑑𝑖𝑟]+(𝑤2∗)𝑇[0 𝑀𝑝𝑒𝑑 0]+(𝑤3∗)𝑇[𝑀𝑣𝑢𝑒 0 0]= =(𝑝󰇗∗)[𝐺4𝑇[00 𝑀𝑑𝑖𝑟]+𝐺2𝑇[0 𝑀𝑝𝑒𝑑 0]+𝐺3𝑇[𝑀𝑣𝑢𝑒 00]] (3.15) Dinámica de la bicicleta 16 El vector Qaccion se puede calcular como: 𝑄𝑎𝑐𝑐=𝑀𝑑𝑖𝑟 𝐺4(3,:)𝑇+𝑀𝑝𝑒𝑑 𝐺2(2,:)𝑇+𝑀𝑣𝑢𝑒 𝐺 3(1,:)𝑇 (3.16) Donde Gi (j, :) representa la fila j de la matriz Gi. El vector Qacc se puede expresar de una forma más compacta como: 𝑄=[𝐺4(3,:)𝑇 𝐺2(2,:)𝑇 𝐺3(1,:)𝑇][𝑀𝑑𝑖𝑟 𝑀𝑝𝑒𝑑 𝑀𝑣𝑢𝑒] (3.17) donde la matriz S contiene las filas de las matrices Gi y el vector T= [Mdir Mped Mvue] contiene los tres pares de acción del ciclista. 3.2.3 Fuerza generalizada de resistencia aerodinámica Se va a aproximar la resistencia aerodinámica como proporcional al cuadrado de la velocidad de avance de la bicicleta y que está aplicada exclusivamente en el centro de gravedad del sólido 3, que forman el cuadro de la bicicleta y el ciclista, ya que al ser éste el sólido más voluminoso, ofrece la mayoría de la resistencia al aire. En ese caso, el módulo del vector de resistencia aerodinámica resulta: |𝐹𝑎𝑒𝑟𝑜3|=𝑐𝑣(𝑣3𝐺)𝑇𝑣3𝐺 (3.18) y su dirección será la contraria a la de la velocidad del centro de gravedad del sólido 3, por tanto: 𝐹𝑎𝑒𝑟𝑜3=−|𝐹𝑎𝑒𝑟𝑜3|∗𝑣3𝐺 |𝑣3𝐺|=−𝑐𝑣|𝑣3𝐺|𝑣3𝐺 (3.19) La potencia virtual que desarrolla la fuerza aerodinámica resulta: 𝑊󰇗𝑎𝑒𝑟𝑜 ∗=(𝑣3𝐺)𝑇𝐹𝑎𝑒𝑟𝑜3=(𝑝󰇗∗)𝑇𝐻3𝑇 𝐹𝑎𝑒𝑟𝑜3 (3.20) Por tanto, el vector Qaero queda entonces: 𝑄𝑎𝑒𝑟𝑜=𝐻3𝑇 𝐹𝑎𝑒𝑟𝑜3=−𝑐𝑣|𝑣3𝐺| 𝐻3𝑇 𝑣3𝐺 (3.21) Se ha supuesto un coeficiente aerodinámico de 0.5. Por tanto, las ecuaciones de Newton-Euler de la bicicleta se pueden escribir de la forma: 𝑀𝑝󰇘+𝐷𝑇𝜆=𝑄𝑎𝑝+𝑄𝑣 (3.22) Agrupando las ecuaciones de Newton-Euler con las ecuaciones de restricción a nivel de aceleración quedan las ecuaciones de movimiento de la bicicleta en forma de ecuaciones diferenciales-algebraicas (DAE) como sigue: [𝑀 𝐷𝑇 𝐷 0][𝑝󰇘𝜆]=[𝑄𝑎𝑝+𝑄𝑣 −𝐷󰇗𝑝󰇗 ] (3.23) 4 RESOLUCIÓN NUMÉRICA omo ya se ha comentado en reiteradas ocasiones, para la resolución del problema se va a acudir a un programa de cálculo numérico, MATLAB. Para poder resolver el problema dinámico hay que seguir un determinado orden, el cual se detalla en el siguiente esquema: Antes de proceder a explicar como el procedimento de resolución, resulta interesante realizar un inciso sobre las coordenadas y los grados de libertad. Tal y como se acaba de deducir el problema de la bicicleta consta de 9 coordenadas y 3 grados de libertad, es decir, 3 coordenadas a partir de las cuales se podría obtener los valores de las demás. A priori, se podrían elegir 3 variables al azar, pero la presencia de restricciones holónomas y no holónomas hace que el problema cinemático exija una elección de coordenadas a través de un método de cálculo, el cual facilite la elección de coordenadas, provocando con ello que el problema tenga una resolución más simple. Como se podría leer en los documentos de referencia [2] y [3], lo más recomendable es separar las variables en dinámicas, cinemáticas y dependientes. Las variables dinámicas serán las independientes, mientras que las coordenadas restantes, se dividen en cinemáticas, para las coordenadas asociadas a las variables con restricciones no holónomas y las dependientes asociadas a las restricciones holónomas. De esta manera que quedarían los siguientes vectores de variables 4.1 Simbólico Este apartado consiste en definir en simbólico las restricciones, y todo lo que sea necesario para definir estas, como son las matrices de giro, velocidades angulares y velocidades del centro de gravedad de los sólidos y también se definen las matrices necesarias para cálculos posteriores como son el jacobiano de las restricciones, la derivada de este jacobiano. Además, se definen las matrices que intervienen en las ecuaciones para resolver la dinámica, matriz de masas, matrices de restricciones y las matrices de fuerzas generalizadas. Una vez que se tenga este fichero, se compila creando unos archivos que contienen las matrices y restricciones anteriormente comentadas. Esta acción tiene un alto coste de computación por lo que también es alto el tiempo de compilación, pero solo será necesario realizarlo una vez, ya que una vez compilado, creará los archivos a partir de los cuales se trabaja. 4.2 Dinámica Cuando se tengan todos los ficheros creados por simbolicobici.m se creará otro con el nombre de dinámicadirecta. Este fichero tendrá como objetivo resolver las ecuaciones de la dinámica y a partir del cual C 𝑞𝑑𝑖𝑛=[𝜃 𝜓 𝛾] 𝑞𝑐𝑖𝑛=[𝑥𝑐 𝑦𝑐 𝜑 𝜀] 𝑞𝑑𝑒𝑝=[𝛽𝜉] (4.1) SIMBÓLICO Restricciones Matrices necesarias DINÁMICA Ec. Newton-Euler FICHEROS ADICIONALES Newton Raphson Transf. velocidad SIMULABICI Condiciones iniciales 17 Resolución numérica 18 poder integrar las ecuaciones, por lo que este fichero será una función, que tendrá como parámetros de entrada el tiempo y el vector de estado y, siendo este vector: 𝑦=[𝑞𝑖𝑛𝑑 𝑞󰇗𝑖𝑛𝑑 𝑞𝑐𝑖𝑛] (4.2) O desarrollando este vector: 𝑦= [ 𝜃𝛾 𝜓 𝜃󰇗𝛾󰇗𝜓󰇗 𝑥𝑐 𝑦𝑐 𝛽𝜀 ] (4.3) La salida de esta función será el vector dy, es decir la derivada con respecto al tiempo del vector definido en la ecuación 5.1, resultando el siguiente vector. 𝑑𝑦=𝑦󰇗=[𝑞󰇗𝑖𝑛𝑑 𝑞󰇘𝑖𝑛𝑑 𝑞󰇗𝑐𝑖𝑛] (4.4) Por la tanto, el fichero llamado dinamicadirecta.m tendrá como objetivo definir el vector dy. Para la correcta obtención de este vector hay que seguir un minucioso orden. El vector qind ya se tendría definido al tratarse de las variables independientes, cuyos valores serán definidos antes de iniciar el procedimiento. El primer paso será definir el vector de estado inicial, esto es, el vector “y” en el punto inicial, para esto solo es necesario declarar las variables de entrada como el vector y en el adecuado orden. Como se puede apreciar en la ec 5.2, no aparecen las variables que se asoció a las restricciones holónomas, beta y chi, por lo que habría que estimar su valor inicial. El siguiente paso sería obtener las coordenadas dependientes mediante el uso del algoritmo de Newton-Raphson. Este paso permitirá obtener las velocidades del vector 𝑞󰇗𝑑𝑒𝑝 mediante la matriz de transformación de velocidades, ya que para poder obtener esta matriz se necesitan el valor de las coordenadas dependientes. El vector 𝑞󰇗𝑑𝑒𝑝, el cual está compuesto de: 𝑞󰇗𝑑𝑒𝑝=[𝑞󰇗𝑐𝑖𝑛 𝑞󰇗𝑑𝑒𝑝] (4.5) Siendo el primer término (𝑞󰇗𝑐𝑖𝑛) el vector de velocidades correspondiente a las coordenadas cinemáticas y el segundo término es el asociado a las velocidades de las coordenadas dependientes, 𝑞󰇗𝑑𝑒𝑝 solo será necesario para poder definir las ecuaciones de Newton-Euler y por este motivo deben ser calculada El último paso es definir el vector de aceleraciones de las coordenadas independientes. Para ello se utilizarán, como ya se adelantó anteriormente, las ecuaciones de Newton-Euler. Estas ecuaciones necesitan depende de la matriz de masas, las matrices C, Cq y dCq, además de las fuerzas generalizadas, estas matrices ya estaban 19 Análisis de la influencia del ciclista en el comportamiento dinámico de la bicicleta calculadas en el fichero simbólico, por lo que solo habrá que insertar las parámetros y variables de las que dependan para que el programa devuelva los valores de las matrices y así poder formar el sistema de ecuaciones tal y como se muestran en las ecuaciones 4.17. Una vez se tenga todo definido lo único que queda es despejar el sistema para calcular el vector ddq, el cual contiene el vector de aceleraciones y las fuerzas de reacción, pero solo se necesitan las aceleraciones del vector de variables independientes. Como una vez que se termine este fichero, el objetivo es integrarlo, es conveniente añadir los parámetros de Baumgarten ya que estos métodos de resolución suelen tener imprecisiones y estos parámetros ayudan a que las imprecisiones no sean significantes o incluso que no aparezcan, algo que es bastante complejo que ocurra. La inclusión de estos parámetros en las ecuaciones de Newton Euler es de la siguiente manera: [𝑀 𝐷𝑇 𝐷 0][𝑝󰇘𝜆]=[𝑄𝑎𝑝+𝑄𝑣 −𝐷󰇗𝑝󰇗+𝐸𝑠𝑡𝐵𝑎𝑢𝑚] (4.6) Donde 𝐸𝑠𝑡𝐵𝑎𝑢𝑚 es: 𝐸𝑠𝑡𝐵𝑎𝑢𝑚=−2𝛼1(𝐷󰇗𝑝󰇗)−𝛼22𝐶 (4.7) Y con esto ya se puede dar por finalizado el archivo dinamicadirecta.m 4.3 Ficheros adicionales Para conseguir que el fichero donde están definidas las ecuaciones de Newton-Euler se necesitan dos ficheros adicionales que deben de ser creados por separados. 4.3.1 Newton Raphson Para el cálculo de las variables asociadas a las restricciones holónomas, se hace uso del algoritmo de Newton- Raphson, basado en un método iterativo en el que cuantas más veces se itera, más cerca se estará del valor exacto. Para este trabajo se han realizado 100 iteraciones. Una opción válida podría ser elegir una tolerancia mínima e iterar hasta que estemos dentro de esta tolerancia, pero el programa que se ha usado para el cálculo numérico a veces no permite tener la precisión que se requiere y puede dar lugar a error o también podría tardar un tiempo de compilación excesivo. Como se han realizado un número finito de iteraciones, para comprobar que el resultado de este algoritmo es correcto existe una verificación que consiste en evaluar en las ecuaciones de restricción los valores calculados, si el resultado es igual a cero, el valor calculado será el exacto. En este trabajo se ha conseguido un valor del orden de 10−16, por lo que se da por valido las 100 iteraciones sin necesidad de continuar el proceso de iteración. Este algoritmo se basa en la siguiente expresión: 𝑞𝑑𝑒𝑝𝑖+1=𝑞𝑑𝑒𝑝𝑖−(𝜕𝐶𝑑𝑒𝑝 𝜕𝑞 )−1∗𝐶𝑑𝑒𝑝 (4.8) Siendo 𝑞𝑑𝑒𝑝𝑖+1 el valor de la nueva iteración, 𝑞𝑑𝑒𝑝𝑖 el valor de la iteración anterior,la matriz que aparece en el tercer término que se encuentra es la inversa de la matriz resultante de hacer el jacobiano a las restricciones holónomas y después evaluadas con el valor actual de qdepi y por ultimo Cdep son las restricciones evaluadas en el punto de la iteración. Cabe destacar que para realizar este procedimiento es necesario un estimador inicial a partir del cual el resultado inicial podría variar. Además, si se hubieran puesto unas condiciones iniciales nulas, quedaría un jacobiano completo de ceros una vez evaluado en este punto, por lo que sería imposible calcular la inversa, y por tanto no se podría aplicar este método de resolución. Resolución numérica 20 En este caso, los valores iniciales son: 𝑞𝑑𝑒𝑝𝑖𝑛𝑖=[𝛽𝜉]=[0 𝜋2 ⁄] (4.9) 4.3.2 Matriz de transformación de velocidades Para calcular las velocidades de las coordenadas se utiliza la matriz de trasformación de velocidades, para llegar hasta ella es necesario definir los siguientes conceptos: Para empezar, se crea el vector de restricciones C: 𝐶(𝑝,𝑝󰇗)=[𝐶𝑐𝑜𝑛(𝑝) 𝐶𝑟𝑜𝑑(𝑝,𝑝󰇗)] (4.10) El objetivo es encontrar una forma lineal de obtener las velocidades: 𝐷∙𝑝󰇗=0 (4.11) Para las restricciones no holónomas: 𝐶󰇗(𝑝)=𝜕𝐶𝑐𝑜𝑛 𝜕𝑡 =𝜕𝐶𝑐𝑜𝑛 𝜕𝑝 𝜕𝑝 𝜕𝑡=𝜕𝐶𝑐𝑜𝑛 𝜕𝑝 ∙𝑝󰇗=𝐷∙𝑝󰇗 (4.12) Para las restricciones holónomas: 𝐶󰇗(𝑝,𝑝󰇗)=𝜕𝐶𝑟𝑜𝑑 𝜕𝑝󰇗∙𝑝󰇗=𝐷∙𝑝󰇗 (4.13) Agrupando queda: 𝐷= [ 𝜕𝐶𝑐𝑜𝑛 𝜕𝑝 𝜕𝐶𝑟𝑜𝑑 𝜕𝑝󰇗 ] (4.14) Separando entre variables independientes y dependientes: 𝐷𝑖𝑛𝑑∙𝑝󰇗𝑖𝑛𝑑+𝐷𝑑𝑒𝑝∙𝑝󰇗𝑑𝑒𝑝=0 (4.15) Ddep depende de los valores de posición de las variables dependientes, por lo que habrá que suponer unos valores iniciales para poder realizar los cálculos. Los valores iniciales son los siguientes: 𝑞𝑘𝑖𝑛0=[0 0 0 0] (4.16) 21 Análisis de la influencia del ciclista en el comportamiento dinámico de la bicicleta Cuando se tengan las matrices totalmente calculadas se despeja la ecuación: Aquí se tiene una ecuación lineal donde se puede despejar fácilmente el vector de velocidades p󰇗dep. 4.4 Simula bicicleta El último fichero que se creará se llama simulabici.m y a partir de éste se podrá obtener todos los diferentes casos que se quieran. Simulabici consistirá en integrar el archivo dinámicadirecta y poder obtener el vector de coordenadas en cada momento, lo que permitiría una infinidad de posibilidades como por ejemplo ver la evolución de las coordenadas con respecto al tiempo, también se puede crear un fichero (animaBici), que teniendo como entrada el vector de coordenadas en cada momento, muestra una animación del movimiento del ciclista, algo que resulta muy útil para ver la influencia de las condiciones iniciales o de la acción del ciclista. También se puede calcular las restricciones y tener una idea del error que tiene este método, o incluso si el error es demasiado grande ayudaría a saber que algo se está haciendo erróneamente. Para conseguir todo esto lo primero que hay que hacer es definir un intervalo de tiempo en el que se quiere que se integre, luego habrá que definir el vector de estado en el punto inicial. Ya con esto definido se puede resolver estas ecuaciones con un integrador de tipo ODE (Ordinary Differential Equations). Del cual se obtendrá el valor del vector y en cada momento. Hay que resaltar que solo con el vector no se podría realizar ningún análisis de los comentados anteriormente, ya que para tener la posición en cada instante se necesitaría las coordenadas dependientes y para las restricciones, las velocidades de la mayoría de las coordenadas. Para estos casos se actúa de la misma manera que en anteriores capítulos o incluso como el procedimiento del archivo de dinamicadirecta.m. Para las coordenadas dependientes se utilizará el método de Newton Raphson y para las velocidades el método de la matriz de transformación de velocidades. Una vez se han realizado todos los pasos anteriores ya se puede empezar a simular, pero para los casos en los que el momento de giro sea distinto de cero, habrá que insertar un controlador en el fichero de dinámica de manera que cuando la aceleración en theta, es decir, en el vuelco, sea mayor de un 5% de 2𝜋 el valor de ésta será 0. De esta manera se asegurará el comportamiento correcto cuando la bicicleta gire. 𝐷𝑑𝑒𝑝∗𝑝󰇗𝑑𝑒𝑝=−𝐷𝑖𝑛𝑑∗𝑝󰇗𝑖𝑛𝑑 (4.17) Resolución numérica 22 23 5 RESULTADOS n este capítulo se analizan todos los casos que se han tratado resultantes de asignar distintos valores a los parámetros de entrada con el archivo de dinámica directa, para ello se puede ir cambiando los valores iniciales del vector y o cambiar el valor de los momentos aplicados en la bicicleta, incluso se puede cambiar ambas opciones a la vez. Para el primer caso se ha impuesto una condición inicial en 𝜓󰇗 de valor 20 rad/s, ya que teniendo en cuenta que el radio de la rueda trasera es de 0,3 m, al introducir esa condición inicial, se está imponiendo una velocidad de avance en la rueda trasera de unos 6 m/s. Además, hay que conocer la masa del ciclista (incorporada en el sólido 3), dato que se define en el archivo Parámetros que se puede encontrar en el anexo, para esta simulación se ha definido una masa de 75 kg. Al tener definido ya el vector de estado en el instante inicial como todo cero menos el parámetro 𝜉, cuyo valor era 𝜋 2, la bicicleta queda en una posición recta, tal y como se aprecia en la figura 4. Figura 4. Posición inicial de la bicicleta Con estas condiciones la bicicleta comienza a avanzar y va perdiendo velocidad hasta que termina parándose, si ahora además de la velocidad inicial en psi, se aplica un momento en el pedal, se obtendría un movimiento rectilíneo continuo. Para esta simulación se han definido unos parámetros de Baumgarten de 20 cada uno, ya que de no poner estos parámetros la inestabilidad del sistema provocaría que beta, parámetro que define la inclinación, varíe, provocando que en ciertos puntos no se cumplan las restricciones de contacto y que la rueda delantera pierda el contacto con el suelo. Realizando un breve análisis de la dinámica de este sistema cuando está en línea recta, se obtienen los siguientes resultados: En la figura 5a se puede garantizar que el problema se ha resuelto correctamente, ya que esta figura muestra el valor de las restricciones en cada instante, se recuerda que para que las restricciones se cumplan sus valores tienen que ser cero en todo momento. Una vez se ha explicado esto, es el momento ideal para ver el efecto que tiene los parámetros de Baumgarten en la estabilidad de la integración. Para ello basta con asignar un valor nulo a los parámetros y obtener los resultados. E Resultados 30 Figura 14g Evolución temporal de 𝜸󰇗 Figura 14 Evolución temporal de 𝝍󰇗 Al contrario que en el apartado anterior, no se aprecia una diferencia importante en la velocidad de psi ni en el avance en el eje x, pero sí en el eje y, cuya diferencia entre los casos extremos varía algo más de 4 m. Por lo que se podría decir que variar la inercia es más crítico a la hora de tomar una curva que variar la masa, ya que esta variación afecta en mayor medida al avance en el eje y. 7.3 Comportamiento al variar la posición del centro de gravedad Hasta ahora, la figura del ciclista se había dejado a un lado, ya que solo se había analizado la influencia de la masa, y la de la inercia, que aunque depende de la masa y de la geometría, la masa adquiere más importancia en su valor, aunque la geometría sea cuadráticamente proporcional a la inercia, es díficil que la posición de el centro de gravedad varíe mucho, en cambio en cuanto a la masa, aunque es directacmente proporcional a la inercia, es más frecuente que el valor de ésta varíe mucho. En cambio al variar la posición del centro de gravedad, el análisis se centra solo en la figura del ciclista, ya sea en la posición en la que esté pedaleando, como la geometría del tren superior de éste. Este apartado se dividirá en dos, en el primero se comparará se aumentará la posición del centro de gravedad para una misma posición, y la segunda consistirá en variar el ángulo del ciclista a dos posiciones más erguidas. Para cambiar la posición del centro de gravedad se añade un sólido 6 (solidario al sólido 3), que será el tren superior del ciclista, con una posición fija. Este sólido modifica la posición del centro de gravedad de la siguiente forma: 𝐺𝑥=0∗ 𝑚3 𝑚3+𝑚6+𝑥𝑠∗𝑚6 𝑚3+𝑚6 (5.3) Donde 𝐺𝑥 es la distancia en el eje x entre la coordenada del centro de gravedad de 3 y el centro de gravedad del conjunto y 𝑥𝑠 es la coordenada del centro de gravedad del ciclista. Una vez se ha obtenido 𝐺𝑥 se calcula la nueva posición del centro de gravedad 3: 𝑥3′=𝑥3+𝐺𝑥 (5.4) Para el eje z se procede de la misma manera. Pero al desplazar el centro de gravedad también hay que desplazar el valor de las inercias, para ello se utiliza el teorema de Steiner. h 31 Análisis de la influencia del ciclista en el comportamiento dinámico de la bicicleta 𝐼𝑂=𝐼𝐺+𝑚𝑑2 (5.5) Siendo 𝐼𝑂la inercia en el punto nuevo, 𝐼𝐺 la inercia en el centro de gravedad antiguo, m la masa del sólido y d la distancia entre los dos puntos (O y G). Cuando se obtenga la matriz de inercias en el punto nuevo tanto del sólido 3 como del sólido 6, la matriz inercia nueva resultará de la siguiente forma: 𝐼3𝑂′=𝐼3𝑂+𝐼6𝑂 (5.6) Donde 𝐼3𝑂 es la matriz de inercias del sólido 3 en el punto nuevo, al igual que 𝐼6𝑂 para el sólido 6. 7.3.1 Comportamiento al cambiar el centro de gravedad del ciclista. En este subapartado se va a comparar la posición normal, la posición del centro de gravedad del conjunto bicicleta ciclista y la posición si el centro de gravedad estuviera 10 cm más arriba en el eje z. Figura 15a Evolución temporal de 𝜸 Figura 15b Evolución temporal de 𝝍 Figura 15c Evolución temporal de 𝒙𝒄 Figura 15d Evolución temporal de 𝒚𝒄 Resultados 32 En este apartado prácticamente se observan los mismos resultados que en el anterior, ya que se ha variado la inercia también al modificar el centro de gravedad. Pero hay un detalle que llama la atención en la figuras f y g. Como se aprecia en esta s figuras, el giro es menor en el caso en el que se ha calculado el nuevo centro de gravedad que cuando el centro estaba en el sillín, pero si se sigue subiendo la posición del centro de gravedad, solo en el eje z, este giro aumenta en vez de disminuir. 7.3.2 Comportamiento al variar la posición del ciclista. Para este subapartado se puede proceder de varias maneras, como por ejemplo, definir un ángulo de inclinación del ciclista y a partir de éste analizar su influencia asignando distintos valores a este ángulo. Pero se puede hacer una aproximación sin tener que definir este ángulo y consiste en variar el centro de gravedad si la posición es más erguida o bajarlo si la posición es más horinzontal, además si la posición es más vertical será necesario aumentar el coeficiente de resistencia aerodinámica, ya que el área afectada por el viento aumenta. Por lo que para este análisis solo se va a aumentar el coeficiente de resistencia aerodinámica, ya que se acaba de demostrar que un cambio en el centro de gravedad apenas afectará a los resultados. Para las dos simulaciones se tomarán un coeficiente de 0.5 y 0.8 obteniendo los siguientes resultados: Figura 15e Evolución temporal de 𝝋 Figura 15 Evolución temporal de 𝜺 Figura 15g Evolución temporal de 𝜸󰇗 Figura 15 Evolución temporal de 𝝍󰇗 h f 33 Análisis de la influencia del ciclista en el comportamiento dinámico de la bicicleta Figura 16a Evolución temporal de 𝜸 Figura 16b Evolución temporal de 𝝍 Figura 16c Evolución temporal de 𝒙𝒄 Figura 16d Evolución temporal de 𝒚𝒄 Figura 16e Evolución temporal de 𝝋 Figura 16 Evolución temporal de 𝜺 f Resultados 34 Figura 16g Evolución temporal de 𝜸󰇗 Figura 16 Evolución temporal de 𝝍󰇗 A difencia del subapartado anterior, en estas simulaciones sí se pueden apreciar leves diferencias en prácticamente todas las gráficas, como por ejemplo la diferencia de la velocidad es de unos 3 rad/s, siendo en el caso de cv=0,5 la que tiene mayor velocidad, también se aprecia que en el eje y hay un desplazamiento de más de 3 m para el caso con menos cv. Unos resultados que era obvio, ya que cuanto mayor sea el coeficiente mayor será la fuerza que tendrá que vencer y por lo tanto no avanzará tanto. A priori parece que los resultados no son lo suficientes diferentes como para decir que el coeficiente aerodinámico no es muy importante en la dinámica del ciclista, pero se recuerda que esta simulación se ha tomado una velocidad inicial de aproximadamente de 10 km/h y la fuerza de resistencia aerodinámica es proporcional cuadráticamente, por lo que a una velocidad de 30 km/h el coeficiente aerodinámico es un factor clave en la dinámica del ciclista. h 6 CONCLUSIONES Tras analizar las diferentes situaciones se podría afirmar que para la dinámica de la bicicleta cuando está girando los factores a tener más en cuenta son la inercia y la fuerza de resistencia aerodinámica, ya que estas dos son las únicas de provocar grandes diferencias de resultados, mientras la masa, aunque se aprecian diferencias notables, la diferencia del valor de la masa para que estas se puedan apreciar tiene que ser muy elevada. Por último, destacar la influencia del centro de gravedad que, aunque parte de su influencia se debe al cambio en los valores de la matriz de inercia, si solo se aumenta en el z la posición mayor cuando se había visto que al aumentar la posición el giro era menor. Para próximos estudios resultaría interesante el estudio de otro método de resolución en el cual no haga falta la implementación de un controlador para el correcto funcionamiento de la dinámica de la bicicleta en una curva, ya que en este trabajo se le ha impedido el vuelco. Además, también se podría realizar una simulación donde el ciclista tuviera movilidad, como puede ser mover los brazos o inclinarse en las curvas y así contrarrestrar los efectos de la fuerza centrífuga. 35 Conclusiones 36 7 CÓDIGO MATLAB En el código a la variable chi se le ha llamado e2 7.1 Simbolicobici.m syms xc yc phi theta psi beta gamma epsilon e2 real syms dxc dyc dphi dtheta dpsi dbeta dgamma depsilon de2 real syms ddxc ddyc ddphi ddtheta ddpsi ddbeta ddgamma ddepsilon dde2 real syms x3 z3 x4 z4 x5 z5 Rt Rd cv real syms m2 m3 m4 m5 m6 I2xx I2yy I2zz I3xx I3yy I3zz I3xz real syms I4xx I4yy I4zz I4xz I5xx I5yy I5zz real syms Mdir Mpedal Mvuel real q=[xc yc phi theta psi beta gamma epsilon ]; dq=[dxc dyc dphi dtheta dpsi dbeta dgamma depsilon ]; ddq=[ddxc ddyc ddphi ddtheta ddpsi ddbeta ddgamma ddepsilon ]; %DEFINICION MATRICES DE GIRO Atheta=[ 1 0 0 0 cos(theta) -sin(theta) 0 sin(theta) cos(theta)]; Agamma=[ cos(gamma) -sin(gamma) 0 sin(gamma) cos(gamma) 0 0 0 1]; Apsi= [ cos(psi) 0 sin(psi) 0 1 0 -sin(psi) 0 cos(psi)]; Abeta= [ cos(beta) 0 sin(beta) 0 1 0 -sin(beta) 0 cos(beta)]; Aepsilon=[ cos(epsilon) 0 sin(epsilon) 0 1 0 -sin(epsilon) 0 cos(epsilon)]; Aphi= [ cos(phi) -sin(phi) 0 sin(phi) cos(phi) 0 0 0 1]; %%MATRICES DE SOLIDOS At1=Aphi; At2=Aphi*Atheta; A2=Aphi*Atheta*Apsi; A3=Aphi*Atheta*Abeta; A4=A3*Agamma; A5=A4*Aepsilon; 37 Código matlab 38 % Para definir velocidades angulares diff(A)*A'=antisimetrico(w) dAt2=diff(At2,phi)*dphi+diff(At2,theta)*dtheta; dA2=diff(A2,phi)*dphi+diff(A2,theta)*dtheta+diff(A2,psi)*dpsi; dA3=diff(A3,phi)*dphi+diff(A3,theta)*dtheta+diff(A3,beta)*dbeta; dA4=diff(A4,phi)*dphi+diff(A4,theta)*dtheta+diff(A4,beta)*dbeta... +diff(A4,gamma)*dgamma; dA5=diff(A5,phi)*dphi+diff(A5,theta)*dtheta+diff(A5,beta)*dbeta... +diff(A5,gamma)*dgamma+diff(A5,epsilon)*depsilon; dAt2=simplify(dAt2); dA2=simplify(dA2); dA3=simplify(dA3); dA4=simplify(dA4); dA5=simplify(dA5); %VELOCIDADES ANGULARES w2_sk=simplify(dA2*A2'); w3_sk=simplify(dA3*A3'); w4_sk=simplify(dA4*A4'); w5_sk=simplify(dA5*A5'); w2=[w2_sk(3,2) w2_sk(1,3) w2_sk(2,1)]'; w3=[w3_sk(3,2) w3_sk(1,3) w3_sk(2,1)]'; w4=[w4_sk(3,2) w4_sk(1,3) w4_sk(2,1)]'; w5=[w5_sk(3,2) w5_sk(1,3) w5_sk(2,1)]'; w2local=simplify(A2'*w2); w3local=simplify(A3'*w3); w4local=simplify(A4'*w4); w5local=simplify(A5'*w5); %%POSICION CENTROS GRAVEDAD Rg2=[xc yc 0]'+At2*[0 0 Rt]'; ug3=[x3 0 z3]'; ug4=[x4 0 z4]'; ug5=[x5 0 z5]'; rg3=Rg2+A3*ug3; rg4=Rg2+A3*ug4; rg5=rg4+A4*ug5; %RESTRICCIONES NO HOLONOMAS Rg2=[xc yc 0]'+At2*[0 0 Rt]'; rg2c=At2*[ 0 0 -Rt]'; vg2=[ dxc dyc 0]'+Rt*[cos(phi)*sin(theta)*dphi+sin(phi)*cos(theta)*dtheta; sin(phi)*sin(theta)*dphi-cos(phi)*cos(theta)*dtheta; -sin(theta)*dtheta]; vc=vg2+cross(w2,rg2c); vc=simplify(vc); 39 Análisis de la influencia del ciclista en el comportamiento dinámico de la bicicleta drg5=diff(rg5,xc)*dxc+diff(rg5,yc)*dyc+diff(rg5,phi)*dphi... +diff(rg5,theta)*dtheta+diff(rg5,psi)*dpsi+diff(rg5,beta)*dbeta... +diff(rg5,gamma)*dgamma+diff(rg5,epsilon)*depsilon; vg5=drg5; vd=vg5+cross(w5,A5*Rd*[cos(e2) 0 -sin(e2)]'); vd=simplify(vd); %RESTRICCIONES HOLONOMAS % restriccion de contacto p=[q e2]; dp=[dq de2]; ddp=[ddq dde2]; rp5l=Rd*[cos(e2) 0 -sin(e2)]'; rp5=rg5+A5*rp5l; % coordenada z (rp5(3)) igual a cero es la restriccion rp5=simplify(rp5); tp5l=diff(rp5l,e2); tp5=A5*tp5l; %% coordenada z otra restriccion simplify(rp5); %FICHEROS UTILES qind=[theta psi gamma]; qkin=[xc yc phi epsilon]; qdep=[beta e2]; Ch=[rp5(3);tp5(3)]; Cnh=[vc(1);vc(2);vd(1);vd(2)]; C=[Cnh;Ch]; Cqh=jacobian(Ch,qdep); B=[jacobian(Cnh,dp); jacobian(Ch,p)]; B=simplify(B); dB=jacobian(B*dp',p); matlabFunction(dB,'file', 'dB_bici') matlabFunction(B,'file','B_bici') matlabFunction(C,'file','Cbici') matlabFunction(Cqh,'file','Cqdeph') matlabFunction(Ch,'file','Ch2') %ARCHIVOS PARA LA DINAMICA I2(1,:) = [I2xx 0 0 ]; I2(2,:) = [ 0 I2yy 0 ]; I2(3,:) = [ 0 0 I2zz]; I3(1,:) = [I3xx 0 I3xz]; I3(2,:) = [ 0 I3yy 0 ]; I3(3,:) = [I3xz 0 I3zz]; I4(1,:) = [I4xx 0 I4xz]; I4(2,:) = [ 0 I4yy 0 ]; I4(3,:) = [I4xz 0 I4zz]; I5(1,:) = [I5xx 0 0 ]; I5(2,:) = [ 0 I5yy 0 ]; I5(3,:) = [ 0 0 I5zz]; %Vectores de velocidad V=Hdq a=Hddq+h*dq Código matlab 46 theta,x3,z3); C=Cbici(Rd,Rt,beta,dbeta,depsilon,dgamma,dphi,dpsi,dtheta,dxc,dyc,e2... ,epsilon,gamma,phi,theta,x4,x5,z4,z5); %MONTAR LAS ECUACIONES Mfinal=[ M Cq'; Cq zeros(6)]; dq=[dxc dyc dphi dtheta dpsi dbeta dgamma depsilon de2]; Mdir=0;Mpedal=18;Mvuel=0; alfa1=20; alfa2=20; T=[Mdir Mpedal Mvuel]'; %Est_Baum=-2*alfa1*(dCq*dq')-(alfa2^2)*C; V=[Qg+Qv+Qaero+S*T;-dCq*dq'];%+Est_Baum]; q1=Mfinal\V; ddqind1=q1([4 7 5]); ddtheta=ddqind1(1); ddgamma=ddqind1(2); ddpsi=ddqind1(3); if abs(ddtheta)>=0.05*2*pi ddtheta=0; end ddqind=[ddtheta ddgamma ddpsi]'; dy=[dqind; ddqind; dqkin]; 47 Análisis de la influencia del ciclista en el comportamiento dinámico de la bicicleta 7.5 simulabici.m clear clc global Param Param=Parametros; Rt=Param(3); Rd=Param(4); x3=Param(5); z3=Param(6); x4=Param(7); z4=Param(8); x5=Param(9); z5=Param(10); m2=Param(11); m3=Param(12); m4=Param(13); m5=Param(14); g=Param(29); cv=Param(30); I2xx=Param(15); I2yy=Param(16); I3xx=Param(17); I3yy=Param(18); I3zz=Param(19); I3xz=Param(20); I4xx=Param(21); I4yy=Param(22); I4zz=Param(23); I4xy=Param(24); I4xz=Param(25); I4yz=Param(26); I5xx=Param(27); I5yy=Param(28); Código matlab 48 %%INTEGRAR PARA OBTENER EL VECTOR Y t=[0 5]; q0=[0 0 0 0 0 20 0 0 0 0]'; qdepini=[0 pi/2]'; [t,q]=ode45(@dinamicadirecta,t,q0); % CREAR VECTORES QIND QKIN DQIND qind=q(:,1:3); qkin=q(:,7:10); theta=qind(:,1); gamma=qind(:,2); psi=qind(:,3); xc=qkin(:,1); yc=qkin(:,2); phi=qkin(:,3); epsilon=qkin(:,4); beta=qdepini(1); e2=qdepini(2); %CALCULO DE QDEP qdep=newtonraphsonbici(theta(1),gamma(1),psi(1),xc(1),yc(1),phi(1)... ,epsilon(1),beta,e2); beta(2)=qdep(1); e2(2)=qdep(2); for j=1:length(t)-1 i=j+1; qdep=newtonraphsonbici(theta(i),gamma(i),psi(i),xc(i),yc(i),phi(i)... ,epsilon(i),beta(j),e2(j)); beta(i)=qdep(1); e2(i)=qdep(2); end beta=beta'; e2=e2'; 49 Análisis de la influencia del ciclista en el comportamiento dinámico de la bicicleta qnueva=[xc yc phi theta psi beta gamma epsilon e2]; %%%%%%%%%%%%%%%%%%%%%%% dqind=q(:,4:6); dtheta=dqind(:,1); dgamma=dqind(:,2); dpsi=dqind(:,3); %CALCULO DE DQKIN for i=1:length(beta) dqdep=cinbici(beta(i),epsilon(i),gamma(i),phi(i),psi(i),theta(i),... e2(i),dtheta(i),dgamma(i),dpsi(i)); dxc(i)=dqdep(1); dyc(i)=dqdep(2); dphi(i)=dqdep(3); depsilon(i)=dqdep(4); dbeta(i)=dqdep(5); de2(i)=dqdep(6); end %ANIMACION BICI for i=1:length(t) animaBici(t(i),[xc(i) yc(i) phi(i) theta(i) psi(i) beta(i) gamma(i)... epsilon(i) e2(i)]) end Código matlab 50 %COMPROBACION RESTRICCIONES for i=1:length(t) C(i,:)=Cbici(Rd,Rt,beta(i),dbeta(i),depsilon(i),dgamma(i),dphi(i),... dpsi(i),dtheta(i),dxc(i),dyc(i),e2(i),epsilon(i),gamma(i),phi(i)... ,theta(i),x4,x5,z4,z5); end figure (2) plot(t,C) figure (3) plot (t,qnueva) figure (4) plot (t,q) 8 REFERENCIAS [1] ESCALONA FRANCO, José Luis, Cinemática de la bicicleta [2] ESCALONA FRANCO, José Luis, Dinámica inversa de una bicicleta a partir de datos experimentales [Consulta 11/09/2019] En: https://www.academia.edu/12809356/DIN%C3%81MICA_INVERSA_DE_UNA_BICICLETA_A_PA RTIR_DE_DATOS_EXPERIMENTALES [3] SCHWAB A.L. and MEIJAARD, J.P., Dynamics of Flexible Multibody Systems with Non-Holonomic Constraints: A Finite Element Approach. En: Multibody System Dynamics 10. Holanda: Kluwer Academic Publisher. 2003, 107–123 [4] SHABANA, AHMED A., Computational dynamics New York etc. John Wiley and Sons. 1994 ISBN : 0-471-30551-0 [5] ESCALONA FRANCO, José Luis, Dinámica de sistemas multicuerpo ca.2018 4 Apuntes de la asignatura de Cinemática y dinámica de máquinas. Grado en Ingeniería en Tecnologías Industriales [6] PÉREZ ANAYA, MANUEL, Estudio de la dinámica de la bicicleta con el mínimo número de coordenadas.2018. Trabajo fin de grado ineditado. 51