Herramienta de simulación computacional de maquinaria hidráulica

i Equation Section (Next) Equation Chapter 1 Section 1 Trabajo Fin de Grado Grado en Ingeniería de las Tecnologías Industriales Herramienta de simulación computacional de maquinaria hidráulica Autor: Ignacio Ortega Coronado Tutor: Daniel García Vallejo Dep. Ingeniería Mecánica y Fabricación Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2018 ii iii Trabajo Fin de Grado Grado en Ingeniería de las Tecnologías Industriales Herramienta de simulación computacional de maquinaria hidráulica Autor: Ignacio Ortega Coronado Tutor: Daniel García Vallejo Profesor titular Dep. de Ingeniería Mecánica y Fabricación Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2018 iv v Proyecto Fin de Carrera: Herramienta de simulación computacional de maquinaria hidráulica Autor: Ignacio Ortega Coronado Tutor: Daniel García Vallejo El tribunal nombrado para juzgar el Proyecto arriba indicado, compuesto por los siguientes miembros: Presidente: Vocales: Secretario: Acuerdan otorgarle la calificación de: Sevilla, 2018 El Secretario del Tribunal vi vii A mi familia A mis compañeros viii ix Agradecimientos Han sido cuatro años muy duros de grado, en los que sin el continuo apoyo y afecto de mis compañeros no habría llegado hasta aquí. Por ello agradecerles enormemente su compañía en el día a día, en los malos y los buenos momentos. Agradecer también a mi familia y a mis profesores, con los que he ido aprendiendo en este camino. Ignacio Ortega Coronado Sevilla, 2018 xvi 4 Actuadores hidráulicos 41 4.1 Módulo de compresibilidad B 41 4.2 Teoría de las propiedades concentradas 43 4.2.1 Ecuación de continuidad 44 4.2.2 Relación flujo -presión 45 4.3 Modelo cilindro 46 4.3.1 Volúmenes de control 46 4.3.2 Módulo de compresibilidad 46 4.3.3 Ecuaciones diferenciales 47 4.3.4 Modelado de elementos 47 4.4 Modelo actuador rotatorio 50 4.5 Parámetros 52 4.5.1 General 52 4.5.2 Aceite hidráulico 52 4.5.3 Válvulas 52 4.5.4 Cilindros 53 4.5.5 Actuador rotatorio 54 5 Sistema de control 55 5.1 Fundamento 55 5.2 Referencia 56 5.3 Controlador PD 56 5.4 Sintonización 57 6 Cálculo 59 6.1 Sistema de ecuaciones diferenciales 59 6.2 Condiciones iniciales 60 6.3 Programa MatLab 61 6.3.1 “Simbolico.m” 61 6.3.2 “Restriccionessymb.m” 61 6.3.3 “constantes.m” 62 6.3.4 “Simulacion.m” 62 6.3.5 “funcion.m” 62 6.3.6 Funciones hidráulicas 62 6.3.7 “control.m” 62 6.3.8 Funciones de referencia 63 6.3.9 Funciones de postproceso 63 7 Resultados 65 7.1 Multicuerpo 65 7.1.1 Equilibrio 65 7.1.2 Simulación 66 7.2 Bucle cerrado 69 7.2.1 Condición inicial 69 7.2.2 Parámetros 70 7.2.3 Cilindros 71 7.2.4 Controlador 77 7.2.5 Trayectorias 79 7.2.6 Dinámica 83 7.2.7 Restricciones 85 xvii 8 Conclusiones 89 8.1 Conclusiones 89 8.2 Desarrollo futuro 90 8.2.1 Esquemas de co-simulación 90 8.2.2 Programación del integrador 90 8.2.3 Leyes de control 91 8.2.4 Restricciones no holónomas 91 8.2.5 Reducción de simplificaciones 91 8.2.6 Joystick 91 8.2.7 Realidad virtual 91 Referencias 93 xviii xix Índice de Figuras Figura 1-1 Dolmen en Burren, Irlanda 1 Figura 1-2 Nueva York 2 Figura 1-3 Pala Otis 2 Figura 1-4 Pala Bucyrus 3 Figura 1-5 Excavadora hidráulica de Poclain 3 Figura 1-6 CAT 336 DL 4 Figura 1-7 Cuchara 5 Figura 1-8 Riper 5 Figura 1-9 Martillo 5 Figura 1-10 Pulpo 5 Figura 1-11 Pinza de demolición 5 Figura 1-12 Cizalla 5 Figura 1-13 Excavadora con ruedas 6 Figura 1-14 Excavadora frontal 6 Figura 1-15 Retroexcavadora 7 Figura 2-1 Retroexcavadora 9 Figura 2-2 Sólidos 10 Figura 2-3 Actuadores hidráulicos 11 Figura 2-4 Posición sólido 2 12 Figura 2-5 Ángulos de Cardan del sólido 2 13 Figura 2-6 Coordenadas del sistema 14 Figura 2-7 Mecanismo de 4 barras en el extremo del brazo 17 Figura 2-8 Par de revolución entre los sólidos 6 y 8 17 Figura 2-9 Sólido 2 19 Figura 2-10 Esquema sólido 2 19 Figura 2-11 Sólido 3 20 Figura 2-12 Esquema sólido 3 21 Figura 2-13 Sólido 4 22 Figura 2-14 Esquema sólido 4 22 Figura 2-15 Sólido 5 23 Figura 2-16 Esquema sólido 5 23 xx Figura 2-17 Sólido 6 25 Figura 2-18 Esquema sólido 6 25 Figura 2-19 Sólido 7 26 Figura 2-20 Esquema sólido 7 27 Figura 2-21 Sólido 8 28 Figura 2-22 Esquema sólido 8 28 Figura 3-1 Fuerzas en un par de revolución 37 Figura 4-1 Variación de volumen provocado por una fuerza 41 Figura 4-2 Esquema de un cilindro hidráulico 42 Figura 4-3 Balance de materia 44 Figura 4-4 Caudal en función de la diferencia de presiones 45 Figura 4-5 Esquema del sistema hidráulico 46 Figura 4-6 Esquema válvula estranguladora 47 Figura 4-7 Posiciones de la válvua direccional 49 Figura 4-8 Fuerzas actuantes en el cilindro 49 Figura 4-9 Circuito hidráulico del actuador rotatorio 50 Figura 5-1 Control con realimentación negativa 55 Figura 5-2 Señal de referencia 56 Figura 7-1 Posiciones de equilibrio 66 Figura 7-2 Evolución de las coordenadas del brazo con condiciones iniciales de equilibrio 66 Figura 7-3 Evolución de las coordenadas del brazo sin condiciones iniciales de equilibrio 67 Figura 7-4 Restricciones del brazo 67 Figura 7-5 Restricciones de la base 68 Figura 7-6 Restricciones del brazo con mayor valor de los coeficientes de Baumgarte 68 Figura 7-7 Posición inicial 69 Figura 7-8 Extensión del cilindro de la pluma 71 Figura 7-9 Velocidad de extensión del cilindro de la pluma 71 Figura 7-10 Evolución de las presiones en el cilindro de la pluma 72 Figura 7-11 Extensión del cilindro del balancín 72 Figura 7-12 Velocidad de extensión del cilindro del balancín 73 Figura 7-13 Evolución de las presiones en el cilindro del balancín 73 Figura 7-14 Extensión del cilindro de la cuchara 74 Figura 7-15 Velocidad de extensión del cilindro de la cuchara 74 Figura 7-16 Evolución de las presiones en el cilindro de la cuchara 75 Figura 7-17 Velocidad giro de la cabina 75 Figura 7-18 Giro de la cabina 76 Figura 7-19 Oscilaciones de la velocidad de giro de la cabina 76 Figura 7-20 Evolución de la presión del actuador rotatorio 77 Figura 7-21 Evolución de las señales de control de los cilindros 77 xxi Figura 7-22 Evolución de la señal de control del actuador rotatorio 78 Figura 7-23 Posición final 79 Figura 7-24 Comparación entre la posición inicial y final 79 Figura 7-25 Trayectoria 80 Figura 7-26 Ángulo de giro de la pluma 80 Figura 7-27 Ángulo de giro del balancín 81 Figura 7-28 Ángulo de giro del mecanismo de cuatro barras 81 Figura 7-29 Giro de la cabina 82 Figura 7-30 Fuerza resultante de los cilindros 83 Figura 7-31 Momento resultante del actuador rotatorio 84 Figura 7-32 Fuerzas generalizada aplicadas sobre la coordenada θ4 84 Figura 7-33 Fuerza generalizada resultante sobre la coordenada θ4 85 Figura 7-34 Restricciones del brazo 85 Figura 7-35 Fuerza de reacción en el par de revolución 86 Figura 7-36 Restricciones de la base 86 Figura 7-37 Fuerzas de reacción en la base 87 Figura 7-38 Momentos de reacción en la base 87 xxii xxiii Notación <𝑋,𝑌,𝑍> Sistema de referencia global <𝑋𝑎,𝑌𝑎,𝑍𝑎> Sistema de referencia intermedio girado una vez <𝑋𝑏,𝑌𝑏,𝑍𝑏> Sistema de referencia intermedio girado dos veces <𝑥𝑖,𝑦𝑖,𝑧𝑖> Sistema de referencia local del sólido 𝑖 𝟏𝟑×𝟑 Matriz identidad A Área 𝐚 Vector de aceleraciones 𝑎𝐺𝑖 Aceleración del c.d.g. del sólido 𝑖 𝐀𝜃𝑖 Matriz de giro respecto al eje 𝑦𝑖 𝐀𝜑𝐢 Matriz de giro respecto al eje 𝑥𝑖 𝐀𝜓𝑖 Matriz de giro respecto al eje 𝑧𝑖 𝐚𝐆𝑖 Aceleración lineal del c.d.g. del sólido i expresada en coordenadas globales 𝐀𝑖 Matriz de giro del sólido 𝑖 𝐀𝑗𝑖 Matriz de giro que permite pasar del sistema de referencia local del sólido j al del sólido 𝑖 𝐛 Vector en el sistema de referencia global B Módulo de compresibilidad 𝐛 Vector en el sistema de referencia local 𝐵0 Módulo de compresibilidad del fluido 𝐵𝑐𝑖 Módulo de compresibilidad del contenedor 𝑖 𝐵𝑒 Módulo efectivo de compresibilidad 𝐵𝑒𝑖 Módulo efectivo de compresibilidad del volumen de control 𝑖 𝐂 Vector de restricciones xxiv 𝐶𝑑 Coeficiente de descarga 𝐶ℎ Capacidad hidráulica 𝐂𝐪 Matriz jacobiana de 𝐂 respecto a 𝐪 𝐶(𝑡) Ley de control 𝐂𝑡 Derivada parcial de 𝐂 respecto a 𝑡 𝐂𝑡𝑡 Derivada parcial de 𝐂t respecto a 𝑡 𝐶𝑣 Coeficiente semi-empírico de flujo 𝐶𝑣1 Coeficiente semi-empírico de flujo de la válvula direccional 𝐶𝑣2 Coeficiente semi-empírico de flujo de la válvula estranguladora 𝑑 Diámetro 𝐸 Módulo de elasticidad 𝑒(𝑡) Error del controlador 𝐅ℎ𝑦𝑑 Vector de fuerzas producidas por la acción hidráulica 𝐹𝑗 Fuerza 𝑗 𝐅𝑖𝑗 Vector de fuerza que ejerce i sobre j 𝐹𝑠 Módulo de la fuerza resultante en el cilindro hidráulico 𝐹𝜇 Módulo de la fuerza de fricción en el cilindro hidráulico 𝑮 𝑖 Matriz jacobiana de 𝛚i respecto a 𝐪󰇗 𝐠𝑖 Matriz jacobiana de 𝛚i respecto a 𝐪 g.d.l. Grado(s) de libertad ℎ Extensión máxima del cilindro hidráulico 𝐇𝑖 Matriz jacobiana de 𝐫i respecto a 𝐪 𝐡𝑖 Matriz jacobiana de 𝐯Gi respecto a 𝐪 𝐼𝐺𝑖 Tensor de inercia del sólido i respecto a su c.d.g. 𝐈𝐺𝑖 Tensor de inercia del sólido i respecto a su c.d.g. en coordenadas locales 𝐾𝑑 Ganancia derivativa 𝐾𝑝 Ganancia proporcional 𝐾𝑡 Coeficiente de la válvula de seguridad del actuador rotatorio xxv 𝐋 Matriz que incluye las matrices 𝐇i y 𝐆 i de los sólidos del sistema 𝐿ℎ𝑖 Longitud de la manguera 𝐥 Matriz que incluye las matrices 𝐡i y 𝐠i de los sólidos del sistema 𝐋𝑖 Momento cinético del sólido i expresada en coordenadas locales 𝐿 󰇍  𝑖 Momento cinético del sólido i 𝐌 Matriz de masa generalizada 𝑀 Masa de fluido 𝑚󰇗 Caudal másico 𝐌  Matriz de masa 𝑀 󰇍 󰇍  𝑘 Momento 𝑘 𝑀2 Momento generado por el motor hidráulico 𝐦𝑖 Matriz diagonal cuyos términos son 𝑚𝑖 𝑚𝑖 Masa del sólido 𝑖 𝑛1 Velocidad de rotación de la bomba 𝑛2 Velocidad de giro del motor 𝑝0 Presión de la línea de baja presión del circuito del actuador rotatorio 𝑝ℎ Presión en la línea de alta presión del circuito del actuador rotatorio 𝑝𝑖 Presión del volumen de control 𝑖 𝑝𝑝 Presión de la bomba 𝑝𝑠𝑒𝑡 Presión de seguridad 𝑝𝑇 Presión del tanque 𝐐 Vector de fuerzas generalizadas 𝐪 Vector de coordenadas generalizadas 𝑄 Caudal volumétrico 𝐪󰇗 Vector de velocidades generalizadas 𝐐  Vector de fuerzas externas 𝐪󰇘 Vector de aceleraciones generalizadas 𝐐 𝑟𝑒𝑎𝑐 Vector de fuerzas de reacción Introducción 4 Esto abrió el camino al uso de accionamiento hidráulico en maquinaria de este tipo y el desarrollo que tenemos en la actualidad. En todas las obras de construcción civil es usual encontrar una excavadora, y además el número de fabricantes se ha multiplicado los últimos 50 años. Uno de los fabricantes más reconocidos a nivel mundial por su importancia es Caterpillar. Concretamente en este trabajo se ha emulado la retroexcavadora de gran tamaño CAT 336 DL, Manual Caterpillar [2], (Figura 1-6). Figura 1-6 CAT 336 DL 1.2 Generalidades 1.2.1 Excavadoras La excavadora hidráulica es una maquinaria pesada ampliamente usada en construcción, capaz de desplazarse y girar un determinado ángulo su efector final, que puede ser variado como se verá a continuación. Como su nombre indica, su principal función está en la excavación. No obstante, numerosas excavadoras, como la CAT 336 DL [2] usada como referencia para el trabajo, cuentan con otras funciones según de los accesorios finales empleados: • Cucharones: también denominados cucharas, están destinados a la labor de excavación. En el trabajo se ha tenido en cuenta este elemento final para los cálculos (Figura 1-7). • Riper: herramienta que permite escarificar y descortezar terrenos. Es muy útil para romper y arrastrar roca (Figura 1-8). • Martillos: funcionan por repetición en el sistema hidráulico. Son muy usados en aplicaciones de demolición y construcción (Figura 1-9). • Pulpos: herramientas de agarre que permiten la manipulación de materiales pesados en las construcciones, desplazamiento de vehículos, o transporte de rocas (Figura 1-10). • Pinzas de demolición: permiten la rotación total de las mismas a izquierda y derecha, pudiendo extraer clasificar, manipular y cargar elementos (Figura 1-11). • Cizallas: usadas para el tratamiento de chatarra y corte de estructuras (Figura 1-12). 5 Herramienta de simulación computacional de maquinaria hidráulica Figura 1-7 Cuchara Figura 1-8 Riper Figura 1-9 Martillo Figura 1-10 Pulpo Figura 1-11 Pinza de demolición Figura 1-12 Cizalla En general durante el funcionamiento del brazo de la excavadora, la base móvil de la misma se encuentra parada, para no influir en el desarrollo de la tarea. Si la zona de trabajo se encuentra fuera del alcance del brazo, entonces se procede al desplazamiento hacia una nueva posición, donde ejecuta la operación necesaria una vez estacionada la máquina. Introducción 6 De forma habitual se va a referir de forma indistinta a la máquina como excavadora o retroexcavadora, aunque esta última denominación es más específica puesto que hace referencia a un tipo concreto de excavadora. Para entender la distinción, se pueden identificar las diferentes partes de una máquina de este tipo: • Chasis: permite la traslación de la máquina y es el punto de apoyo del resto de la excavadora. Básicamente se pueden encontrar dos tipos de traslación, mediante orugas o cadenas (Figura 1-6), y mediante ruedas (Figura 1-13). El primero es más habitual para grandes cargas, ya que distribuye mejor la fuerza al terreno al disponer mayor superficie de contacto, pero generalmente son más lentas. La traslación mediante ruedas se caracteriza por tener una baja capacidad de carga, al tener una transmisión de fuerza sobre una superficie mucho menor, pero permite mayor velocidad de desplazamiento. • Estructura: en este trabajo habitualmente se le llamará cabina, aunque también incluye el contrapeso, el motor y la transmisión. Está apoyada sobre la corona de giro de la excavadora, y gira alrededor de un eje vertical sobre el chasis, generalmente es capaz de girar 360𝑜. • Brazo articulado: formado generalmente por la pluma, el brazo o balancín, y el efector final, en este caso la cuchara. En algunos casos, el eslabón intermedio o balancín se omite. Figura 1-13 Excavadora con ruedas Una descripción más precisa de los elementos que componen la excavadora será realizada en el Capítulo 2. Teniendo en cuenta estos elementos, se pueden definir, a grosso modo, dos tipos fundamentales de excavadoras: • Excavadora frontal o pala de empuje: tienen la cuchara hacia arriba. Al cargar hacia arriba tiene mayor altura de descarga, por lo que es útil para cargar el frente de una cantera, por ejemplo (Figura 1-14). • Retroexcavadora: la pala está orientada hacia abajo, por lo que la cuchara avanza moviéndose hacia la cabina, retrocediendo, de ahí su nombre. Son útiles para trabajar por debajo del nivel de la excavadora, por lo que es ideal para la excavación de zanjas o cimentaciones (Figura 1-15). Figura 1-14 Excavadora frontal 7 Herramienta de simulación computacional de maquinaria hidráulica Figura 1-15 Retroexcavadora De igual forma, existen muchas configuraciones y tipos particulares en el mercado, enfocados a tareas específicas, o con unas prestaciones concretas. 1.2.2 Accionamiento hidráulico El sistema de accionamiento proporciona al brazo de la retroexcavadora la energía capaz de moverlo, y determina la velocidad de los movimientos, la fuerza a la que están sometidos los cuerpos y su rendimiento dinámico. Se pueden encontrar en general tres tipos de accionamiento para los sistemas mecánicos: • Eléctrico: se caracteriza por una mayor precisión, ser silenciosos y fácilmente controlables. Tienen la desventaja en este caso de tener una potencia limitada. • Neumático: emplea aire comprimido, y tiene movimientos rápidos, pero de poca precisión. • Hidráulico: emplea aceite para transmitir la potencia. Este último basa su funcionamiento en el principio de Pascal, mediante el cual, la presión ejercida sobre un fluido se transmitirá por todo el dominio fluido y en todas las direcciones. Provocando diferencias de presiones en los cilindros, se genera una fuerza resultante que es la que actúa sobre los diferentes elementos de la máquina. Las principales ventajas de este tipo de accionamiento son: • Relativa rapidez, aunque depende principalmente de la diferencia de presiones generada. • Alta relación potencia peso. De hecho, en este trabajo se desprecia el peso de los cilindros en relación con el peso de los diferentes sólidos sobre los que actúan. • Autolubricación debido al aceite, lo que ayuda a que tengan una mayor durabilidad. • Estabilidad a cargas estáticas. • Alta capacidad de carga. Pueden trabajar a más de 100 𝑏𝑎𝑟, lo que implica una altísima fuerza de accionamiento. Sin embargo, también existen algunas desventajas, entre las que se destaca: • Difícil mantenimiento • Requiere de una instalación compleja. • Las fugas son frecuentes. • En general son equipos caros, ya que requieren de una fabricación precisa con tolerancias relativamente estrechas para reducir la fricción y las fugas. A pesar de estas desventajas, constituyen el accionamiento adecuado para maquinaria pesada considerada. Introducción 8 1.3 Objeto y motivación El objetivo de este trabajo es realizar la simulación de una retroexcavadora. La mayoría de software dedicado a este tipo de asuntos son simuladores específicos que precisan un potente hardware. Se pretende tener una herramienta sencilla de cálculo numérico desarrollada en MatLab, y observar el comportamiento conjunto de un sistema mecánico con otro mucho más rígido como es el hidráulico. Para ello se ha enfocado el problema desde un punto de vista paramétrico, utilizando el lenguaje simbólico. De este modo, para realizar una nueva simulación de otra excavadora, o máquina similar, sólo habría que cambiar los parámetros. Además, permitirá ver la influencia de esos diferentes parámetros, realizando diferentes simulaciones con cambios en ellos. Además de todo esto, se pretende analizar cómo funciona la conjunción de un código destinado a la simulación de un sistema multicuerpo con la simulación de la acción hidráulica, ya que constituyen dos problemas de dinámicas muy diferentes. Por ello se considera interesante comprobar la evolución de las presiones y las coordenadas y la relación que guardan entre sí. 1.4 Estructura del trabajo Para llevar a cabo el trabajo, en primer lugar, hay que definir el modelo 3D, consiguiendo los parámetros geométricos y másicos de los sólidos. También habrá que definir el sistema, las coordenadas que se van a utilizar, y las posibles restricciones a las que esté sometido. Estas cuestiones se desarrollarán el capítulo 2. En el capítulo 3, se llegan a las ecuaciones de la Dinámica de Sistemas Multicuerpo que modelan el movimiento de los sólidos del sistema. Se explica brevemente su origen y el resultado final que se deberá implementar en el programa de cálculo. Por su parte, el capítulo 4 explica lo relacionado con la teoría de las propiedades concentradas para modelar el accionamiento hidráulico de la máquina. Finalmente se plantea el modelo de los diferentes actuadores a implementar y los distintos parámetros hidráulicos que se han seleccionado en función de catálogos y funcionalidad. Una vez que ya se tiene la ecuación a integrar, en el capítulo 5 se plantea la posibilidad de incluir un controlador que manipule la extensión de los cilindros, de manera que se asemeje a la realidad de un operario con los mandos de la excavadora. Tras haber modelado matemáticamente las dos facetas del trabajo, se procede a unir ambos programas en el capítulo 6. El resultado final de este capítulo es la formulación de la ecuación a integrar en MatLab; y se explica brevemente los programas desarrollados en MatLab y la función que lleva a cabo cada uno. En el capítulo 7 se plantean y se analizan los resultados de diferentes simulaciones y los resultados obtenidos, para llegar a las conclusiones del capítulo 8. . 9 2 SISTEMA MULTICUERPO En este capítulo se describe el problema clásico de sistemas multicuerpo al que se aplican las ecuaciones de carácter general para poder obtener el movimiento. También se va a realizar la particularización para el caso concreto de la retroexcavadora, así como las restricciones y simplificaciones realizadas 2.1 Sistema En esta sección se define el sistema y los sólidos que lo componen, así como los actuadores, y las uniones entre los cuerpos (Figura 2-1). Figura 2-1 Retroexcavadora Sistema Multicuerpo 10 2.1.1 Sólidos El sistema está constituido por una retroexcavadora (Figura 2-2). Esta máquina se ha divido en 7 sólidos diferentes. La denominación de sólido 1 se usa para el sistema de referencia global, denominado barra fija. Por tanto, se empieza a numerar a partir del sólido 2, el cual se corresponde con la base móvil o chasis de la excavadora, formada por el tren de rodaje para el desplazamiento de la máquina. El sólido 3 se corresponde con el grueso de la excavadora, que incluye la cabina del operario, que gira alrededor del eje vertical, este conjunto ha recibido el nombre de estructura en el apartado 1.2.1. Los sólidos 4 y 5 hacen referencia al brazo y antebrazo, respectivamente, a los que se les llamará pluma y balancín. El sólido 6 es la cuchara de la máquina. Finalmente, los sólidos 7 y 8, en cuya unión actúa un cilindro hidráulico, constituyen junto con el sólido 6, un mecanismo de 4 barras cuya función es hacer más efectivo el momento generado en el sólido 6. Podremos tener una idea mucho más exacta en la sección 2.5. Figura 2-2 Sólidos 11 Herramienta de simulación computacional de maquinaria hidráulica 2.1.2 Actuadores Teniendo en cuenta la configuración típica de una máquina de este tipo, habrá una actuación gracias a un motor de combustión en el sólido 2, para efectuar el desplazamiento de la máquina en sí y dar energía al sistema hidráulico. Puesto que no es el objeto del trabajo, se va a considerar inicialmente que la máquina está parada, evitando entrar en el modelado de dicho motor de combustión y poder centrarse en la combinación de los actuadores hidráulicos y la Dinámica de Sistemas Multicuerpo. El giro de la cabina se debe a la actuación de un motor hidráulico que ejerce un par de rotación, mientras que el giro de los diferentes sólidos del brazo articulado se realizará mediante cilindros hidráulicos, cuya actuación se modela en el capítulo 4. Destacar que de cara al análisis del sistema multicuerpo no se va a tener en cuenta estos actuadores como sólidos adicionales, puesto que complicaría enormemente el problema sin aportar ningún cambio conceptual en el análisis de la máquina. Su efecto se tendrá en cuenta como una fuerza aplicada en puntos concretos de los sólidos. Figura 2-3 Actuadores hidráulicos En la Figura 2-3 aparecen resaltados los actuadores hidráulicos a considerar. El motor hidráulico no aparece en la figura, pero funcionará como un par sobre el cuerpo principal de la máquina. 2.1.3 Uniones A excepción del sólido 2 con la barra fija, todas las uniones del sistema las constituyen pares de revolución que solo permiten el giro relativo entre los diferentes sólidos. El sólido 2 y su comportamiento se explica en las secciones 2.2.1 y 2.3.1. Sistema Multicuerpo 12 2.2 Coordenadas El objetivo fundamental en esta sección es describir las coordenadas que definen la geometría del problema que se va a resolver. Un aspecto fundamental de estas coordenadas es que éstas no pueden tomar cualquier valor, sino que están sometidas a unas ecuaciones de restricción, que se definen en la sección 2.3. De este modo, los posibles valores están relacionados. 2.2.1 Selección de coordenadas Para la descripción de los sistemas y las diferentes relaciones hay diferentes tipos de coordenadas que se pueden usar con sus ventajas e inconvenientes. Debido a la configuración de cadena abierta del brazo de la retroexcavadora, y que a partir de la situación de la base de esta solo hay pares de revolución, se ha optado por elegir coordenadas relativas que añaden un giro por cada cuerpo adicional. De este modo el resultado de resolver las ecuaciones es directamente la rotación de la junta. Además, los actuadores hidráulicos actúan entre un cuerpo y el siguiente de manera relativa. No obstante, habría otras opciones, como usar coordenadas lagrangianas, de referencia o naturales; pero como se ha dicho, las coordenadas relativas tienen una relación directa con el problema. Realmente, se trata de un problema 3D, a pesar de que el brazo se mueva en un solo plano por su configuración geométrica. Para el posicionamiento del brazo se debe posicionar la base de la retroexcavadora, y para ello tendremos que usar coordenadas 3D. Para la definición de la base hay que elegir 6 coordenadas que tengan en cuenta los 6 grados de libertad que tiene un sólido en el espacio. Es directo suponer que la posición del centro de gravedad (c.d.g) de la base, a la que llamaremos sólido 2, se define con las coordenadas cartesianas 𝑥2,𝑦2 y 𝑧2 (Figura 2-4). Por su parte, la orientación de un sólido en el espacio es algo más compleja. 𝐫2=[𝑥2 𝑦2 𝑧2] (2–1) Figura 2-4 Posición sólido 2 13 Herramienta de simulación computacional de maquinaria hidráulica Es importante definir la matriz de rotación de un sistema de referencia. Teniendo un sistema de referencia global <𝑋,𝑌,𝑍>, y un sistema de referencia local <𝑥,𝑦,𝑧>, para obtener las componentes en el sistema de referencia global de un vector expresado en el sistema de referencia local es necesaria la matriz de rotación. 𝐛=𝐀 𝐛 (2–2) De este modo la matriz 𝐀 da cuenta del giro necesario entre los sistemas de referencia, y se usará para definir las posiciones de los diferentes sólidos. De este modo, los ángulos que se usan para definir la matriz de giro son los mismos que describen la orientación del sólido en el espacio. De nuevo, son varias las opciones que permiten expresar la orientación de un sólido en el espacio. En este caso se va a usar una variante de los ángulos de Euler. Los ángulos de Euler se basan en realizar 3 giros al sistema de referencia global hasta hacerlo coincidir con el local. Estos giros son sucesivos y dependiendo del orden en que se hagan y respecto a qué ejes, surgen diferentes opciones. Una de ellas es la que se considera a continuación, en los cuales primero se gira respecto al eje 𝑍, luego respecto al eje 𝑋𝑎, girado ya una vez, y por último respecto al eje 𝑌𝑏, ya girado dos veces. Toda combinación de tres coordenadas utilizadas para determinar la orientación de un sistema de referencia puede tener una fácil comprensión física, no obstante, tienen la particularidad de que en ciertas configuraciones geométricas existe una singularidad que conlleva errores numéricos. La variante de los ángulos de Euler utilizada no es una excepción, si bien para alcanzar esa singularidad en un vehículo (que es el típico caso de uso de estos ángulos), es difícil llegar a ella a menos que estemos en condiciones extremas. Una vez explicado todo eso llamaremos 𝜓2 al giro respecto a 𝑍, 𝜑2 al giro respecto a 𝑋𝑎, y 𝜃2 al giro respecto al 𝑌𝑏 (Figura 2-5). En lo que sigue se usará 𝜑, 𝜃 y 𝜓 para denotar giros respecto a los ejes 𝑥, 𝑦 y 𝑧 respectivamente. Figura 2-5 Ángulos de Euler del sólido 2 Sistema Multicuerpo 20 Los parámetros másicos del sólido 2 son: 𝑚2=13019.89 𝑘𝑔 (2–33) 𝐈𝐺2=[18921.239 0 0 0 21975.593 0 0 0 38119.584] 𝑘𝑔 𝑚2 (2–34) La posición del punto A viene definida por: 𝐫2𝐴=𝐫𝟐+𝐀𝟐𝐮 2𝐴 (2–35) En esta ecuación: 𝐫𝟐=[𝑥2 𝑦2 𝑧2] 𝑚 𝐮 2𝐴=[ 0 0 0.630] 𝑚 (2–36) 2.5.1.2 Sólido 3 El sólido 3, donde se ubica la cabina del operario tiene una geometría no simétrica. El bloque situado sobre el cuerpo es un contrapeso que tiene como objetivo garantizar que el c.d.g de la máquina no esté demasiado adelantado en condiciones de trabajo, impidiendo el vuelco la excavadora. Los puntos de interés en este caso son el mismo punto A anterior, el punto B donde se apoya el sólido 4, y el punto PCP3 donde se apoya el cilindro que impulsa la pluma (Figura 2-12). Se va a hacer una simplificación puesto que realmente el c.d.g. del sólido 3 se encuentra desviado en 𝑦3 4 𝑐𝑚 de la mitad de la anchura por la asimetría. Como esta medida es despreciable respecto a las medidas de la máquina se supone que el c.d.g. se encuentra en el centro geométrico. Figura 2-11 Sólido 3 21 Herramienta de simulación computacional de maquinaria hidráulica 4 Figura 2-12 Esquema sólido 3 Los parámetros másicos del sólido 3 son: 𝑚3=17506.977 𝑘𝑔 (2–37) 𝐈𝐺3=[15918.865 −1068.892 −1490.856 −1068.892 30432.515 476.218 −1490.856 476.218 40501.375] 𝑘𝑔 𝑚2 (2–38) La posición del c.d.g. del sólido 3 se obtiene de: 𝐫𝟑=𝐫2𝐴−𝐀𝟑𝐮 3𝐴 (2–39) Siendo 𝐮 3𝐴 según el cuerpo modelado en la Figura 2-12: 𝐮 3𝐴=[1.538 0 −0.638] 𝑚 (2–40) La posición del punto B y del punto PCP3 vienen definidas por: 𝐫3𝐵=𝐫𝟑+𝐀𝟑𝐮 3 𝐵 (2–41) 𝐮 3 𝐵=[1.394 0 0.182] 𝑚 (2–42) 𝐫3𝑃𝐶𝑃3=𝐫𝟑+𝐀𝟑𝐮 3 𝑃𝐶𝑃3 (2–43) 𝐮 3 𝑃𝐶𝑃3=[2.834 0 −0.240] 𝑚 (2–44) Sistema Multicuerpo 22 2.5.1.3 Sólido 4 El sólido 4 (Figura 2-13) se corresponde con la pluma de la excavadora. Es el primer elemento del brazo hidráulico de la máquina. Se tienen tres puntos importantes. En primer lugar, el punto B de unión al sólido 3. El punto C sirve de unión entre el sólido 4 y 5. También están los puntos de aplicación de las acciones hidráulicas, que son PCP4 y PCB4 (Figura 2-14). Figura 2-13 Sólido 4 Figura 2-14 Esquema sólido 4 Los parámetros másicos del sólido 4 son: 𝑚4=3029.87 𝑘𝑔 (2–45) 𝐈𝐺4=[494.064 0 −125.700 0 10285.189 0 −125.700 0 10008.860] 𝑘𝑔 𝑚2 (2–46) La posición del c.d.g. del sólido 4 se obtiene de: 𝐫𝟒=𝐫3𝐵−𝐀𝟒𝐮 4 𝐵 (2–47) Siendo 𝐮 4 𝐵 según el cuerpo modelado en la Figura 2-14: 𝐮 4 𝐵=[−3.036 0 −0.547] 𝑚 (2–48) 23 Herramienta de simulación computacional de maquinaria hidráulica La posición del punto C, del punto PCP4 y PCB4 vienen definidas por: 𝐫4𝐶=𝐫𝟒+𝐀𝟒𝐮 4 𝐶 (2–49) 𝐮 4 𝐶=[3.341 0 −0.425] 𝑚 (2–50) 𝐫4𝑃𝐶𝑃4=𝐫𝟒+𝐀𝟒𝐮 4 𝑃𝐶𝑃4 (2–51) 𝐮 4 𝑃𝐶𝑃4=[−0.677 0 0.357] 𝑚 (2–52) 𝐫4𝑃𝐶𝐵4=𝐫𝟒+𝐀𝟒𝐮 4 𝑃𝐶𝐵4 (2–53) 𝐮 4 𝑃𝐶𝐵4=[0.413 0 0.690] 𝑚 (2–54) 2.5.1.4 Sólido 5 El sólido 5 (Figura 2-15) es el balancín, que ejerce de unión entre la pluma y la cuchara, aumentando el alcance de la excavadora. En él se encuentran el punto C, que lo une al sólido 4, el punto D, que lo une a la cuchara, y el punto E, que lo une al sólido 7. A esto hay que añadir los puntos de actuación de los cilindros hidráulicos PCB5 y PCC7 (Figura 2-16). Figura 2-15 Sólido 5 Figura 2-16 Esquema sólido 5 Sistema Multicuerpo 24 Los parámetros másicos del sólido 5 son: 𝑚5=1274.076 𝑘𝑔 (2–55) 𝐈𝐺5=[161.665 0 −219.003 0 1766.645 0 −219.003 01687.027] 𝑘𝑔 𝑚2 (2–56) La posición del c.d.g. del sólido 5 se obtiene de: 𝐫𝟓=𝐫4𝐶−𝐀𝟓𝐮 5 𝐶 (2–57) Siendo 𝐮 5 𝐶 según el cuerpo modelado en la Figura 2-16: 𝐮 5 𝐶=[−0.801 0 −0.306] 𝑚 (2–58) La posición del punto D, E, PCB5 y PCC5 vienen definidas por: 𝐫5𝐷=𝐫𝟓+𝐀𝟓𝐮 5 𝐷 (2–59) 𝐮 5 𝐷=[2.243 0 −0.297] 𝑚 (2–60) 𝐫5𝐸=𝐫𝟓+𝐀𝟓𝐮 5 𝐸 (2–61) 𝐮 5 𝐸=[1.804 0 −0.242] 𝑚 (2–62) 𝐫5𝑃𝐶𝐵5=𝐫𝟓+𝐀𝟓𝐮 5 𝑃𝐶𝐵5 (2–63) 𝐮 5 𝑃𝐶𝐵5=[−1.969 0 0.109] 𝑚 (2–64) 𝐫5𝑃𝐶𝐶5=𝐫𝟓+𝐀𝟓𝐮 5 𝑃𝐶𝐶5 (2–65) 𝐮 5 𝑃𝐶𝐶5=[−0.407 0 0.580] 𝑚 (2–66) 25 Herramienta de simulación computacional de maquinaria hidráulica 2.5.1.5 Sólido 6 El sólido 6 (Figura 2-17) se corresponde con la cuchara, el elemento terminal de la excavadora. Tiene como puntos a determinar D de unión al sólido 5, y F de unión al sólido 8 (Figura 2-18). Figura 2-17 Sólido 6 Figura 2-18 Esquema sólido 6 Sistema Multicuerpo 26 Los parámetros másicos del sólido 6 son: 𝑚6=1139.118 𝑘𝑔 (2–67) 𝐈𝐺6=[388.548 0 −85.207 0387.406 0 −85.207 0304.194] 𝑘𝑔 𝑚2 (2–68) La posición del c.d.g. del sólido 6 se obtiene de: 𝐫𝟔=𝐫5𝐷−𝐀𝟔𝐮 6 𝐷 (2–69) Siendo 𝐮 6 𝐷 según el cuerpo modelado en la Figura 2-18: 𝐮 6 𝐷=[−0.837 0 0.196] 𝑚 (2–70) La posición del punto F viene determinada por: 𝐫6𝐹=𝐫𝟔+𝐀𝟔𝐮 6 𝐹 (2–71) 𝐮 6 𝐹=[−0.512 0 0.611] 𝑚 (2–72) 2.5.1.6 Sólido 7 El sólido 7 (Figura 2-19) es una barra en cuyo extremo (común al punto inicial del sólido 8) actúa el cilindro hidráulico de la cuchara. Los puntos significativos de este sólido son E, de unión al sólido 5, y PCC8 de acción del actuador (Figura 2-20). Figura 2-19 Sólido 7 27 Herramienta de simulación computacional de maquinaria hidráulica Figura 2-20 Esquema sólido 7 . Los parámetros másicos del sólido 7 son: 𝑚7=96.438 𝑘𝑔 (2–73) 𝐈𝐺7=[4.571 0 0 0 6.318 0 0 0 10.422] 𝑘𝑔 𝑚2 (2–74) La posición del c.d.g. del sólido 7 se obtiene de: 𝐫𝟕=𝐫5𝐸−𝐀𝟕𝐮 7 𝐸 (2–75) Siendo 𝐮 7 𝐸 según el cuerpo modelado en la Figura 2-20: 𝐮 7 𝐸=[−0.487 0 0] 𝑚 (2–76) La posición del punto PCC8 viene determinada por: 𝐫7𝑃𝐶𝐶8=𝐫𝟕+𝐀𝟕𝐮 7 𝑃𝐶𝐶8 (2–77) 𝐮 7 𝑃𝐶𝐶8=[0.324 0 0] 𝑚 (2–78) Sistema Multicuerpo 28 2.5.1.7 Sólido 8 Por último, el sólido 8 (Figura 2-21) conecta a 7 y 6, y a través de este la cuchara sufre la acción del cilindro hidráulico. Los puntos significativos de este sólido son PCC8 y F (Figura 2-22), que es el punto de revolución donde se impone la restricción en las ecuaciones (2–29) y (2–30). Figura 2-21 Sólido 8 Figura 2-22 Esquema sólido 8 . Los parámetros másicos del sólido 8 son: 𝑚8=255.775 𝑘𝑔 (2–79) 𝐈𝐺8=[6.882 0 0 022.575 0 0 0 26.464] 𝑘𝑔 𝑚2 (2–80) La posición del c.d.g. del sólido 8 se obtiene de: 𝐫𝟖=𝐫7𝐸−𝐀𝟖𝐮 8 𝑃𝐶𝐶8 (2–81) 29 Herramienta de simulación computacional de maquinaria hidráulica Siendo 𝐮 8 𝑃𝐶𝐶8 según el cuerpo modelado en la Figura 2-22: 𝐮 8 𝑃𝐶𝐶8=[−0.428 0 0] 𝑚 (2–82) La posición del punto F viene determinada por: 𝐫8𝐹=𝐫𝟖+𝐀𝟖𝐮 8 𝐹 (2–83) 𝐮 8 𝐹=[0.285 0 0] 𝑚 (2–84) Por tanto, se aprecia aquí que hay dos caminos para llegar al punto F y que la restricción comentada anteriormente surge de garantizar que lo sólidos 6 y 8 están conectados en ese punto. Ecuación del movimiento 36 Teniendo (3–15) y (3–26) en cuenta, se puede expresar 𝐚 en función de las coordenadas generalizadas: 𝐚= [ 𝐚𝐺2 ⋮ 𝐚𝐺𝑛𝑐 𝛂2 ⋮ 𝛂𝑛𝑐 ] = [ 𝐇2 ⋮ 𝐇𝑛𝑐 𝐆 2 ⋮ 𝐆 𝑛𝑐 ] 𝐪󰇘+ [ 𝐡2 ⋮ 𝐡𝑛𝑐 𝐠2 ⋮ 𝐠𝑛𝑐 ] 𝐪󰇗 (3–27) Denominando de forma diferente esas matrices: 𝐋= [ 𝐇2 ⋮ 𝐇𝑛𝑐 𝐆 2 ⋮ 𝐆 𝑛𝑐 ] 𝐥= [ 𝐡2 ⋮ 𝐡𝑛𝑐 𝐠2 ⋮ 𝐠𝑛𝑐 ] (3–28) La expresión (3–27) pasa a ser: 𝐚=𝐋𝐪󰇘+𝐥𝐪󰇗 (3–29) Sustituyendo la expresión (3–29) en las ecuaciones de Newton-Euler (3–13) para sistemas multicuerpo, se llega a: 𝐌 (𝐋𝐪󰇘+𝐥𝐪󰇗)=𝐐 +𝐐 𝑣 (3–30) Si se sigue manipulando esta expresión se obtiene la expresión general del equilibrio de un sistema multicuerpo en coordenadas generalizadas: 𝐌𝐪󰇘=𝐐+𝐐𝑣 (3–31) Donde: 𝐌=𝐋𝑇𝐌 𝐋 𝐐=𝐋𝑇𝐐 =𝐐𝑎𝑝𝑙+𝐐𝑟𝑒𝑎𝑐 𝐐𝑣=𝐋𝑇(𝐐𝑣−𝐌 𝐥𝐪󰇗) (3–32) 37 Herramienta de simulación computacional de maquinaria hidráulica 3.1.3.1 Términos de las ecuaciones de Newton-Euler en coordenadas generalizadas En la ecuación (3–31) aparecen una serie de términos que se explican a continuación: • −𝐌𝐪󰇘: es el término referente a las fuerzas de inercia generalizadas, proporcionales a las aceleraciones generalizadas. La matriz 𝐌 de masa generalizada no es constante, puesto que depende de 𝐪, es decir cambia para cada posición del sistema aún a pesar de que la matriz 𝐌  sí que era constante en el tiempo. Esta variabilidad está propiciada por la matriz 𝐋, que incluye los jacobianos de la posición y orientación de los sólidos que sí dependen de 𝐪. • 𝐐𝑣: son las fuerzas cuadráticas en velocidad, que incluyen 𝐐 𝑣∝−𝛚𝑖×(𝐈𝐺2𝛚𝑖) y 𝐌 𝐥𝐪󰇗. Estos son los términos asociados a 𝐪󰇗. Las aceleraciones de los sólidos dependen de las velocidades generalizadas debido a las fuerzas centrífugas, que por definición dependen de la velocidad al cuadrado, de ahí el nombre de fuerzas cuadráticas en velocidad. Estos términos son computacionalmente más difíciles de procesar. • 𝐐: este término se corresponde con las fuerzas externas generalizadas. Incluye tanto a las fuerzas aplicadas (𝐐𝑎𝑝𝑙), como a las fuerzas de reacción (𝐐𝑟𝑒𝑎𝑐). Las fuerzas aplicadas se pueden deber a diversas causas, como muelles, la gravedad, aerodinámicas,… En este caso las fuerzas aplicadas son la gravedad, y las fuerzas debidas al sistema hidráulico. Las fuerzas aplicadas dependen de las coordenadas y sus derivadas, por lo que se pueden conocer de manera más o menos sencilla. Por su parte las fuerzas de reacción están asociadas a las restricciones impuestas. Son las fuerzas que sufren los sólidos por el hecho de tener que cumplir las restricciones durante su movimiento. Estas fuerzas no se pueden conocer a priori y son difíciles de calcular. A continuación, se usará el Teorema de los multiplicadores de Lagrange para tratar este problema. 3.1.4 Teorema de los multiplicadores de Lagrange En numerosa bibliografía, como A.A. Shabana [6] o [3], se trata este teorema enfocado a sistemas multicuerpo. Sean dos sólidos (𝑖 y 𝑗) entre los que existe una restricción, el sólido 𝑗 ejercerá una fuerza sobre el sólido 𝑖, a la que llamaremos 𝐅𝑗𝑖, y por acción y reacción, la fuerza 𝐅𝑖𝑗 que sufrirá 𝑗 por la acción de 𝑖, cumplirá: 𝐅𝑖𝑗=−𝐅𝑗𝑖 (3–33) Figura 3-1 Fuerzas en un par de revolución Ecuación del movimiento 38 Llamando: 𝐅𝑗𝑖=𝛌 𝐅𝑖𝑗=−𝛌 (3–34) Si la restricción a la que están asociadas esas fuerzas es 𝐂, se puede demostrar que: 𝐐𝑟𝑒𝑎𝑐=𝐋𝑇𝐐 𝑟𝑒𝑎𝑐=𝐂𝐪𝛌 (3–35) De manera que se obtienen las fuerzas de reacción generalizadas. Sustituyendo (3–35) en (3–31): 𝐌𝐪󰇘=𝐐𝑎𝑝𝑙+𝐐𝑣+𝐂𝐪𝛌 (3–36) 3.2 Sistema de ecuaciones diferenciales-algebraicas (DAE) Una vez definidas las restricciones geométricas que debe cumplir el sistema (2–31), y las ecuaciones de equilibrio desarrolladas en coordenadas generalizadas (3–36), en esta sección se van a unir ambas condiciones que formulan el problema que resolver, dando lugar a un sistema de ecuaciones diferenciales-algebraicas. 3.2.1 Ecuaciones del movimiento La ecuación (2–31), implica que las derivadas en el tiempo de 𝐂(𝐪) deben ser igual a cero. Por ello desarrollando la expresión teniendo en cuenta que en este caso no se va a depender explícitamente del tiempo: 𝐂󰇗(𝐪,𝐪󰇗)=𝑑𝐂(𝐪) 𝑑𝑡 =𝜕𝐂(𝐪) 𝜕𝐪 𝐪󰇗+𝜕𝐂(𝐪) 𝜕𝑡 =𝐂𝐪𝐪󰇗+𝐂𝑡=𝟎 𝐂𝑡=𝟎 (3–37) Derivando de nuevo respecto al tiempo: 𝐂󰇘(𝐪,𝐪󰇗,𝒒󰇘)=𝑑𝐂󰇗(𝐪) 𝑑𝑡 =𝐂𝐪𝐪󰇘+𝐂󰇗𝐪𝐪󰇗+𝐂󰇗𝑡=𝟎 𝐂󰇗𝑡=𝐂𝑡𝑡=𝟎 (3–38) En (3–38), el término 𝐂󰇗𝐪𝐪󰇗 se puede obtener como se explica en la bibliografía [4]: 𝐂󰇗𝐪𝐪󰇗=𝜕(𝐂𝐪𝐪󰇗) 𝜕𝐪 𝐪󰇗=(𝐂𝐪𝐪󰇗)𝐪𝐪󰇗 (3–39) La ecuación (3–38) queda, por tanto: 𝐂𝐪𝐪󰇘=−(𝐂𝐪𝐪󰇗)𝐪𝐪󰇗−𝐂𝑡𝑡 𝐂𝐪𝐪󰇘=−(𝐂𝐪𝐪󰇗)𝐪𝐪󰇗 (3–40) 39 Herramienta de simulación computacional de maquinaria hidráulica 3.2.2 Ecuaciones de equilibrio En la ecuación (3–36) se puede mover al otro miembro y cambiar de signo 𝐂𝐪𝛌, de manera que ese término representaría el opuesto de las fuerzas de reacción generalizadas, y quedaría el equilibrio de fuerzas: 𝐌𝐪󰇘+𝐂𝐪𝛌=𝐐𝑎𝑝𝑙+𝐐𝑣 (3–41) 3.2.3 Sistema de ecuaciones diferenciales-algebraicas (DAE) Uniendo (3–40) y (3–41) se llega a: [𝐌 𝐂𝐪 𝑇 𝐂𝐪𝟎][𝐪󰇘𝛌]=[𝐐𝑎𝑝𝑙+𝐐𝑣 −(𝐂𝐪𝐪󰇗)𝐪𝐪󰇗]=[𝐐 𝐐𝑑] (3–42) Este sistema es lo que se llama un sistema de ecuaciones diferenciales algebraicas. La parte de ecuaciones diferenciales se debe a que el término 𝐪󰇘 está en función de términos que dependen de 𝐪 y de 𝐪󰇗. Por su parte el nombre de algebraica se debe a la aparición de los términos 𝛌, los cuales representan una simple ecuación algebraica en la que, una vez resuelta la ecuación diferencial, se pueden obtener los multiplicadores de Lagrange. En este caso se tiene un DAE de orden uno, lo que significa que operando con las matrices se puede desacoplar el sistema de manera que se obtenga de forma independiente: 𝐪󰇘=𝒇(𝐪,𝐪󰇗,𝑡) (3–43) Esta ecuación es una ecuación diferencial no lineal que se resolverá más adelante usando los métodos matemáticos pertinentes. También se pueden desarrollar las matrices de modo que se obtengan los multiplicadores de Lagrange: 𝛌=𝐠(𝐪,𝐪󰇗,𝑡) (3–44) Estas ecuaciones son, por tanto, las que modelan el funcionamiento de la excavadora. Como se puede observar, el número de restricciones es menor que el número de coordenadas, teniendo 4 g.d.l. En algunos mecanismos se incluyen las ya mencionadas restricciones de movilidad, de modo que se determina, a priori, la evolución de las coordenadas. Cuando se incluyen tantas restricciones de movilidad como g.d.l. tiene el sistema, se puede obtener las coordenadas generalizadas en cada instante, e ir obteniendo las velocidades, aceleraciones y fuerzas con las ecuaciones definidas en los apartados 3.2.1 y 3.2.2, de manera que el sistema sólo se resolvería mediante ecuaciones algebraicas. Este tipo de simulación se conoce como dinámica inversa. La realidad no funciona así, la verdadera simulación es la dinámica directa, que integra directamente las DAE descritas anteriormente. El número de coordenadas es mayor que el número de restricciones. El movimiento es resultado de las fuerzas aplicadas junto a las restricciones impuestas, por lo que es necesario integrar las DAE. Éste es el tipo de análisis que se está llevando a cabo en este trabajo. Realizando el desarrollo comentado anteriormente, manipulando la ecuación (3–42): 𝐪󰇘=𝐌−1(𝐐+𝐂𝐪 𝑇(𝐂𝐪𝐌−1𝐂𝐪 𝑇)−1(𝐐𝑑−𝐂𝐪𝐌−1𝐐)) (3–45) 𝛌=−(𝐂𝐪𝐌−1𝐂𝐪 𝑇)−1(𝐐𝑑−𝐂𝐪𝐌−1𝐐) (3–46) La ecuación (3–45) es una ecuación diferencial ordinaria no lineal, implementable en un programa de cálculo de forma sencilla. Los multiplicadores de Lagrange de la ecuación (3–46) se obtendrían de forma algebraica una vez resueltas posición, velocidad y aceleración. Ecuación del movimiento 40 3.2.4 Estabilización de Baumgarte En el sistema de ecuaciones (3–42) se han derivado dos veces las ecuaciones de restricción para poder compactar todas las ecuaciones del sistema. Según J. Baumgarte [7], las restricciones del sistema son de tipo 1, holónomas. El hecho de que 𝐂󰇘(𝐪,𝐪󰇗,𝒒󰇘)=𝟎, es una consecuencia de la verdadera relación 𝐂(𝐪)=𝟎. La forma de que se cumpla la verdadera restricción sería dar unas condiciones iniciales que cumplan las restricciones, de modo que al integrar se seguirían cumpliendo. Realmente, a partir de cierto número de pasos de integración, es posible que la computadora cometa un pequeño error 𝐂(𝐪)=𝛅, y 𝐂󰇗(𝐪,𝐪󰇗)=𝛜, de modo que los siguientes pasos arrastrarían ese error siguiendo la ley: 𝐂(𝐪,𝑡)=𝝐𝑡+𝜹 (3–47) Pero al derivar dos veces esa expresión, 𝐂󰇘(𝐪,𝐪󰇗,𝐪󰇘) sigue siendo cero, por lo que para la computadora no hay ningún error, cuando realmente no se están cumpliendo las restricciones. Para evitar la inestabilidad de las restricciones, se fuerza a que se cumplan de alguna forma, por lo que Baumgarte propuso el uso de coeficientes de estabilización, 𝛼 y 𝛽. 𝐂󰇘(𝐪,𝐪󰇗,𝐪󰇘)+2𝛼𝐂󰇗(𝐪,𝐪󰇗)+𝛽2𝐂(𝐪)=𝟎, 𝛼>0 (3–48) Como en las condiciones iniciales del sistema, se va a cumplir 𝐂(𝐪𝟎)=𝟎, y 𝐂󰇗(𝐪𝟎,𝐪󰇗𝟎)=𝟎, se puede concluir que (3–48) es analíticamente equivalente a (3–38). Desde el punto de vista computacional, el uso de 2𝛼𝐂󰇗(𝐪,𝐪󰇗)+𝛽2𝐂(𝐪), garantiza la estabilidad del sistema con el debido cumplimiento de las restricciones. Normalmente se elige 𝛽=𝛼, y cuanto más alto se elijan esos coeficientes, con mayor exactitud se cumplirán las restricciones, a costa de un mayor tiempo de compilación. Desarrollando (3–48), y utilizando las ecuaciones (3–37), (3–38) y (3–39) se llega a: 𝐂𝐪𝐪󰇘=−(𝐂𝐪𝐪󰇗)𝐪𝐪󰇗−2𝛼𝐂𝐪𝐪󰇗−𝛽2𝐂=𝐐𝑑 (3–49) Por lo que se ha modificado el término 𝐐𝑑, para garantizar la estabilidad del sistema y el cumplimiento de las restricciones. 41 4 ACTUADORES HIDRÁULICOS En este capítulo se van a desarrollar las ecuaciones que permitirán obtener las fuerzas resultantes en los cilindros, para lo que se usa la teoría de las propiedades concentradas. Finalmente, se realizará el cálculo para un circuito hidráulico que es el que va a alimentar los diferentes cilindros, describiendo las propiedades concretas de los diferentes componentes del sistema. Para el modelo de la hidráulica, se ha tomado como base A. Mikkola [8]. 4.1 Módulo de compresibilidad 𝑩 Las máquinas hidráulicas emplean aceite hidráulico que es compresible para grandes esfuerzos. Esta compresibilidad hace que el fluido funcione como un muelle, afectando a la dinámica del sistema. Figura 4-1 Variación de volumen provocado por una fuerza La variación de volumen provocada por la fuerza implicará un cambio en la presión del fluido, y la relación entre la presión y la variación de volumen es el llamado módulo de compresibilidad. 𝐵=− 𝑑𝑝 𝑑𝑉 𝑉 ⁄=−𝑉 𝑑𝑝 𝑑𝑉 (4–1) El signo negativo se debe a que una reducción de volumen implica un aumento de presión, y un aumento de volumen una reducción de la presión en el fluido. Aunque el módulo de compresibilidad de un fluido es dependiente de la presión y la temperatura, en la práctica se puede suponer constante. Un valor típico es 𝐵=1860 𝑀𝑃𝑎. Actuadores hidráulicos 42 A esta rigidez, además, hay que añadir la rigidez del contenedor del fluido, y la posible compresibilidad de aire disuelto, que afecta severamente al valor de 𝐵. La unión de todos los factores da lugar al módulo de compresibilidad efectivo 𝐵𝑒. Típicamente para considerar el volumen se tienen en cuenta las mangueras y el contenedor de fluido que en este caso es el cilindro hidráulico (Figura 4-2). Figura 4-2 Esquema de un cilindro hidráulico En la Figura 4-2, el volumen total de fluido sería el del volumen del cilindro más el de la manguera: 𝑉𝑇=𝑉1+𝑉2=𝑉𝑐1+𝑉𝑐2 (4–2) Donde 𝑉𝑐1 y 𝑉𝑐2 son los volúmenes de los contenedores. La variación de volumen de fluido vendrá dada por la compresión de fluido y por el aumento de volumen de los contenedores debido al aumento de presión. 𝑑𝑉𝑇=𝑑𝑉1+𝑑𝑉2−𝑑𝑉𝑐1−𝑑𝑉𝑐2 (4–3) Usando la ecuación (4–1), se puede expresar: −𝑑𝑉𝑡=𝑉𝑡𝑑𝑝 𝐵𝑒 (4–4) −𝑑𝑉1=𝑉1𝑑𝑝 𝐵0 (4–5) −𝑑𝑉2=𝑉2𝑑𝑝 𝐵0 (4–6) 𝑑𝑉𝑐1=𝑉𝑐1𝑑𝑝 𝐵𝑐1 (4–7) 𝑑𝑉𝑐2=𝑉𝑐2𝑑𝑝 𝐵𝑐2 (4–8) En las ecuaciones (4–7) y (4–8), 𝐵𝑐𝑖 hace referencia a la flexibilidad de los contendores. 43 Herramienta de simulación computacional de maquinaria hidráulica Para contenedores metálicos, es habitual: 𝐵𝑐1≅𝑡𝐸 𝑑 (4–9) Donde • 𝑡 es el espesor de la pared. • 𝑑 es el diámetro. • 𝐸 es el módulo elástico del material, que en el caso de acero es de 𝐸=210 𝐺𝑃𝑎. La ecuación (4–3) puede expresarse como: 𝑉𝑡𝑑𝑝 𝐵𝑒=𝑉1𝑑𝑝 𝐵0+𝑉2𝑑𝑝 𝐵0+𝑉𝑐1𝑑𝑝 𝐵𝑐1 +𝑉𝑐2𝑑𝑝 𝐵𝑐2 (4–10) Si tenemos en cuenta (4–2), queda: 1 𝐵𝑒=1 𝐵0+𝑉𝑐1𝑑𝑝 𝑉𝑇𝐵𝑐1+𝑉𝑐2𝑑𝑝 𝑉𝑇𝐵𝑐2 (4–11) Esta expresión se puede generalizar de la forma: 1 𝐵𝑒=1 𝐵0+∑ 1 𝐵𝑐𝑖𝑉𝑐𝑖 𝑉𝑇 𝑛 𝑖=1 (4–12) El módulo de compresibilidad efectivo para las mangueras es de los que más influencia tiene dada la flexibilidad de estas. 4.2 Teoría de las propiedades concentradas Esta teoría es muy usada para modelar sistemas hidráulicos. Para poder aplicarla los efectos de las ondas de presión han de ser insignificantes, para lo que es necesario que los tramos de las canalizaciones sean cortos, lo cual se cumple puesto que es una máquina; y que los ciclos de trabajo sean lentos, que también se cumple ya que los movimientos de una retroexcavadora no se caracterizan por su alta velocidad. Básicamente, la teoría de las propiedades concentradas divide el sistema hidráulico en volúmenes donde la presión se puede suponer constante, de modo que se puede plantear una ecuación diferencial para cada volumen. Estos volúmenes estarán separados por válvulas, cuyo caudal se calcula por la diferencia de presiones entre uno y otro volumen. Actuadores hidráulicos 44 4.2.1 Ecuación de continuidad Sea un volumen de control al que entra un caudal 𝑄𝑖𝑛 y sale un caudal 𝑄𝑜𝑢𝑡. Se puede expresar: 𝑚󰇗𝑖𝑛=𝜌𝑖𝑛𝑄𝑖𝑛 (4–13) 𝑚󰇗𝑜𝑢𝑡=𝜌𝑜𝑢𝑡𝑄𝑜𝑢𝑡 (4–14) Figura 4-3 Balance de materia La masa del volumen de control es: 𝑀=𝜌𝑉 (4–15) Por conservación de la materia, la única forma de que varíe la masa en el volumen de control está en que los caudales másicos a la entrada y salida sean diferentes: 𝑑𝑀 𝑑𝑡=𝑚󰇗𝑖𝑛−𝑚󰇗𝑜𝑢𝑡 (4–16) Teniendo todo esto en cuenta, y suponiendo la densidad igual a la salida y entrada: 𝑑𝜌 𝑑𝑡𝑉 𝜌+𝑑𝑉 𝑑𝑡=𝑄𝑖𝑛−𝑄𝑜𝑢𝑡 (4–17) Suponiendo que la masa en el volumen de control es también constante: 𝑑𝑀=0=𝜌𝑑𝑉+𝑉𝑑𝜌 (4–18) Esta ecuación lleva a: 𝑑𝜌 𝜌=−𝑑𝑉 𝑉 (4–19) Atendiendo a la ecuación (4–1): 𝑑𝜌=𝑑𝑝𝜌 𝐵𝑒 (4–20) Sustituyendo (4–20) en (4–17) se llega a la ecuación: 𝑑𝑝 𝑑𝑡=𝐵𝑒 𝑉[𝑄𝑖𝑛−𝑄𝑜𝑢𝑡−𝑑𝑉 𝑑𝑡] (4–21) Habitualmente, el término 𝑑𝑉 𝑑𝑡 se refiere al movimiento del pistón del cilindro, que será función de las coordenadas del sistema sobre el que actúe. 45 Herramienta de simulación computacional de maquinaria hidráulica 4.2.2 Relación flujo -presión En la ecuación (4–21), es necesario expresar los términos 𝑄𝑖𝑛 y 𝑄𝑜𝑢𝑡 en función de la presión. Sea un flujo turbulento (𝑅𝑒>2400), en el que las pérdidas de presión crecen exponencialmente con la velocidad, el caudal a través de una válvula sigue la ecuación del orificio: 𝑄=𝐶𝑑𝐴 √2 𝜌𝑑𝑝 (4–22) Donde 𝐶𝑑 es el coeficiente de descarga. De forma semi-empírica, se puede expresar: 𝑄=𝐶𝑣√|Δ𝑝|Δ𝑝 |Δ𝑝| (4–23) Donde 𝐶𝑣=𝐶𝑑𝐴√2 𝜌, y Δ𝑝 |Δ𝑝| indica el sentido del caudal, hacia la zona de menor presión. Sin embargo, a pesar de que se suponga un flujo turbulento, en el caso de diferencias de presiones cercanas a cero, la pendiente de la curva en la ecuación (4–23) tiende a infinito. Por ello se hace necesario modelar la parte laminar, de evolución lineal, para evitar problemas numéricos. Para ello se usará un spline de manera que se eviten discontinuidades en la curva, lo que también implicaría errores numéricos (Figura 4-4). El polinomio interpolador se ha calculado teniendo en cuenta que en 𝑅𝑒=2400, el punto y la pendiente son las del equivalente siguiendo la ecuación (4–23); y que, para una diferencia de presiones igual a cero, el caudal es cero con la pendiente de la recta que une el origen con ese punto de transición laminar-turbulento. Figura 4-4 Caudal en función de la diferencia de presiones Actuadores hidráulicos 52 4.5 Parámetros Atendiendo al modelo de la excavadora [2] y a catálogos comerciales de válvulas, Parker [10], se utilizan los siguientes parámetros para el circuito hidráulico. Hay que destacar que son valores orientativos que mantienen el orden de magnitud. Habría que tener acceso a datos mucho más específicos de las características de la máquina y del sistema hidráulico. 4.5.1 General En cuanto a los parámetros de carácter general destacar las presiones de bomba (4–48) y tanque (4–49), el diámetro (4–50) y módulo de compresibilidad (4–51) de las mangueras, el módulo de compresibilidad de las cámaras del cilindro hidráulico (4–52), y el número de Reynolds crítico (4–53) mencionado en 4.2.2 𝑝𝑝=36 𝑀𝑃𝑎 (4–48) 𝑝𝑇=0.1 𝑀𝑃𝑎 (4–49) 𝐷ℎ=0.01 𝑚 (4–50) 𝐵ℎ=350 𝑀𝑃𝑎 (4–51) 𝐵0= 21000 𝑀𝑃𝑎 (4–52) 𝑅𝑒𝑐𝑟𝑖𝑡=2400 (4–53) 4.5.2 Aceite hidráulico Para el aceite del circuito se han elegido unas propiedades típicas de este tipo de fluidos: 𝜌=870𝑘𝑔 𝑚3 (4–54) 𝐵=1860 𝑀𝑃𝑎 (4–55) 𝜈=32𝑚2 𝑠 (4–56) 4.5.3 Válvulas Como se ha comentado en el apartado 4.2.2, las válvulas se modelan con la ecuación del orificio. Atendiendo al caudal que circula por la máquina, las presiones típicas, y los catálogos, se tiene: 𝐶𝑣1=5∗10−6 𝑚3 𝑠√𝑃𝑎 (4–57) 𝐶𝑣2=2∗10−6 𝑚3 𝑠√𝑃𝑎 (4–58) 53 Herramienta de simulación computacional de maquinaria hidráulica 4.5.4 Cilindros Se detallan las extensiones máximas, los diámetros de las cámaras y las longitudes de las mangueras estimadas para los cilindros de la retroexcavadora. 4.5.4.1 Pluma ℎ=1.7 𝑚 (4–59) 𝐿ℎ1=4 𝑚 𝐿ℎ2=2 𝑚 𝐿ℎ3=8 𝑚 (4–60) 𝐷𝑝𝑙𝑢𝑚𝑎=0.3 𝑚 (4–61) 𝐷𝑟𝑜𝑑=0.2 𝑚 (4–62) 𝜂=0.7 (4–63) 4.5.4.2 Balancín ℎ=1.7 𝑚 (4–64) 𝐿ℎ1=6 𝑚 𝐿ℎ2=4 𝑚 𝐿ℎ3=10 𝑚 (4–65) 𝐷𝑏𝑎𝑙𝑎𝑛𝑐í𝑛=0.2𝑚 (4–66) 𝐷𝑟𝑜𝑑=0.14 𝑚 (4–67) 𝜂=0.7 (4–68) 4.5.4.3 Cuchara ℎ=1.7 𝑚 (4–69) 𝐿ℎ1=8 𝑚 𝐿ℎ2=6 𝑚 𝐿ℎ3=14 𝑚 (4–70) 𝐷𝑐𝑢𝑐ℎ𝑎𝑟𝑎=0.2 𝑚 (4–71) 𝐷𝑟𝑜𝑑=0.14 𝑚 (4–72) 𝜂=0.7 (4–73) Actuadores hidráulicos 54 4.5.5 Actuador rotatorio Los parámetros del actuador rotatorio han sido escogidos con el orden de magnitud de [9]: 𝐶ℎ=8.125∗10−13𝑚3 𝑃𝑎 (4–74) 𝐾𝑡=1.85∗10−9𝑚3 𝑃𝑎 (4–75) 𝑝𝑠𝑒𝑡=30 𝑀𝑃𝑎 (4–76) 𝑉1 𝑚𝑎𝑥 =71 𝑐𝑚3 (4–77) 𝑛1=1680 𝑟.𝑝.𝑚. (4–78) 𝑥1 𝑚𝑎𝑥=43.62 𝑚𝑚 (4–79) 𝑉2=90 𝑐𝑚3 (4–80) 𝑝0=2 𝑀𝑃𝑎 (4–81) 55 5 SISTEMA DE CONTROL Se desarrolla ahora el modelo del controlador implementado para garantizar las posiciones del cilindro mediante la regulación de la válvula direccional de flujo en cada uno de los actuadores considerados. 5.1 Fundamento Una vez definidas las funciones que modelan el comportamiento físico de la máquina, se procede a introducir un controlador que regule la señal 𝑈 de las ecuaciones (4–35) a (4–37). Para ello habrá que definir una posición deseada de la extensión del cilindro, y mediante una realimentación negativa en lazo cerrado, el controlador emitirá la señal 𝑈 adecuada al sistema. Figura 5-1 Control con realimentación negativa Esta realimentación negativa, compara la posición deseada, 𝑟(𝑡), con la posición actual, 𝑥(𝑡), lo que da lugar a un error 𝑒(𝑡), en función del cual el controlador, 𝐶(𝑡), calcula la señal 𝑈(𝑡), de entrada al sistema, manejando así la extensión o retracción del cilindro. La razón por la que se pretende controlar la extensión del cilindro en lugar de las propias coordenadas del sistema se debe a la forma de manejar la máquina. El operario que trabaja desde la cabina utiliza unos mandos con los que controla los diferentes actuadores. Cuando se actúa sobre los mandos, se está realizando un cambio en las posiciones de referencia de los cilindros, por lo que se producirá una actuación hasta que se alcance la nueva posición deseada. En el caso de uso normal de una retroexcavadora no se sabe con certeza el ángulo que ha girado un sólido con respecto al otro, lo que interesa es saber si hace falta inclinar más o menos la pala. Por ello el control directamente sobre el cilindro hidráulico supone un mayor acercamiento a la realidad. Además, de cara a posibles problemas computacionales, es necesario evitar que la extensión del cilindro supere los límites físicos del propio cilindro, por lo que actuando directamente sobre él se tiene un mayor control de este posible problema. Sistema de control 56 5.2 Referencia Es importante tener en cuenta que, si la señal de referencia entrara como un escalón, ese cambio brusco inicial de la señal produciría una sobreoscilación que haría que, aunque la posición de referencia se encuentre en los límites físicos, el pico de la sobreoscilación los supere. Por tanto, se ha optado por modelar la posición de referencia, así como su velocidad y aceleración para suavizar la señal que será emitida y mitigar en cierta forma ese efecto indeseado. De este modo, la señal de referencia se formaría tal y como se muestra en la Figura 5-2. Se le asigna un tiempo de funcionamiento a la señal para evitar ese impulso. Figura 5-2 Señal de referencia 5.3 Controlador PD Una vez definida la referencia, se procede ahora a definir el controlador. La parte integral del controlador para minimizar el posible error en régimen permanente no se ha tenido en cuenta por diversos motivos. El primero es que, debido al bajo rango de movimiento respecto a cualquier posición inicial, el error que se produce con un controlador proporcional derivativo (PD), es despreciable. Debido a esto, el incluir la parte integral incrementaría el coste computacional a cambio de una muy pequeña mejora. El segundo motivo se debe a la propia operatividad de la máquina. El operario, si necesita una mayor amplitud en la extensión del cilindro, únicamente necesita aumentar más la misma con los mandos. Por ello no es necesario un control de extrema precisión en la excavadora. De este modo, definimos: 𝑒(𝑡)=𝑟(𝑡)−𝑥(𝑡) (5–1) 𝑈(𝑡)=𝐾𝑝(𝑒(𝑡)+𝑇𝑑𝑑𝑒(𝑡) 𝑑𝑡 ) (5–2) De este modo, la ecuación (5–2) es la ley de control 𝐶(𝑡) definida en la Figura 5-1. 𝐾𝑝 es la ganancia de la parte proporcional del controlador y hay que elegirla adecuadamente para cada uno de los cilindros. Un aumento en la ganancia implica aumentar la velocidad del sistema, ya que se amplifica el error y la señal que llega al sistema, pero puede inestabilizar el resultado. Por su parte 𝑇𝑑 es el tiempo derivativo, que aumenta la estabilidad relativa, es decir, reduce las sobreoscilaciones del sistema, pero puede llegar a inestabilizar el sistema si 𝑇𝑑 es demasiado grande. 57 Herramienta de simulación computacional de maquinaria hidráulica También es habitual expresar la ley de control (5–2) de otra forma: 𝑈(𝑡)=𝐾𝑝𝑒(𝑡)+𝐾𝑑𝑑𝑒(𝑡) 𝑑𝑡 𝐾𝑑=𝐾𝑝 𝑇𝑑 (5–3) En el caso del actuador rotatorio, por el modelo de la bibliografía [9], se ha controlado solamente en velocidad. Además, es la única forma de la que se ha conseguido un control estable del mismo. 𝑥𝑠=𝐾𝑑𝑑𝑒(𝑡) 𝑑𝑡 (5–4) 5.4 Sintonización Debido a la complejidad del sistema estudiado, la sintonización de los parámetros del controlador para que funcione de manera rápida y correcta, evitando la inestabilidad, es una tarea laboriosa. Se deben usar métodos heurísticos para obtener parámetros adecuados. La manera de proceder para la sintonización del controlador es parecida a la del método de Zieger-Nichols. Se ha ido sintonizando cada uno de los controladores realizando simulaciones con la ganancia derivativa 𝐾𝑑 igualada a cero, e incrementando poco a poco 𝐾𝑝 hasta que se inestabilizara el sistema. Cuando esto ocurría se reducía a la mitad esa ganancia proporcional y se realizaba el mismo proceso con la ganancia derivativa hasta obtener el mejor resultado posible. A pesar de ello, el método no ha sido demasiado útil, y únicamente se ha usado como referencia para analizar el orden de magnitud de los parámetros del controlador, siendo realizada manualmente. En los ensayos realizados en el capítulo 7, se indicará con qué parámetros está actuando el controlador para los ensayos realizado. Sistema de control 58 59 6 CÁLCULO En este capítulo se va a explicar cómo se van a unir la dinámica de sistemas multicuerpo, y la hidráulica en la integración. 6.1 Sistema de ecuaciones diferenciales Primero se van a recuperar las ecuaciones diferenciales desarrolladas en capítulos anteriores. En primer lugar, del capítulo 3 se obtuvo: 𝐪󰇘=𝐌−1(𝐐+𝐂𝐪 𝑇(𝐂𝐪𝐌−1𝐂𝐪 𝑇)−1(𝐐𝑑−𝐂𝐪𝐌−1𝐐)) (6–1) El vector 𝐐 de fuerzas generalizadas se puede definir como 𝐐=𝐐(𝐪,𝐪󰇗,𝐩,𝑡). En primer lugar, la dependencia de 𝐪,𝐪󰇗, se explica directamente con el término 𝐐𝑣(𝐪,𝐪󰇗), dependiente de la posición y la velocidad. Por su parte, 𝐐𝑔𝑟𝑎𝑣𝑣(𝐪), solo depende de la posición. La matriz 𝐂𝐪(𝐪), 𝐌(𝐪), y 𝐐𝑑(𝐪,𝐪󰇗), también dependen de esas variables. En (6–1), se ha incluido una posible dependencia explícita del tiempo, ya que la señal 𝑈(𝑡) va a variar con el tiempo 𝑡 de una manera concreta. El vector 𝐩, es un vector que contiene las 3 presiones de cada uno de los 3 cilindros, más la presión del actuador rotatorio para el giro de la cabina. La dependencia de 𝐐 de estas variables se basa en: 𝐐𝑎𝑝𝑙=𝐋𝑇𝐅ℎ𝑦𝑑 (6–2) Donde 𝐅ℎ𝑦𝑑(𝐪,𝐪󰇗,𝐩) constituye el vector de fuerzas y momentos resultantes obtenidos de la ecuación (4–38) y (4–60), para cada uno de los actuadores descritos. Teniendo esto en cuenta, se puede decir por tanto que: 𝐪󰇘=𝒇(𝐪,𝐪󰇗,𝐩,𝑡) (6–3) Además, de las ecuaciones (4–31) a (4–33), se concluye que: 𝐩󰇗=𝒉(𝐪,𝐪󰇗,𝐩,𝑡) (6–4) Si ahora se llama: 𝐪󰇗=𝐯 (6–5) Esto implica que: 𝐯󰇗=𝐪󰇘=𝒇(𝐪,𝐯,,𝐩,𝑡) (6–6) 𝐩󰇗=𝒉(𝐪,𝐯,𝐩,𝑡) (6–7) Cálculo 60 De este modo, se plantea un sistema de ecuaciones diferenciales no lineales, en el que las incógnitas son los vectores 𝐪, 𝐯, y 𝐩. Las ecuaciones diferenciales de (6–5) a (6–7), se reducen a: 𝐗󰇗=𝐅(𝐗,𝑡) (6–8) Donde: 𝐗=[𝐪 𝐯 𝐩] (6–9) 𝐗󰇗=[𝐪󰇗𝐯󰇗𝐩󰇗] (6–10) 𝐅(𝐗,𝑡)=[ 𝐯 𝒇(𝐪,𝐯,𝐩,𝑡) 𝒉(𝐪,𝐯,𝐩,𝑡)] (6–11) Este problema es mucho más abordable computacionalmente que el problema inicial y es la manera en que se va a implementar en el programa de cálculo mediante algún algoritmo de integración. Destacar que, al realizar simulaciones en bucle cerrado con un controlador, la forma del sistema no cambia. Lo único que varía es la señal 𝑈(𝑡) según una ley conocida que va cambiando en función del tiempo y de la posición. 6.2 Condiciones iniciales En este sistema de ecuaciones, se necesita incluir unas condiciones iniciales que impliquen que se cumplan tanto las restricciones, como el equilibrio de fuerzas. Además, debe ser una posición de reposo. Por ello: 𝐅(𝐗0,0)=𝟎 (6–12) Que no es más que resolver un sistema de ecuaciones algebraicas no lineales formada por: [𝒇(𝐪0,𝟎,𝐩0,0) 𝒉(𝐪0,𝟎,𝐪0,0)]=𝟎 (6–13) El sistema (6–13) tiene infinitas soluciones dependiendo de donde se parta para resolver el sistema de forma iterativa. Incluyendo una posición inicial orientativa a la deseada, y unas presiones del orden de magnitud apropiado, el sistema llega a una solución en la que las presiones iniciales se ajustan para mantener esa posición inicial del brazo de la excavadora. Una vez obtenida la posición de equilibrio, ya tenemos el punto de partida para llevar a cabo todas las simulaciones deseadas, tanto en bucle abierto como en bucle cerrado. 61 Herramienta de simulación computacional de maquinaria hidráulica 6.3 Programa MatLab En esta sección se van a explicar de manera breve los programas principales desarrollados en MatLab para la simulación de la máquina. 6.3.1 “Simbolico.m” Este archivo contiene el programa simbólico que permite desarrollar las ecuaciones en función de parámetros y variables declarados de manera simbólica. De esta manera, las derivadas, los jacobianos y los productos de matrices quedan en función de las variables en lugar de realizarse la aproximación numérica de los mismos. Una vez realizadas todas las operaciones es posible generar funciones de MatLab que reciban los datos como entrada. Este programa es de compilación muy lenta, pero sólo hay que ejecutarlo una vez para generar las funciones, agilizando la posterior integración de las ecuaciones. De la ejecución de este archivo, obtenemos: • “pos###.m”: es una función que proporciona la posición del punto indicado en “###” en función de las variables y los parámetros. • “jacobiano###.m”: determina el jacobiano del punto “###” en función de las variables y los parámetros. • “poscilindropluma.m” y “velcilindropluma.m”: esta función determina el vector que une los extremos del cilindro de la pluma y la diferencia de velocidad entre esos puntos, respectivamente. • “poscilindrobalancin.m” y “velcilindrobalancin.m”: esta función determina el vector que une los extremos del cilindro del balancín y la diferencia de velocidad entre esos puntos, respectivamente. • “poscilindrocuchara.m” y “velcilindrocuchara.m”: esta función determina el vector que une los extremos del cilindro de la cuchara y la diferencia de velocidad entre esos puntos, respectivamente. • “matrizL.m”: en función de las variables y los parámetros geométricos proporciona la matriz 𝐋. • “matrizl_.m”: en función de las variables, sus derivadas y los parámetros geométricos proporciona la matriz 𝐥. • “matrizMgorro.m”: es una función de la que se obtiene la matriz 𝐌 , en función de los parámetros másicos del sistema. • “Qgravgorro.m”: se obtiene el sumando correspondiente a la gravedad del vector 𝐐 𝑎𝑝𝑙, de fuerzas aplicadas, en función de las masas de los sólidos y de la gravedad. • “Qvgorro.m”: es una función de la que se obtiene 𝑸 𝑣 con las variables, sus derivadas y los parámetros geométricos y másicos de los sólidos como entrada. • “fuerzaaplicadagorro.m”: de esta función se obtiene el sumando del vector de fuerza 𝐐 𝑎𝑝𝑙 debido a los cilindros aplicados a los sólidos, para lo que se le introducen las variables, los parámetros geométricos y las fuerzas resultantes de la actuación de cada cilindro, que provienen de otros programas. 6.3.2 “Restriccionessymb.m” Al igual que en el apartado anterior, se hace uso del lenguaje simbólico para obtener el vector de restricciones y las diferentes derivadas consideradas en la sección 3.2.1. Las funciones obtenidas de este archivo son: • “restricciones.m”: de aquí se obtiene el valor de las ecuaciones que definen el vector 𝐂 de restricciones, para lo que es necesario dar como entradas las coordenadas generalizadas y los parámetros geométricos del sistema. • “jacobiano.m”: en este caso obtenemos la matriz 𝐂𝐪 en función también de las coordenadas generalizadas y de los parámetros geométricos. • “jacobianopunto.m”: por último, se obtiene la función 𝐂󰇗𝐪=(𝐂𝐪𝐪󰇗)𝐪, la derivada respecto al tiempo del jacobiano de las restricciones, para lo que se requiere como entrada, además de los datos anteriores, la derivada de las coordenadas generalizadas. Resultados 68 incluida en las condiciones iniciales de equilibrio para provocar el movimiento del sistema. Por su parte, en la Figura 7-5 se aprecia el cumplimiento de las restricciones asociadas al sólido 2, la base. En este caso el cumplimiento es mucho más preciso, ya que no se ha producido una perturbación en sus condiciones iniciales. Todos estos resultados se han obtenido con unos coeficientes de Baumgarte de 𝛼=𝛽=100. Figura 7-5 Restricciones de la base Si aumentamos en un orden de magnitud los coeficientes de estabilización de Baumgarte (𝛼=𝛽=1000), vemos que las restricciones del mecanismo de 4 barras, aunque inicialmente existe esa perturbación inicial, en general se cumplen también con un orden de magnitud mayor de precisión, del orden de 10−3. Figura 7-6 Restricciones del brazo con mayor valor de los coeficientes de Baumgarte 69 Herramienta de simulación computacional de maquinaria hidráulica 7.2 Bucle cerrado El objetivo de esta sección es comprobar el buen funcionamiento del simulador, e ilustrar la cantidad de resultados y análisis que se puede hacer a partir de ellos. 7.2.1 Condición inicial Se ha realizado una simulación en bucle cerrado de todas las variables controladas, para realizar un análisis de los resultados obtenidos. Se parte de la posición de equilibrio: 𝐪0= [ 𝑥20 𝑦20 𝑧20 𝜑20 𝜃20 𝜓20 𝜓30 𝜃40 𝜃50 𝜃60 𝜃70 𝜃80 ] = [ 0 0 0 0 0 0 0 −0.5521 1.4518 0.0045 −1.3268 1.9714 ] 𝐩𝑝𝑙𝑢𝑚𝑎0=[15.54 15.54 25.73] 𝑀𝑃𝑎 𝐩𝑏𝑎𝑙𝑎𝑛𝑐í𝑛0=[13.90 13.90 29.46] 𝑀𝑃𝑎 𝐩𝑐𝑢𝑐ℎ𝑎𝑟𝑎0=[15.04 15.04 30.13] 𝑀𝑃𝑎 𝐩ℎ0=2 𝑀𝑃𝑎 (7–1) La posición de equilibrio determinada por 𝐪0 se corresponde con la Figura 7-7. Figura 7-7 Posición inicial Resultados 70 Por su parte, en los equilibrios de presiones de los tres cilindros, se ve cómo 𝑝1=𝑝2<𝑝3. Esto tiene sentido ya que no debe haber flujo a través de la válvula estranguladora, y para ello las presiones a la entrada y salida deben ser iguales. Por otro lado, la presión del lado de la barra del cilindro, al tener menor área, necesita más presión para que se produzca un equilibrio de fuerzas en los cilindros, por lo que el hecho de que su valor sea del orden del doble no es nada sorprendente. El valor inicial de 𝑝ℎ coincide con la presión del lado de baja presión del actuador rotatorio, de manera que inicialmente no se genera momento. Las posiciones de 𝑥 de los respectivos cilindros en el instante inicial es: 𝑥𝑝𝑙𝑢𝑚𝑎0=0.7307 𝑚 𝑥𝑏𝑎𝑙𝑎𝑛𝑐í𝑛0=1.5027 𝑚 𝑥𝑐𝑢𝑐ℎ𝑎𝑟𝑎0=0.7072 𝑚 (7–2) 7.2.2 Parámetros El objetivo de esta simulación es llevar las extensiones de los cilindros a: 𝑥𝑝𝑙𝑢𝑚𝑎𝑓=1.3000 𝑚 𝑥𝑏𝑎𝑙𝑎𝑛𝑐í𝑛𝑓=1.000 𝑚 𝑥𝑐𝑢𝑐ℎ𝑎𝑟𝑎𝑓=0.1000 𝑚 (7–3) Por su parte, el control en velocidad buscará que el giro alcance una velocidad determinada de giro para luego parar. Esta velocidad es: 𝜓3 󰇗= 0.4 𝑟𝑎𝑑 𝑠 (7–4) Para lograr esas posiciones finales los parámetros de control elegidos son: 𝐾𝑝𝑝𝑙𝑢𝑚𝑎=3 𝐾𝑑𝑝𝑙𝑢𝑚𝑎=0.2 𝐾𝑝𝑏𝑎𝑙𝑎𝑛𝑐í𝑛=1.5 𝐾𝑑𝑏𝑎𝑙𝑎𝑛𝑐í𝑛=0.08 𝐾𝑝𝑐𝑢𝑐ℎ𝑎𝑟𝑎=5 𝐾𝑑𝑐𝑢𝑐ℎ𝑎𝑟𝑎=0.2 𝐾𝑑=200 (7–5) Como se comentó en el Capítulo 5, la elección de estos parámetros es fruto de la prueba y error en el modelo, debido a la complejidad que entraña su control. 71 Herramienta de simulación computacional de maquinaria hidráulica 7.2.3 Cilindros Como resultado de esta simulación, obtenemos las siguientes evoluciones de la extensión de los cilindros y sus respectivas velocidades. En el caso del cilindro de la pluma, los resultados son: Figura 7-8 Extensión del cilindro de la pluma Figura 7-9 Velocidad de extensión del cilindro de la pluma Como se puede apreciar en las Figura 7-8 y Figura 7-9 la dinámica del cilindro evoluciona de manera muy lenta respecto a la referencia. Las referencias introducidas corresponden a velocidades del orden de las de trabajo habitual de una retroexcavadora, pero la pluma se demora unos 5 segundos en alcanzar la posición deseada. Esto se aprecia también en la velocidad, que no llega a la velocidad establecida por la referencia. La solución a este problema se encontraría en aumentar las constantes del controlador, de manera que las variables evolucionaran de manera más rápida hacia la referencia. El inconveniente es la inestabilidad que aparece en el régimen permanente a partir de ciertos valores de 𝐾𝑝 y 𝐾𝑑, habiendo sido seleccionados los máximos valores posibles. Es importante indicar que, a pesar del retardo de la respuesta, la posición se alcanza perfectamente sin error alguno en régimen permanente, justificando la elección de un controlador sin parte integral. La evolución de las presiones (Figura 7-10) está relacionada con la de la velocidad. Recordar, que el término 𝑝1, es la presión antes de la válvula estranguladora, 𝑝2 es la presión en la sección mayor del pistón del cilindro, y 𝑝3 es la presión en la cámara de menor sección. Resultados 72 Inicialmente, el controlador abre la válvula direccional de manera que pueda entrar aceite en el cilindro y se produzca la extensión. Por ello, la presión antes del estrangulamiento aumenta rápidamente a valores cercanos al de la bomba (36 𝑀𝑃𝑎), por estar conectada directamente a ella, provocando el flujo a través de la válvula. La diferencia entre las presiones a ambos lados del cilindro se va haciendo más pequeña conforme aumenta la velocidad, de modo que la fuerza resultante favorezca la extensión del cilindro. Cuando la velocidad tiende a cero, entorno al segundo 7, las presiones han alcanzado un nuevo equilibrio para la posición deseada, de modo que 𝑝1=𝑝2, evitando que haya flujo en la válvula estranguladora; y 𝑝2 y 𝑝3 tienen un valor tal que la fuerza resultante es nula en el cilindro. Figura 7-10 Evolución de las presiones en el cilindro de la pluma En la Figura 7-11 y Figura 7-12, se muestra la retracción del cilindro del balancín. Aquí se puede comprobar el mismo problema del retardo. No obstante, en este caso la velocidad de respuesta es mayor, siendo menor la diferencia entre referencia y extensión real. Esto se debe a que hay menos masa sobre la que actuar con este cilindro, y se permite mayor velocidad sin llegar a la inestabilidad. De nuevo, el error en régimen permanente es prácticamente cero. Figura 7-11 Extensión del cilindro del balancín 73 Herramienta de simulación computacional de maquinaria hidráulica Figura 7-12 Velocidad de extensión del cilindro del balancín En la Figura 7-13 se puede ver la evolución de las presiones en la retracción del cilindro del balancín. Como en este caso el cilindro se está retrayendo, la cámara de presión 𝑝3 es la que se comunica directamente con la bomba a través de la válvula direccional, por lo que es la que alcanza las altas presiones. Al estar siendo comprimida la cámara de presión 𝑝2, esta aumenta, siendo mayor que 𝑝1, comunicada con el tanque, lo que explica su bajada inicial. Finalmente, tiende a igualarse a 𝑝2 cuando deja de haber flujo de aceite. Figura 7-13 Evolución de las presiones en el cilindro del balancín Resultados 74 En la retracción del cilindro de la cuchara ilustrado en la Figura 7-14 y Figura 7-15, se sigue la evolución de los dos cilindros anteriores. Al ser el que menos masa ha de poner en movimiento, la respuesta tiene mucha más velocidad. La velocidad de retracción prácticamente es capaz de alcanzar la velocidad de referencia, y el retardo en alcanzar la posición es de apenas 1 segundo. Igual que en el resto de los cilindros, el error es nulo en régimen permanente. Figura 7-14 Extensión del cilindro de la cuchara Figura 7-15 Velocidad de extensión del cilindro de la cuchara La evolución de presiones en la retracción de la cuchara se puede observar en la Figura 7-16. Al igual que el cilindro del balancín, al ser una retracción tiene el mismo comportamiento. En este cilindro, no obstante, se puede apreciar mucho mejor la estrecha relación entre la velocidad de extensión de los cilindros y la evolución de las presiones. El pico en el segundo 3 apreciable en la velocidad, también se aprecia en ese mismo instante en las presiones, cuando las presiones tienden al nuevo estado de equilibrio del final del movimiento. Esto se debe a que, como se explicó en la sección 4.2.2, la diferencia de presiones y el caudal a través de una válvula, se relacionan con una expresión del tipo 𝑄∝√Δ𝑝. Una consecuencia de esto se aprecia directamente en las gráficas: si la diferencia de presiones entorno a la válvula es cero, no habrá caudal a través de ella, y, por tanto, no habrá movimiento del pistón. Cuando 𝑝1=𝑝2, se para la extensión o retracción del cilindro. 75 Herramienta de simulación computacional de maquinaria hidráulica Figura 7-16 Evolución de las presiones en el cilindro de la cuchara Por último, el control en velocidad del actuador rotatorio ha resultado ser el más efectivo como se puede apreciar en la Figura 7-17, ya que permite elevar muchísimo 𝐾𝑑, sin inestabilizar el sistema, por lo que referencia y velocidad real son prácticamente idénticas a excepción del inicio. En la Figura 7-18 también se muestra el giro de la cabina correspondiente a 𝜓3 como consecuencia del actuador rotatorio. Aunque se muestre la referencia de posición para el giro, como para el control del actuador rotatorio solo se ha usado control en velocidad, no es de gran interés de análisis. Sin embargo, se puede explicar el error en régimen permanente, ya que no hay ningún factor que obligue a que el giro sea el mismo que la referencia, y esa diferencia inicial en las velocidades origina esa diferencia en la posición. Notar que la máquina ha girado casi 180 𝑜 en 10 𝑠, que indica una velocidad de giro cercana a la de operación. Figura 7-17 Velocidad giro de la cabina Resultados 76 Figura 7-18 Giro de la cabina Si se hace un zoom sobre la parte superior de la Figura 7-17, podemos ver unas oscilaciones (Figura 7-19) que se van amortiguando, que tendrán sus consecuencias sobre la presión 𝑝ℎ (Figura 7-20) que provoca el par de giro (Figura 7-31). Estas oscilaciones son muy pequeñas, pero los picos de presión serán relativamente grandes. Figura 7-19 Oscilaciones de la velocidad de giro de la cabina Si se representa la presión del actuador rotatorio de la línea de alta presión (𝑝ℎ), el resultado es la Figura 7-20. En el instante en que la velocidad de giro alcanza el valor máximo y oscila de manera amortiguada, la presión sufre unas oscilaciones de las mismas características del orden de 10 𝑀𝑃𝑎 que se van amortiguando. Igual sucede cuando la velocidad cae a cero para detener el giro. Aparecen oscilaciones en la presión, aunque de menor amplitud, y en la velocidad de giro apenas se notan. 77 Herramienta de simulación computacional de maquinaria hidráulica Inicialmente la presión del actuador rotatorio aumenta hasta que se alcanza una velocidad estable, entonces cae, y se mantiene oscilando en torno a un valor cercano a cero, ya que la velocidad no varía. La presión se hace negativa para hacer que la velocidad disminuya hasta cero. El hecho de que la presión sea negativa, sólo quiere decir que la presión de la línea de alta se hace menor que la de baja para frenar la rotación, ya que realmente interesa la diferencia de presiones. El problema radica en que al tener que frenar el actuador, los papeles de las dos líneas se invierten, habiendo mayor presión en el otro ramal. Esto requeriría un modelado del actuador mucho más elaborado que tuviera en cuenta ese cambio. Figura 7-20 Evolución de la presión del actuador rotatorio 7.2.4 Controlador Relacionado con la sección anterior, en la Figura 7-21 se muestra el desplazamiento de la señal 𝑈(𝑡) y 𝑥𝑠(𝑡), que indica la posición de la válvula direccional del cilindro y del actuador rotatorio, respectivamente. Figura 7-21 Evolución de las señales de control de los cilindros Resultados 84 Por su parte, el momento ejercido por el actuador rotatorio, al ser directamente proporcional a la presión de la línea de alta presión, tiene su misma forma (Figura 7-31). Este momento se hace próximo a cero cuando la velocidad de giro se hace igual a la de referencia, ya que como no hay amortiguamiento considerado en el modelo de este elemento, no hace falta aplicar una fuerza que acelere el giro del sólido. Igual que 𝑝ℎ, 𝑀 se hace negativo para decelerar el giro. Figura 7-31 Momento resultante del actuador rotatorio Si se analizan ahora las fuerzas generalizadas, es decir, las que actúan directamente sobre las coordenadas generalizadas del sistema, se demuestra que las fuerzas cuadráticas en velocidad 𝑄𝑣, son de un orden de magnitud menor que las fuerzas aplicadas, 𝑄𝑎𝑝𝑙. Como ejemplo de esto se ha graficado las fuerzas aplicadas sobre la coordenada 𝜃4 (Figura 7-32). Ello quiere decir que, a las velocidades de trabajo normales de la retroexcavadora, las fuerzas cuadráticas en velocidad serían despreciables. Otra consecuencia de esto es que, como ya se comentó antes, la fuerza resultante de los cilindros en el equilibrio sólo debe vencer la fuerza gravitatoria. Figura 7-32 Fuerzas generalizada aplicadas sobre la coordenada 𝜃4 85 Herramienta de simulación computacional de maquinaria hidráulica Utilizando esta misma coordenada, se verifica que las fuerzas generalizadas resultantes (Figura 7-33) en el inicio y en el final sean nulas, para que el sistema esté en equilibrio. Además, es claro que inicialmente disminuye la fuerza para provocar la elevación de la pluma, aumentando finalmente la fuerza para parar el sólido. Figura 7-33 Fuerza generalizada resultante sobre la coordenada 𝜃4 7.2.7 Restricciones Para acabar de analizar la simulación en bucle cerrado, se muestra el cumplimiento de las restricciones y las fuerzas de reacción que aparecen debido a esas restricciones. Los coeficientes de Baumgarte empleados para garantizar el cumplimiento de las restricciones son 𝛼=𝛽=100. Empezando por el brazo, las restricciones correspondientes se muestran en la Figura 7-34. Se aprecia el buen cumplimiento de las restricciones, ya que los movimientos son relativamente lentos, y solo fuerzan ligeramente las restricciones durante el movimiento del cilindro de la cuchara. Figura 7-34 Restricciones del brazo Resultados 86 Las fuerzas necesarias para que se cumplan se muestran en la Figura 7-35. Como todas las fuerzas mencionadas, finalmente se establece un valor de equilibrio para garantizar el par de revolución. Figura 7-35 Fuerza de reacción en el par de revolución Por otro lado, se analizan las restricciones que hacen que la base esté inmóvil. Como ya ocurrió en la sección 7.1.2, las restricciones de la base se cumplen de manera más exacta que las del brazo, debido a que están menos solicitadas por el accionamiento. Figura 7-36 Restricciones de la base 87 Herramienta de simulación computacional de maquinaria hidráulica Las fuerzas de reacción se pueden relacionar directamente a las 6 coordenadas del sólido 2, ya que cada una de sus restricciones fijaba una de sus coordenadas. El resultado se muestra en la Figura 7-37 y la Figura 7-38. Estas fuerzas son las que se transmitirán al terreno. La tercera coordenada, es la mayor, ya que tiene que absorber todo el peso de la retroexcavadora completa. Si se calcula el peso total, la fuerza gravitatoria sería de 356.3 𝑘𝑁, justo la cantidad inicial de 𝛌(3). Las otras dos sólo tienen que absorber la componente axial de la fuerza ejercida por los cilindros. Figura 7-37 Fuerzas de reacción en la base Por otro lado, se ha representado también los momentos de reacción, fruto de fijar los ángulos de Euler de la base. 𝛌(6), correspondiente al momento respecto al eje vertical, tiene como resultado una evolución es la misma que el momento de giro aplicado por el actuador rotatorio. Las otras dos componentes simplemente anulan el posible momento aplicado en los otros dos ejes. Figura 7-38 Momentos de reacción en la base Resultados 88 89 8 CONCLUSIONES En este capítulo se van a exponer las conclusiones derivadas de los resultados del capítulo 7, y se van a plantear posibles desarrollos futuros de este trabajo. 8.1 Conclusiones De los resultados obtenidos, se puede concluir que el modelo matemático y la programación tienen la solidez necesaria para poder analizar los movimientos de una retroexcavadora, junto con la evolución de las presiones del circuito hidráulico. Se pueden extraer gran cantidad de resultados del modelo, como presiones, trayectorias, evoluciones temporales o posiciones de equilibrio, así como las aceleraciones, las fuerzas aplicadas o las fuerzas de reacción. Con el modelo se pueden cambiar diferentes parámetros y comprobar la influencia de los mismos. En este caso se podrían modificar los parámetros hidráulicos para ver los efectos del accionamiento en la excavadora modelada. A pesar de las simplificaciones realizadas los resultados obtenidos son muy coherentes y tienen sentido físico en su comportamiento. La programación del controlador ha dado lugar a resultados de evolución demasiado lenta para lo esperable en una máquina de este tipo. No obstante, la respuesta en bucle cerrado tiene unos tiempos dentro del orden de magnitud, motivo por el que se han incluido esos resultados. Este problema es común como se puede ver en la bibliografía. Por ejemplo, en P. H. Chang y S. J. Lee [11] se habla sobre la enorme dificultad del control de este tipo de máquinas, en gran parte debido a su gran tamaño y grandes masas. Por tanto, el sistema se caracteriza por ser fuertemente no lineal, tanto por el cambio de posición, como por las no linealidades de la hidráulica. También se comenta que la no linealidad aumenta con el tamaño de la máquina, por manejar mayores masas y mayores presiones. En el caso estudiado, la retroexcavadora es de gran tamaño, lo que explica esa lenta evolución. La mayoría de los estudios realizados sobre el control de estas máquinas, se han realizado en bajas condiciones de velocidad de operación, lejos de la velocidad real esperada. Destacar del modelo su adaptabilidad a otros tamaños de máquina. Se pueden cambiar únicamente las masas y geometría en el archivo de los datos, por lo que se convierte en una programación versátil. Además, es fácilmente reconvertible a otro tipo de máquinas con pequeños cambios, modificando algunas cuestiones de restricciones o posición de puntos, pero el resto de la programación del simulador serían las mismas. Todo el desarrollo del trabajo se ha realizado para la actuación y movimiento de la máquina en vacío; sin embargo, para poder comprobar el funcionamiento de la máquina trabajando, sólo sería necesario aumentar los parámetros másicos de la cuchara cargada o la reacción del terreno sobre la misma en el caso de estudiar la excavación. Por último, el programa podría ser utilizado con fines académicos. En las clases se podría enseñar la evolución e interacción de los sólidos con las animaciones, y conseguir una mejor comprensión de la influencia de los parámetros. También tiene interés de cara a tener una aplicación directa de la hidráulica. Conclusiones 90 8.2 Desarrollo futuro Durante el desarrollo del trabajo, se han realizado diversas simplificaciones y no se han tenido en cuenta diversos factores, que podrían tener una influencia significativa en los resultados obtenidos. Por ello se plantean diversas vías de mejora que se podrían tener en cuenta en un futuro desarrollo de este software. 8.2.1 Esquemas de co-simulación Como se ha comentado numerosas veces en el trabajo, el modelo hidráulico constituye un sistema mucho más rígido que el de las ecuaciones del sistema multicuerpo. Por tanto, para una correcta integración del modelo global sería más coherente utilizar un paso de tiempo de integración mucho más pequeño en el sistema hidráulico que el multicuerpo. Es aquí donde entraría en juego la co-simulación. La bibliografía al respecto es numerosa, pero en este trabajo habría que destacar «Esquema de co-simulación del proceso de sinterización», D. García Vallejo [12], donde se desarrolla un procedimiento matemático para calcular la porosidad y la temperatura, siendo estas dos interdependientes. Utilizando un esquema predictor-corrector se realiza la simulación, con incrementos de tiempos de integración diferentes para porosidad y temperatura. A grosso modo, el predictor consiste en estimar la porosidad en un instante de tiempo, conocida en el instante anterior. Suponiendo la porosidad constante media entre esos dos instantes de tiempo, se puede calcular la temperatura. La parte del corrector ajusta la estimación de la porosidad hecha en el predictor usando la nueva temperatura, y de nuevo se calcula una porosidad media con la que calcular también una temperatura corregida. Se hacen tantas correcciones como sean necesarias para cumplir las tolerancias impuestas por el usuario. En este esquema de co-simulación, se utiliza un incremento de tiempo ℎ𝜃 para estimar la porosidad, y otro diferente ℎ𝑇 para estimar la temperatura. La idea sería aplicar algún esquema de este tipo a nuestro sistema para conseguir una integración más exacta y unos resultados más fiables. 8.2.2 Programación del integrador De igual forma que lo planteado en el apartado anterior, se podría estudiar la programación propia del bucle integrador del sistema. Esto implicaría un mayor control de las tolerancias los pasos de integración y, con ello, poder modificar la velocidad de integración del programa. El principal inconveniente de esto es la posible falta de exactitud. La rigidez de la hidráulica complica la integración de manera habitual, por lo que sería necesario ser capaz de modificar el paso de integración y ajustar las evaluaciones de manera más adaptativa. En este trabajo esta tarea no se ha realizado por exceder los límites de los objetivos fijados. Sería por tanto interesante plantear la programación del integrador utilizando algún algoritmo más sofisticado, como los propuestos por R. A. Ortega Aguilera, en «Formulación conservativa de restricciones en sistemas multicuerpo», [13]. Aquí, se explican otros métodos para la integración de sistemas con las ecuaciones diferenciales algebraicas como el método del coeficiente penalizador, o el método de Lagrange aumentado. Además, se incurre en métodos de integración como las iteraciones simultáneas e integraciones anidadas. Aprovechando estas iteraciones simultáneas se podría aprovechar para incluir un esquema de co-simulación como el visto anteriormente y tener un modelo de integración más robusto para un sistema tan complejo como este. 91 Herramienta de simulación computacional de maquinaria hidráulica 8.2.3 Leyes de control Se ha podido comprobar en las simulaciones en bucle cerrado que el controlar un sistema de este tipo es una tarea muy compleja, y en la que intervienen muchos aspectos que no se pueden tener en cuenta, y menos con la ley de control sencilla implementada. Además, la respuesta es más lenta a medida que intervienen más masas debido a una mayor inercia y mayores presiones. Por ello, para un mejor desempeño del programa, se podría utilizar leyes de control más complicadas, y que tengan en cuenta más factores, como el control TDC utilizado en [11] Relacionado con esto interviene la señal de referencia empleada. Se ha visto que los cilindros han sido controlados en posición, mientras que el actuador rotatorio se ha controlado en velocidad, ya que la referencia de la que se ha extraído el modelo recomendaba hacerlo así y el resultado es mucho mejor. Se podría intentar controlar los cilindros también en velocidad de manera más efectiva, y de manera predictiva, incluir el término derivativo de la velocidad, que es la aceleración. Esto sería sencillo de realizar si se hubiera programado un integrador propio, ya que se podría llamar a la función que contiene todo el sistema multicuerpo para calcular la aceleración dentro del bucle. El hecho de haber programado el integrador permitiría además introducir un término integral en la ley de control empleada, ya que se tendría acceso a los instantes de tiempo durante la integración, y se podrían realizar integraciones de tiempo parciales que cuantifiquen el error acumulado. Parece claro que uno de los siguientes pasos a realizar en la mejora del software es la programación de un integrador más adecuado a las necesidades. 8.2.4 Restricciones no holónomas Como se comentó en 2.3.1, la base móvil o chasis de la excavadora íbamos a considerarlo como estacionario. Parece lógico que el siguiente paso sea moverlo. Para ello habría que considerar la restricción no holónoma de la rodadura sin deslazamiento del sistema de desplazamiento de la máquina. Esto implicaría un tratamiento ligeramente diferente de las ecuaciones y la inclusión de más variables para definir el movimiento de cada una de las dos orugas. 8.2.5 Reducción de simplificaciones Como la transmisión se lleva a cabo a través de un motor de gasolina, se podría modelar el mismo y con ese modelo sacar la fuerza de transmisión aplicada. Además, es ese motor del que se saca la potencia para mover las bombas que impulsan el aceite al circuito hidráulico. Se podría por tanto eliminar la suposición de presión de la bomba constante y modelar la dependencia del motor de gasolina. Podría cogerse un modelo más exacto de máquina y con ello conseguir un modelo muy completo que tenga en cuenta muchos parámetros. 8.2.6 Joystick Se plantea también la posibilidad de controlar el simulador de manera interactiva mediante un joystick. Hay algunos artículos al respecto, como en «Virtual excavator simulator featuring HILS and haptic joysticks» [14], de Yun-Joo Nam y Myeong-Kwan Park. En este artículo, se explica cómo introducir el accionamiento del joystick en un software de simulación. El entorno de programación es C++ , y se podría adaptar a lenguaje MatLab, o al revés, adaptar todo el programa al nuevo lenguaje, para observar el rendimiento matemático de las dos maneras. 8.2.7 Realidad virtual Para la visualización de resultados, resulta de interés mostrar el movimiento de una manera más realista. En orden de hacer esto posible, se podrían introducir los sólidos modelados en la interfaz de MatLab de realidad virtual VRML (Virtual Reality Modelling Language), donde se podría simular el entorno incluso, y conseguir tener una mejor idea del movimiento obtenido en las diferentes simulaciones realizadas. Conclusiones 92 93 REFERENCIAS [1] A. López Guinea, Proyecto de diseño y cálculo del brazo articulado de una pala excavadora, Escuela de Ingeniería de Bilbao: Departamento de Ingeniería Mecánica, 2017. [2] Caterpillar, «Excavadoras hidráulicas 33D L 336D LN,» Peoria, Illinois, 2008. [3] J. L. Escalona Franco, «Dinámica de Sistemas Multicuerpo,» Apuntes de la asignatura Cinemática y Dinámica de Máquinas, Universidad de Sevilla, 2017. [4] D. García Vallejo, «Dinámica de sistemas multicuerpo rígido-flexibles en coordenadas absolutas,» (Tesis doctoral) Escuela Técnica superior de Ingeniería, Sevilla, 2006. [5] M. D. Ardema, Newton-Euler Dynamics, Santa Clara, California: Springer, 2005. [6] A. A. Shabana, Computational dynamics, New York: John Wiley and Sons, 1994. [7] J. Baumgarte, «Stabilization of constraints and integrals of motion in dynamical systems,» Computer Methods in Applied Mechanics and Engineering, nº 1, pp. 1-16, 1972. [8] A. Mikkola, Modelling of Hydraulics, Lappeenranta University of Technology. [9] G. Sun y M. Kleeberger, «Dynamic responses of hydraulic mobile crane with consideration of the drive system,» Mechanism and Machine Theory, nº 38, p. 1489–1508, 2003. [10] Parker, «Hydraulic valves industrial standard Catalogue HY11-3500,» de Catalogue HY11-3500, Reino Unido, 2015. [11] P. H. Chang y S.-J. Lee, «A straight-line motion tracking control of hydraulic excavator system,» Mechatronics, nº 12, pp. 119-138, 2002. [12] D. García Vallejo, «Esquema de co-simulación del proceso de sinterización,» Escuela Técnica Superior de Ingeniería, Sevilla. [13] R. A. Ortega Aguilera, «Formulación conservativa de restricciones en sistemas multicuerpo,» (Trabajo de investigación) Universidad Politécnica de Madrid, Madrid, 2006. [14] Y.-J. Nam y M.-K. Park, «Virtual excavator simulator featuring HILS and haptic joysticks,» Journal of Mechanical Science and Technology, vol. 29, nº 1, pp. 397-407, 2015. [15] D. Dopico, A. Luaces, M. González y J. Cuadrado, «Dealing with multiple contacts in a human-in-the- loop,» Multybody Syst Dyn, nº 25, pp. 167-183, 2011.