Repositorio Institucional de Documentos

Proyecto Fin de Carrera Diseño y Control de un Manipulador Flexible Accionado por un Cilindro Hidráulico Tampere University of Technology Proyecto Fin de Carrera Título Diseño y Control de un Manipulador Flexible Accionado por un Cilindro Hidráulico Autor Álvaro Martí Martínez Director Jouni Mattila Ponente Javier Civera Universidad de Zaragoza Tampere University of Technology 2013 1 Proyecto Fin de Carrera Diseño y Control de un Manipulador Flexible Accionado por un Cilindro Hidráulico 2 3 Diseño y Control de un Manipulador Flexible Accionado por un Cilindro Hidráulico RESUMEN El proyecto fin de carrera que se expone a continuación consiste en comparar dos formas de estimación de la posición de un manipulador flexible. Éstos manipuladores flexibles son más ligeros, por lo tanto se pueden desplazar con mayor rapidez y además la carga a desplazar puede ser mayor. El estudio realizado implica el diseño y el modelado del manipulador flexible y del sistema hidráulico actuador. El diseño del manipulador se realiza en MATLAB mediante las ecuaciones de EulerBernoulli para la frecuencia deseada. El modelado del manipulador se desarrolla en el toolbox de MATLAB SimMechanics junto con el sistema hidráulico y todos los sensores necesarios. El sistema hidráulico se compone de una servo válvula o válvula progresiva y un cilindro hidráulico. En el proyecto se desarrolla el modelado de éste servo mecanismo, junto con un controlador y un sistema de lazo cerrado para el control de posición y/o velocidad. La estimación de la posición se realiza mediante dos tipos de sensores distintos, utilizando galgas extensométricas o acelerómetros. Se establecen diferentes puntos de sensado a lo largo del manipulador y con éstos valores se calculan los coeficientes de una función hiperbólica junto con las condiciones de contorno, dicha función hiperbólica aproxima la forma del manipulador. El proyecto también comprende la construcción del prototipo y la comparación de los resultados con los obtenidos en la simulación. Para ello el hardware necesario es una tarjeta de entradas y salidas para el ordenador, encoders, galgas extonsométricas, amplificadores de señal, acelerómetros, válvula progresiva, cilindro hidráulico, transductores de presión... Se desarrolla una interfaz de comunicación mediante el software dSPACE, el cual en base a un archivo de MATLAB Simulink te permite la interacción y visualización de todos los parámetros del programa y así actúa sobre el prototipo. 4 5 Tabla de Contenidos 1 Introducción .............................................................................................................. 9 2 Modelado ................................................................................................................. 12 2.1 Modelo de Barra Flexible ................................................................................ 12 2.1.1 Teoría de Euler-Bernoulli ......................................................................... 12 2.1.2 Diseño de la Barra .................................................................................... 15 2.1.3 Modelo de barra flexible en SimMechanics ............................................. 22 2.2 Modelo del Sistema Servo Hidráulico ............................................................. 24 2.2.1 Dinámica de la Servo Válvula .................................................................. 24 2.2.2 Dinámica del Cilindro Hidráulico ............................................................ 27 2.2.3 Modelo Lineal .......................................................................................... 28 2.2.4 Sistema Hidráulico en SimMechanics ...................................................... 33 2.2.5 Brazo Hidráulico....................................................................................... 34 2.2.6 Frecuencia ................................................................................................. 35 3 Control ..................................................................................................................... 37 3.1 Control del Sistema Servo Hidráulico ............................................................. 37 3.1.1 Servo Válvula ........................................................................................... 37 3.1.2 Cilindro Hidráulico ................................................................................... 38 3.2 Controlador ...................................................................................................... 40 3.3 Estimación de la Forma del Manipulador ........................................................ 43 3.3.1 Galgas Extensométricas ............................................................................ 44 3.3.2 Acelerómetros ........................................................................................... 51 4 Trabajo Experimental .............................................................................................. 56 4.1 Configuración del Hardware ............................................................................ 56 4.1.1 Configuración de las Galgas Extensiométricas ........................................ 59 4.2 Configuración del Software ............................................................................. 61 4.3 Comparación entre Simulink y las Mediciones en el Laboratorio ................... 63 4.3.1 Resultados de la Simulación ..................................................................... 63 4.3.2 Resultados en el Prototipo Experimental ................................................. 64 4.3.3 Comparación Simulación y Prototipo....................................................... 65 4.4 Resultados del Sistema Controlado ................................................................. 66 4.4.1 Inclinación sin Controlador ...................................................................... 66 6 4.4.2 Inclinación con Controlador ..................................................................... 67 5 Conclusiones ........................................................................................................... 70 6 Referencias .............................................................................................................. 72 7 Anexos ..................................................................................................................... 76 7.1 Anexo I: Parametric Study ............................................................................... 76 7.1.1 MATLAB Program .................................................................................. 76 7.1.2 Results ...................................................................................................... 77 7.2 Anexo II: Test Bench, Previous Works ........................................................... 81 7.2.1 Cylinder test bench ................................................................................... 81 7.2.2 Joint construction ...................................................................................... 85 7.2.3 Strength calculations ................................................................................ 85 7.2.4 Joint bearings ............................................................................................ 90 7.3 Anexo III: SimMechanics Model .................................................................... 98 7.4 Anexo IV: Hollow Section Profile ................................................................. 103 7.5 Anexo V: Servo Valve Data Sheet ................................................................ 107 7.6 Anexo VI: MATLAB Files ............................................................................ 110 7.6.1 System Data ............................................................................................ 110 7.6.2 Lumped Method Calculation .................................................................. 111 7.6.3 Hydraulic System Data ........................................................................... 112 7.6.4 Estimation functions ............................................................................... 114 7.7 Anexo VII: Test Bench Drawing ................................................................... 116 7.8 Anexo VIII: Finite Element Analysis with MATLAB .................................. 117 7 Lista de símbolos y abreviaturas  Momento de inercia de la sección del perfil  Módulo de Young () Desplazamiento en dirección  a la distancia  desde el punto fijo  Frecuencia Natural   Frecuencia Natural de la solución n   Frecuencia Natural de la servo válvula   Frecuencia Natural del sistema global ρ Densidad  Densidad de masa () Área de la sección del perfil  Masa del extremo del perfil  Longitud de la barra  Rigidez  Gravedad  Inclinación de la barra 1 Ángulo del punto de anclaje (Cilindro-barra) 2 Ángulo del punto de anclaje (Cilindro-bancada)   Diámetro del carrete de la servo válvula   Posición del carrete   Caudal de carga   Presión de alimentación   Presión de depósito   Presión de cámara A   Presión de cámara B   Ganancia de flujo 8  Coeficiente de presión de flujo   Ganancia de la válvula ! Módulo de compresibilidad efectivo "  Flujo cámara A "  Flujo cámara B #  Amortiguamiento de la servo válvula  $ Área del pistón  $ Posición del vástago del cilindro % & Coeficiente de pérdidas '  Fricción viscosa 1 Introducció n En el campo de la robótica se conoce como bajo el control humano manipula objetos sin contacto directo y que además presenta un comportamiento elástico an cargas o movimientos ( Figura 1.1 La i nvestigación sobre el campo de importante desde hace más de una década a causa de las nuevas aplicaciones de la r obótica en diferentes áreas, tales como las tareas industriales. La tendencia es materiales ligeros en las estructuras de los manipuladores con el fin de mejorar el rendimiento de los robots industriales actuales, que son grandes y pesados Figura 1.1 En aplicaciones aeroespaciales la grandes estructuras de control como grúas lleva a cabo con instrumentos flexibles manipuladores automáticos precisos. Los mani puladores flexibles presentan ventajas importantes en comparación con los rígidos, tales como la reducción de consumo de energía, se necesitan accionadores más n En el campo de la robótica se conoce como manipulador flexible a un dispositivo que bajo el control humano manipula objetos sin contacto directo y que además presenta un comportamiento elástico an te diferentes tipos de situaciones como la aplicación de Figura 1.1 ). nvestigación sobre el campo de los manipuladores flexibles , se ha vuelto más importante desde hace más de una década a causa de las nuevas aplicaciones de la obótica en diferentes áreas, tales como las tareas industriales. La tendencia es materiales ligeros en las estructuras de los manipuladores con el fin de mejorar el rendimiento de los robots industriales actuales, que son grandes y pesados Figura 1.1 Esquema de un Manipulador Flexible. En aplicaciones aeroespaciales la ligereza de los materiales también es importante, grandes estructuras de control como grúas o la cirugía mínimamente invasiva lleva a cabo con instrumentos flexibles y delgado s, en donde son necesarias manipuladores automáticos precisos. puladores flexibles presentan ventajas importantes en comparación con los rígidos, tales como la reducción de consumo de energía, se necesitan accionadores más 9 a un dispositivo que bajo el control humano manipula objetos sin contacto directo y que además presenta un te diferentes tipos de situaciones como la aplicación de , se ha vuelto más importante desde hace más de una década a causa de las nuevas aplicaciones de la obótica en diferentes áreas, tales como las tareas industriales. La tendencia es a utilizar materiales ligeros en las estructuras de los manipuladores con el fin de mejorar el rendimiento de los robots industriales actuales, que son grandes y pesados . ligereza de los materiales también es importante, en la cirugía mínimamente invasiva que se s, en donde son necesarias puladores flexibles presentan ventajas importantes en comparación con los rígidos, tales como la reducción de consumo de energía, se necesitan accionadores más 16 El primer paso es diseñar un apropiado banco de pruebas para obtener los valores correctos. El banco de pruebas está compuesto por un cuerpo flexible con una masa en el extremo y unido a la bancada por una articulación de rotación. Éste prototipo es accionado por un circuito servo hidráulico donde el cilindro está acoplado entre el cuerpo flexible y la bancada. 2.1.2.1 Dimensionado Para comenzar, la frecuencia es una propiedad importante a considerar, teniendo en cuenta las mediciones y la comparación de los datos a obtener. En el diseño de la barra, la frecuencia natural influiría en el dimensionado, pero la frecuencia global será calculada incluyendo el conjunto servo hidráulico más adelante. Por otra parte, además de la frecuencia, la flecha que la barra sufre sobre su extremo libre debe ser considerable con el fin de comparar las estimaciones de la deformación. Por consiguiente, un estudio paramétrico ha sido desarrollado para dimensionar la barra flexible. El objetivo de éste análisis es concluir cuáles son las medidas a seleccionar para el perfil hueco rectangular. Las propiedades incluidas en el estudio son la deformación sobre el extremo libre, la frecuencia natural, la tensión máxima y el pandeo. La fórmula de la flecha sobre el extremo de la barra. VW=' X= 2 3 Obtenida del desarrollo de la ecuación de Euler-Bernoulli. La siguiente fórmula es la que resuelve la tensión máxima, producida sobre el punto de máximo momento flector que es en el extremo de unión a la bancada. Y=Z 2 Donde es el momento producido sobre la barra en dicho punto. YZ/2 es la distancia entre la parte superior de la sección transversal y la linea neutral. La fuerza de pandeo crítica, [wikipedia, Buckling 2013], para la barra en voladizo es: 17 '=\ )  ] ) ≫,S9=]=2 (3.1.2.1.3) La última fórmula es una aproximación para obtener la frecuencia natural de sistemas con masas concentradas y distribuidas al mismo tiempo. _==` a bc?.)2∗d X= 2∗e∗f g 8== h )i (3.1.2.1.4) Se proponen los límites para el estudio paramétrico a través de las propiedades mecánicas de el material como la tensión máxima del acero, la frecuencia máxima o la fuerza de pandeo crítica. • La máxima frecuencia propuesta es 5 Hz. • La máxima tensión del acero es 355 MPa y con un coeficiente de seguridad de 2 queda en 177.5 MPa. • La fuerza de pandeo crítica es igual a la fuerza producida por el peso de la masa del extremo del manipulador. Una vez que los límites han sido definidos, el estudio paramétrico puede ser resuelto. El análisis es ejecutado como una iteración en la que en cada paso se varían los parámetros del manipulador y el peso de la masa del extremo. El programa implementado en MATLAB se muestra en el Anexo 8.1 junto con las gráficas de resultados. Los parámetros a iterar son la altura del perfil, anchura, espesor, longitud del manipulador y peso de la masa acoplada en el extremo. Finalmente, los valores seleccionados de la barra se muestran en la Tabla 3.1, junto con el intervalo de valores admisible para no sobrepasar los límites. Tabla 3.2 Valores seleccionados tras el análisis y variación permitida. Parámetro Valor Variación permitida por los límites Altura (H) 0.08 m 0.0401-0.1 Anchura (D) 0.06 m 0.0101-0.1 Espesor (e) 0.004 m 0.0011-0.01 Longitud (L) 2.5 m >2.2 Masa sobre el extremo (M) 50 Kg 32-165 18 El perfil seleccionado para el manipulador (Ruukki double grade S355J2H) es una sección hueca estructural que cumple los requerimientos del estándar EN 10219. Ha sido seleccionado debido a sus aplicaciones en estructuras en maquinaria, vehículos de transporte y equipos de elevación. Todas las propiedades mecánicas se muestran en el Anexo 8.4. 2.1.2.2 Banco de pruebas Para el desarrollo del prototipo es necesario diseñar un banco de pruebas que permita el correcto movimiento del manipulador, y que favorezca la actuación del cilindro hidráulico. Anteriormente, fue desarrollado un banco de pruebas previo en el departamento de hidráulica de la Universidad de Tampere para estudiar la deformación Anexo 8.2. El problema fue que la barra no era lo suficientemente flexible, por lo tanto, se tiene que desarrollar un nuevo banco de pruebas para el manipulador, aunque basado en el anterior diseño. Algunas partes del anterior banco de pruebas son aprovechadas para el nuevo, como la columna principal y cilindro hidráulico, pero la distancia entre la unión del cilindro y la de la barra, 1, tiene que ser recalculada. La longitud 1 será superior a la anterior debido a las fuerzas soportadas sobre el cilindro. Un análisis de fuerzas es expuesto a continuación con el fin de obtener el valor adecuado para la longitud requerida (1). 19 Figura 3.1 Esquema del banco de pruebas. El ángulo del manipulador , ángulo de la pendiente del manipulador con referencia a la horizontal, es definido por la geometría del banco de pruebas y el recorrido del cilindro. \ 2+=1+2+ _Sk=` K ) + ) ) −2∗ K ∗ ) ∗cos La fuerza aplicada sobre el cilindro es obtenida por la momento de carga y trigonometría. La Figura 3.2 muestra la fuerza ejercida sobre el cilindro hidráulica en la posición de máxima carga. El cilindro tiene las siguientes propiedades: Dc = 32 mm, Carrera = 300 mm, Fuerza Max = 16889 N (210 bar). 20 Figura 3.2 Fuerza a soportar por el cilindro en función de L1 con inclinación nula. De acuerdo con las simulaciones, la longitud 1debe encontrarse entre 250 y 300 mm, si es menor, el sistema será demasiado lento para que se produzca una oscilación en la barra flexible, y si es mayor, el recorrido del conjunto se reduce debido a la carrera del cilindro. La longitud elegida para el banco de pruebas es 280 mm. Entonces a1 es 17.8 grados y a2 es 42.5 grados. La altura entre la parte superior del cilindro y la línea neutral de la barra se establece en 90 mm con el fin de mantener suficiente espacio para la unión cilindro-manipulador. La fuerza mínima necesaria para desplazar el conjunto es 6812.5 N. La fuerza que soporta el cilindro a lo largo de la rotación del manipulador se muestra en la Figura 3.3. La Figura 3.4 expone la posición necesaria de la carrera del cilindro para conseguir la inclinación deseada. Si la carrera del cilindro es 300 mm el ángulo del manipulador está limitado entre 0 y 63 grados 0.2 0.25 0.3 0.35 0.4 0.45 0.5 0.55 0.6 0.65 0 0.5 1 1.5 2 2.5 x 10 4 Applied Force in Maximum load position (Theta=0) Lenght 1 [m] Force [N] 21 Figura 3.3 Fuerza mínima para cada posición con L1=280 mm. Figura 3.4 Posición del cilindro en función de la inclinación. 0 10 20 30 40 50 60 70 80 90 0 1000 2000 3000 4000 5000 6000 7000 Needed force along the movement Theta [deg] Force [N] 0 10 20 30 40 50 60 70 80 90 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 Cylinder extension along the movement Theta [deg] Cyl Position [m] 22 2.1.3 Modelo de barra flexible en SimMechanics 2.1.3.1 Método de los cuerpos agrupados Para la mayoría de los propósitos de ingeniería, un cuerpo flexible real es un medio continuo. En el método de los cuerpos agrupados se aproxima un cuerpo flexible como un conjunto de cuerpos rígidos, junto con muelles y amortiguadores, y en SimMechanics, puede ser implementado por una cadena alterna de cuerpos rígidos y articulaciones. Los muelles y amortiguadores actúan ya sea en los cuerpos o las articulaciones. Los coeficientes de muelle y de amortiguación son funciones de las propiedades del material y de la geometría del manipulador flexible de que se trate. Los métodos cuerpos agrupados son los más adecuados en los modelos con geometrías lineales, tales como vigas, en el que cada elemento fundamental está acoplado a otros dos en una cadena simple, [Shye & Richardson, 1987]. En la siguiente figura se muestra un esquema del conjunto de cuerpos discretos. Figura 3.5 Esquema del método de cuerpos agrupados. El método de cuerpos agrupados sigue los siguientes pasos, [Chudnovsky, Mukherjee, Wendlandt & Kennedy, 2006]: 1. Dividir la viga en elementos discretos. 2. Determinar los grados de libertad de los elementos y modelar apropiadamente uniones primitivas. 3. Aplicar un los coeficientes de muelle y amortiguación a cada junta para encapsular los parámetros. 4. Utilizar la teoría del cuerpo flexible para determinar las constantes elásticas de la geometría, las propiedades del material y las condiciones de contorno. 5. Aplicar la amortiguación necesaria para cada grado de libertad. 6. Construir la viga mediante la soldadura de una cadena de conjuntos elásticos. 23 A continuación se describe uno de dichos conjuntos. Figura 3.6 Elemento de la cadena de cuerpos agrupados. También se agregan los parámetros de muelle y amortiguador de la articulación. La Figura 3.7 muestra la unión del bloque de amortiguación a la unión y la Figura 3.8 muestra en interior del bloque. Figura 3.7 Diagrama Simulink de un elemento. Figura 3.8 Aplicación de la ecuación del momento sobre la unión. Ecuación normalizada del momento: 24 l=+2#9m=+9 ) ==9n=494n:=9 Dónde 9 2 =X/ e  es el momento de inercia. El coeficiente de amortiguamiento 2ξω o es un valor cuasi-empírico que da cuenta de la energía perdida a efectos viscoelásticos. En SimMechanics, el ángulo relativo = y la velocidad de la unión= m se miden en la articulación con un bloque conjunto del sensor. Entonces Xse multiplica por el ángulo, y el coeficiente de amortiguación del material, 2ξωo, se multiplica por la velocidad angular. Se añaden los dos momentos resultantes y se aplican de nuevo a la unión utilizando un actuador de articulación. En la siguiente figura se muestra el resultado de la barra flexible en la simulación en MATLAB junto con el actuador que se modelará a continuación. Figura 3.9 Simulación del manipulador flexible. 2.2 Modelo del Sistema Servo Hidráulico 2.2.1 Dinámica de la Servo Válvula La servo válvula más extendida en el mercado es la válvula carrete o corredera que emplea una bobina para el mando del eje. Las válvulas de carrete están clasificadas por el número de vías que el fluido puede seguir. Todas las válvulas requieren vías de alimentación y de retorno, y al menos una vía para la carga, las válvulas suelen ser de tres o cuatro puertos. 25 El centro de la servo válvula es otro tipo de clasificación (Figura 3.10), se define por la anchura del asiento de la válvula en comparación con el orificio de la vía cuando el carrete se encuentra centrado. Figura 3.10 Clasificación de la válvula según su tipo de carrete centrado. Figura 3.11 Caudal de la válvula en función de la posición y tipo de carrete. El área de flujo (Figura 3.11), para un carrete de diámetro   y con desplazamiento,  se calcula como a continuación, (  )=8  \    Donde f pq es la fracción de la circunferencia de carrete que se está abierta en la ranura. Normalmente éste parámetro está definido en el siguiente intervalo: 0.25≤8  ≤0.5. 32  ¯ =J4! •   # ¯ =( u +%)  $ J! •  +'  2J• ! 2.2.3.1 Modelo de Espacio de Estados Las ecuaciones del espacio de estados del sistema hidráulico son las siguientes, E∗›(E)=∗›(E)+!∗°(E) (E)=%∗›(E)+F∗°(E) Donde ›=(m    ) y °=  = ± ² ² ² ² ² ³ −'   $ − $  −! $ •  −! u •  0 ! $ •  0 ! u •  ´ µ µ µ µ µ ¶ ;!= ± ² ² ² ³ 0 ! † •  −! † •  ´ µ µ µ ¶  %=¸1 0 0¹;F=0 Es momento de calcular la rigidez hidráulica del sistema. El módulo de compresibilidad B se define de la siguiente manera, !=−• ? V V• donde V ? es el volumen inicial de fluido, Δp y ΔV son las variaciones de presión y de volumen. !=−•• ? + $,?  $ ‘V'  $ 1  $ V 33 V'=  $ ) ! •• ? + $,?  $ ‘V  » = $ ) ! •• ? + $,?  $ ‘ De manera similar, la segunda constante de muelle, K ½ž , se calcula para la otra cámara,  » = $ ) ! •• ? +(− $,? ) $ ‘ Y la rigidez del cilindro hidráulico se define con la siguiente aproximación.  » = » + » = $ ) ! •• ? + $,?  $ ‘+ $ ) ! •• ? +(− $,? ) $ ‘ Por tanto, la frecuencia natural se calcula como,   =J » . 2.2.4 Sistema Hidráulico en SimMechanics El circuito hidráulico compuesto por una válvula proporcional con 4 vías y 3 posiciones y varios manómetros para medir el cilindro y la presión de suministro desarrollado en SimMechanics se expone ahora (Figura 3.12). La entrada se adapta a la posición de la válvula de carrete. 34 Figura 3.12 Diagrama del circuito hidráulico en Simulink. En el modelo de Simulink completo se puede observar la unión del conjunto hidráulico con el cilindro y a su vez la unión mecánica con el manipulador flexible. En el Anexo 8.3 muestra el esquema del sistema completo con la parte del manipulador flexible, el sistema hidráulico, cilindro y los bloques de adaptación de la referencia así como los bloques de estimación de la posición. 2.2.5 Brazo Hidráulico La rigidez lineal se puede transformar a la correspondiente rigidez a torsión como se expresa a continuación. V¾=Œ » V, V¾=V'= » V= » V, 35 Œ » = )  » Finalmente, la frecuencia natural se traduce transforma así.   =JŒ » ¿,¿=EE∗ ) 4+À;EE∗ )  2.2.6 Frecuencia El diseño de sistemas hidráulicos a menudo se basa en un modelo de carga simple, representado por parámetros individuales. Este tipo de modelo de carga se puede utilizar si la conexión entre el sistema hidráulico y la carga mecánica es rígida. Sin embargo, en muchas aplicaciones el sistema mecánico, al que los elementos hidráulicos de potencia están conectados, es débil en comparación con la rigidez del sistema hidráulico. Tales estructuras mecánicas causan problemas de resonancia, porque esa frecuencia puede ser más baja que la frecuencia natural hidráulico. Si la resonancia estructural domina la respuesta de frecuencia del sistema de servo, es extremadamente importante tener en cuenta este hecho en el diseño del sistema. Las principales razones son la estabilidad de los sistemas y el ancho de banda que está limitado por la frecuencia natural más baja en el bucle de control. Diagrama de Bode del sistema servo hidráulico en bucle cerrado: Figura 3.13 Diagrama de Bode del sistema Hidráulico. -100 -50 0 50 100 Magnitude (dB) 10 3 10 4 10 5 10 6 -180 -135 -90 -45 0 Phase (deg) Bode Diagram Frequency (rad/s) 36 Frecuencia natural del sistema servo: ω Á =29634rad s,f Á =4716.4Hz Rigidez hidráulica: K ½ =A Ç ) B •V Ÿ? +x Ç,? A Ç ‘+A Ç ) B •V ž? +(L−x Ç,? )A Ç ‘=11176242N Frecuencia natural del brazo hidráulico: Œ » = ) ∗ »   =JŒ » ¿=54.92: E8  =8.74ZË Como se ha explicado anteriormente, la frecuencia natural del sistema hidráulico es más lenta que la frecuencia del sistema servo, por lo tanto, se podría producir un problema de resonancia, pero en este caso es prácticamente imposible debido a que la frecuencia del manipulador es tres órdenes de magnitud menor que la frecuencia del conjunto hidráulico. 37 3 Control El objetivo de ésta sección es desarrollar un controlador simple y adecuado para el sistema servo hidráulico y proporcionar una estimación para la posición del extremo del mediante los valores de los sensores utilizados. En éste proyecto se ha diseñado un sistema de posición servo para alcanzar la inclinación adecuada del manipulador flexible, proporcionando un ángulo de la articulación, sin duda, bien medido y rigidez para comenzar con la estimación de flexión. En ésta sección se proponen dos métodos de estimación, el primero a través de la medición de galgas extensométricas y el segundo mediante acelerómetros. Éstos sensores son acoplados a lo largo del manipulador. 3.1 Control del Sistema Servo Hidráulico Un sistema servo es un sistema de control que mide su propia salida y fuerza a dicha salida a seguir la señal de mando de forma rápida y precisa. De éste modo, el efecto de anomalias en el propio dispositivo de control y en la carga se pueden minimizar, al igual que sucede con las perturbaciones externas. Un sistema servo puede ser diseñado para controlar casi caulquier valor físico medible, como por ejemplo, movimiento, fuerza, presión, temperatura, tensión o corriente. Los parámetros del sistema servo que se va a utilizar en el estudio y en la experimentación se muestran en la siguiente Tabla 4.1: Tabla 4.1 Parámetros del sistema hidráulico. Carrera (L) 0.3 m Diámetro del pistón (D) 0.032 m Diámetro del vástago (Dr) 0.025 m Módulo de compresibilidad (B) 1e6 Masa del manipulador (M) 50 Kg Coeficiente pérdidas (Cp) 0 Fricción viscosa (Fv) 20000 N Flujo Nominal (Qn) 50 l/min Posición máxima válvula 1 Presión Nominal (Pn) 5e5 Pa Parámetro de la válvula (Tau) 30 e-3 Presión de Alimentación (Ps) 210e5 Pa 3.1.1 Servo Válvula La dinámica de la servo válvula se desarrollará de acuerdo con la hoja de características que proporciona el fabricante (Válvula BOSCH 4WRPH 6 Anexo 8.5), el tiempo de 38 respuesta es menor a 10 ms. Y la función de transferencia puede ser aproximada por un sistema de segundo orden con los siguiente parámetros, tal y como se ha explicado en la sección de modelado. #  =0.8−1,À — =0.01 Œ  (E)=  E )  ) +2#  E   +1 De la teoría de los sistemas de segundo orden, se obtiene la fórmula aproximada para el tiempo de respuesta: #  ∗  ∗À — ≅4 por tanto, la frecuencia natural es:   =400:/E  3.1.2 Cilindro Hidráulico La función de transferencia del sistema hidráulico se compone de las funciones de la válvula y del cilindro: Œ(E)=Œ  (E)∗Œ ÍÎ (E)=Œ  (E)∗ ¯ E¦E )  ¯) +2# ¯ E  ¯ +1© =ŒÏ(E)∗2.627n004 E 2 +4.865E ) +3913E  ¯ =J4! •   # ¯ =( u +%)  $ J! •  +'  2J• ! En la siguiente Tabla 4.2 se muestra la localización de los polos del sistema: 39 Tabla 4.2 Localización de los polos del sistema en lazo abierto. LOCALIZACIÓN DE LOS POLOS DEL SISTEMA, Œ ( E ) Real Imaginario 0 -400 -400 -2.43 -2.43 0 0 0 62.5 -62.5 Los polos correspondientes a la válvula se encuentran localizados en -400 sobre el eje real, así que la dinámica de la válvula es mucho más rápida que la del cilindro, como se muestra en el gráfico del lugar de las raíces (Figura 4.1). Por lo tanto, no tiene gran influencia sobre el comportamiento del sistema. Figura 4.1 Lugar de las raíces de la función de transferencia del sistema hidráulico. -400 -300 -200 -100 0 100 -300 -200 -100 0 100 200 300 Root Locus Real Axis (seconds -1 ) Imaginary Axis (seconds -1 ) 40 3.2 Controlador Una aproximación estándar para controlar los sistemas servo hidráulicos es aplicar métodos de control de realimentación lineal. El objetivo del controlador de posición es imitar el comportamiento de la trayectoria, aumentando el amortiguamiento y evitando la inestabilidad. Ésta realimentación influye sobre la localización del par de polos conjugados. Hay muchas opciones de realizar control de posición, e incluso tan simples como el Pcontrol. De acuerdo con [Jelai & Kroll 2004], si se utiliza un controlador PT (primer orden con desplazamiento), se alcanzarán respuestas más rápidas a un escalón en el sistema realimentado. El ratio de amortiguamiento del sistema sin controlador es de 0.039. Œ (E)= 1+EÀ Los parámetros y À son sintonizados para lograr un ratio de amortiguamiento superior a 0.7 y evitar la inestabilidad. À =0.33;  <0.2; Se halla la función de transferencia en lazo cerrado para la entrada › ,Е“ , incluyendo el controlador. Œ Î (E)=› $ (E) ›  (E)=Œ (E)Œ(E) 1+Œ (E)Œ(E) 41 Finalmente, la localización de los polos del sistema queda de la siguiente manera (Tabla 4.3): Tabla 4.3 Localización de los polos del sistema en lazo cerrado. LOCALIZACIÓN DE LOS POLOS DEL SISTEMA, Œ Î ( E ) Real Imaginario -1.5142 -1.5142 -2.4334 -2.4334 1.3344 -1.3344 62.4703 -62.4703 donde, como se muestra en la Figura 4.2, los polos complejos en -2.4334 disminuyen el amortiguamiento global, por lo tanto  podría disminuirse más. Figura 4.2 Lugar de las raíces incluyendo el controlador. Localización final de los polos del sistema mostrada en la Figura 4.3. -6 -4 -2 0 2 4 6 -60 -40 -20 0 20 40 60 0.0160.0340.0540.08 0.115 0.16 0.26 0.45 10 20 30 40 50 60 70 10 20 30 40 50 60 70 0.0160.0340.0540.08 0.115 0.16 0.26 0.45 Root Locus Real Axis (seconds -1 ) Imaginary Axis (seconds -1 ) 48 1. Ô( K )=− ’∗ÍÛÛ(…Ü) ) 2. Ô( ) )=− ’∗ÍÛÛ(…Ý) ) 3. Ô( 2 )=− ’∗ÍÛÛ(…g) ) Los tres puntos de medición han sido seleccionados en el mismo punto donde la simulación de cuerpos agrupados tienen la unión entre ellos. Éstas deformaciones son medidas en tiempo real por un Medio Puente de Wheatstone implementado por MatLab Simulink. Figura 4.8 Diagrama en Simulink de un medio puente de Wheatstone. Ganancia a la salida del puente: Œ;== −2 •  ∗Œ' donde •  es la tensión de alimentación, 24•, y Œ' es el Factor de Galga, 2.11. La medición de la deformación no puede ser posicionada en el punto de apoyo simple debido a las condiciones de contorno en dicho punto, momento igual a 0, aparecería un sistema incompatible. 49 A continuación es necesario el algoritmo que resuelve la ecuación de forma. Se realiza un cálculo matricial para hallar los nueve coeficientes con éstas nueve ecuaciones evaluadas cada una en la distancia correspondiente. ∗%9n88=% (=0)=E;=ℎ()+!E;=ℎ(2)+%∗E;=ℎ(3)+FE;=ℎ(4)+S9Eℎ() +'S9Eℎ(2)+ŒS9Eℎ(3)+ZS9Eℎ(4)+=0  ØØ (=0)=S9Eℎ()+E;=ℎ()+4'S9Eℎ(2)+9ŒS9Eℎ(3) +16ZS9Eℎ(4)+4!E;=ℎ(2)+9%E;=ℎ(3) +16FE;=ℎ(4)=0 (=1)=E;=ℎ()+!E;=ℎ(2)+%E;=ℎ(3)+FE;=ℎ(4)+S9Eℎ() +'S9Eℎ(2)+ŒS9Eℎ(3)+ZS9Eℎ(4)+=0  ØØ (=1)=S9Eℎ()+E;=ℎ()+4'S9Eℎ(2)+9ŒS9Eℎ(3) +16ZS9Eℎ(4)+4!E;=ℎ(2)+9%E;=ℎ(3) +16FE;=ℎ(4)=0  ØØ (=)=S9Eℎ()+E;=ℎ()+4'S9Eℎ(2)+9ŒS9Eℎ(3) +16ZS9Eℎ(4)+4!E;=ℎ(2)+9%E;=ℎ(3) +16FE;=ℎ(4)=0  ØØØ (=)=S9Eℎ()+E;=ℎ()+8!S9Eℎ(2)+27%S9Eℎ(3) +64FS9Eℎ(4)+8'E;=ℎ(2)+27ŒE;=ℎ(3) +64ZE;=ℎ(4)=0 Ô(=1)=−(S9Eℎ())/25−(E;=ℎ())/25−(4'S9Eℎ(2))/25 −(9ŒS9Eℎ(3))/25−(16ZS9Eℎ(4))/25 −(4!E;=ℎ(2))/25−(9%E;=ℎ(3))/25 −(16FE;=ℎ(4))/25 Ô(=2)=−(S9Eℎ())/25−(E;=ℎ())/25−(4'S9Eℎ(2))/25 −(9ŒS9Eℎ(3))/25−(16ZS9Eℎ(4))/25 −(4!E;=ℎ(2))/25−(9%E;=ℎ(3))/25 −(16FE;=ℎ(4))/25 Ô(=3)=−(S9Eℎ())/25−(E;=ℎ())/25−(4'S9Eℎ(2))/25 −(9ŒS9Eℎ(3))/25−(16ZS9Eℎ(4))/25 −(4!E;=ℎ(2))/25−(9%E;=ℎ(3))/25 −(16FE;=ℎ(4))/25 Donde 50 %9n88=(!%F'ŒZ) y %=(000000Ô( K )Ô( ) )Ô( 2 )) Después del cálculo de los coeficientes, la deformación en el extremo del manipulador se puede obtener como: ()=E;=ℎ()+!E;=ℎ(2)+%E;=ℎ(3)+FE;=ℎ(4)+S9Eℎ() +'S9Eℎ(2)+ŒS9Eℎ(3)+ZS9Eℎ(4)+ Los algoritmos de MATLAB están expuesto en el Anexo 8.6. A continuación se expone un ejemplo de la deformación a lo largo de la posición , aplicando unas deformaciones arbitrarias sobre las galgas. Strain1 10e-5 Strain2 5e-5 Strain3 2e-5 Figura 4.9 Resultado de la función de forma del manipulador tras una deformación arbitraria. La implementación en Simulink se muestra en la Figura 4.10, donde se realiza la multiplicación matricial para la obtención de los coeficiente de la ecuación de forma. En cada paso del programa éstos coeficientes son recalculados, así que la ecuación es reajustada a la forma del manipulador en cada iteración. 0 0.5 1 1.5 2 2.5 -3.5 -3 -2.5 -2 -1.5 -1 -0.5 0 0.5 x 10 -3 51 Figura 4.10 Esquema del cálculo de la ecuación de forma en Simulink. 3.3.2 Acelerómetros Un acelerómetro es un dispositivo que convierte la aceleración en una señal eléctrica. Tanto la aceleración estática y dinámica se pueden medir usando un acelerómetro, donde la aceleración dinámica es la aceleración debida a cualquier la fuerza a excepción de la fuerza gravitacional aplicada sobre un cuerpo rígido, y la aceleración estática (o aceleración de la gravedad) es debida a la fuerza gravitatoria. Tres de las técnicas más importantes para medir la aceleración se explican brevemente a continuación, [Naghshineh, Ameri, Zereshki, Krishnan & Abdoli]. • Acelerómetro piezoeléctrico. Este tipo de acelerómetros hace uso del efecto piezoeléctrico. A medida que el elemento piezoeléctrico se comprime debido a la fuerza causada por la aceleración, se genera una señal eléctrica proporcional. Los acelerómetros piezoeléctricos no son adecuados para la medición de las condiciones de aceleración cero (DC), pero son muy adecuadas para vibraciones de alta frecuencia. • Acelerómetro capacitivo. Este tipo de acelerómetro es similar a la del acelerómetro piezoeléctrico excepto que utiliza efecto capacitivo. Dado que este tipo de acelerómetros puede ser micro mecanizado en el silicio, puede ser utilizado en circuitos integrados. • Acelerómetro térmico. Dentro de un acelerómetro térmico hay un calentador para calentar una pequeña burbuja de aire en el interior del circuito integrado. Se cambia la posición de la burbuja de aire caliente a medida que se aplica una fuerza en el acelerómetro. El movimiento de la burbuja calentada se mide por los sensores de temperatura y luego se convierte en una señal eléctrica. 52 Medición de la inclinación mediante dos ejes La utilización de más de un eje para calcular la inclinación produce una solución mas exacta que mediante un solo eje. Considerando trigonometría básica, la inclinación se obtiene combinando … y  Í (Figura 4.11). Figura 4.11 Esquema de cálculo de la inclinación mediante dos ejes.  … =sin,  Í =cos =arctan …  Í  Es importante darse cuenta que la combinación de las aceleraciones es siempre 1g: =` …) + Í) =1 Por eso, ésta forma de medición de la pendiente es válida únicamente para cálculos estáticos, pero el manipulador necesita calcular la posición angular mientras se encuentra en movimiento. Debemos considerar un algoritmo diferente. Adicionalmente, para pasar de la aceleración a posición angular, la señal tiene que ser integrada dos veces. Ésto implica que los errores se acumulan. A me nos que la señal de el acelerómetro se muestree infinitamente rápido, algunas mediciones de la aceleración se perderán. Por otra parte, si hay un ruido o una vibración en la señal, puede afectar a la precisión de los resultados. Todos éstos errores se acumulan, así que cuantas menos integraciones tenga el algoritmo, más precisa será la solución. 53 El problema es que entre el modelo de simulación y las mediciones reales existe una gran diferencia, la gravedad tiene que ser considerada. En el caso del proyecto la gravedad tiene que ser substraída de las mediciones reales. Ésta corrección implica el conocimiento de la inclinación, para modificar los ejes necesarios. Medición de la inclinación mediante la aceleración centrípeta La aceleración centrípeta es una medida proporcional a la fuerza centrípeta. Durante un movimiento circular, la dicha aceleración es la componente radial de la aceleración del objeto. Es la tasa de cambio en la velocidad del objeto en movimiento. Es posible relacionar la velocidad angular con la aceleración centrípeta mediante la siguiente fórmula:  = ) :,=` : La aceleración centrípeta está siempre dirigida hacia el centro del círculo, que forma la trayectoria de giro del objeto. Por consiguiente, si la trayectoria del acelerómetro es constante, como se puede asumir en el modelo del proyecto, la velocidad angular obtenida no cambia de dirección de rotación. Obtener la velocidad angular mediante la aceleración centrípeta evita un proceso de integración. Por tanto, los errores acumulados se reducen. La Figura 4.12 muestra la implementación en Simulink para éste cálculo, Figura 4.12 Diagrama del cálculo de la inclinación mediante la aceleración centrípeta en Simulink. 54 donde el camino superior obtiene la velocidad angular, y el camino inferior le da el sentido de giro a dicha velocidad. La velocidad angular se calcula mediante la aceleración paralela a la dirección del manipulador y la dirección mediante la integración de la dirección perpendicular. Figura 4.13 Dirección de la fuerza centrípeta en función de la trayectoria y la velocidad. Además, las dos señales se filtran mediante un filtro de Kalman que retrasa el resultado. Entonces, integrando el resultado obtenemos la inclinación del manipulador en cada punto de medición. Ésta inclinación en cada punto se utiliza en el algoritmo de estimación después de restarle la inclinación al principio de la barra. De este modo, la pendiente en cada punto es relativa a la propia barra.  …7 = 7 − ×Îáâ–Î Si conocemos la pendiente en cada punto de medición se puede estimar la forma de la barra y con ello la posición final. En éste caso se utiliza una ecuación polinómica de 7 coeficientes para estimar la forma del manipulador. Por tanto se necesitan cuatro condiciones de contorno junto con tres puntos de medición. ()= ã +! ä +% . +F 2 + ) +'+Œ Las condiciones de contorno de una barra en voladizo: 1. La deformación al principio del manipulador es 0, y(0) = 0. 2. La derivada de la deformación al principio del manipulador es 0, y '(0) = 0. 3. Momento en el extremo del manipulador es 0, y '' (L) = 0. 55 4. Esfuerzo cortante en el extremo es 0, y '''(L) = 0. Las tres ecuaciones de medición de la aceleración: 1. tan …K = Ø ( K ) 2. tan …) = Ø ( ) ) 3. tan …2 = Ø ( 2 ) Los tres puntos de medición han sido colocados en los mismos puntos que en el método de las galgas extensométricas. El algoritmo es similar al anterior con la ecuación hiperbólica, cada tiempo de muestreo se calculan los coeficientes mediante una operación matricial que contiene las ecuaciones que se presentan a continuación. ∗%9n88=% (0)=0 ã +!0 ä +%0 . +F0 2 +0 ) +'0+Œ=0 ′(0)=60 ä +5!0 . +4%0 2 +3F0 ) +20+'=0  ØØ ()=30 . +20! 2 +12% ) +6F+2=0 ′′′()=120 2 +60! ) +24%+6=0 tan …K =6 Kä +5! K. +4% K2 +3F K) +2 K +' tan …) =6 )ä +5! ). +4% )2 +3F )) +2 ) +' tan …2 =6 2ä +5! 2. +4% 22 +3F 2) +2 2 +' Donde, %9n88=(!%F'Œ) y %=(0000tan …K tan …) tan …2 ) Después del cálculo de los coeficientes, se puede obtener la desviación en el extremo del manipulador: ()= ã +! ä +% . +F 2 + ) +'+Œ 56 4 Trabajo Experimental Las investigaciones experimentales del proyecto se han desarrollado en el Laboratorio del Departamento de Hidráulica y Automática de la Universidad de Tecnología de Tampere. La configuración del sistema de control se muestra esquemáticamente en la Figura 5.1. Consiste en un manipulador flexible, un actuador y un PC. El software utilizado es MATLAB / Simulink construido sobre dSPACE que es la plataforma de las simulaciones en tiempo real. Figura 5.1 Esquema de las conexiones de Hardware sobre el prototipo. Debido a un problema de fabricación de la viga se construyó más corta que la diseñada, la barra utilizada para la prueba es de dos metros de largo en vez de dos y medio. Esto cambia el modelo de comportamiento, se aumenta la frecuencia y la disminución de la deflexión. Por lo tanto, los parámetros de Simulink se manipulan para comportarse como el cuerpo de prueba. La masa de la punta se incrementa con el fin de tener un sistema más lento. 4.1 Configuración del Hardware La configuración de hardware es la interfaz entre el hardware y el software, ahora se muestra la lista de componentes: • Barra de sección hueca estructural, Ruukki double grade S355J2H • Cilindro hidráulico 32mm 300mm 57 • Servo Válvula Hidráulica Bosch Rexroth AG 4WRPH 6 C4B24L–2X/G24Z4 /M • Tarjeta Amplificadora de la Válvula, Bosch PL10 • Galgas Extensométricas KYOWA KFG-5-120-C1 • Amplificadores para las Galgas Extensométricas, HBM Clip AE101 • Transductores de presión, GE PTX 1400 • Encoder Incremental, HEIDENHAIN ROD-426 • Panel de Conexiones CLP1103 • PC Windows XP • Fuente de alimentación 0-24 V DC La entrada del sistema es generada por el equipo que envía la señal a la CLP1103 que es la tarjeta de adquisición de datos. El Panel de conexiones CLP1103 (Figura + + +) contiene los conectores para veinte entradas analógico-digitales, ocho salidas de señal digital a analógica, y varios conectores que pueden ser utilizados para E / S digitales, Esclavo / DSP I / O, interfaz del encoder incremental, interfaz CAN e interfaces de serie. Figura 5.2 Panel de conexiones CLP1103. 64 20 grados y otro de -20 grados a los dos segundos, se muestra la respuesta del modelo de Simulink, y la comparación entre la estimación de las galgas y de los acelerómetros. Figura 5.7 Comparación de los valores de simulación; Valor real(verde), Valor estimado por acelerómetros (azul) y Valor estimado por las galgas (rojo). En la Figura 5.7, se presenta una comparación entre los modelos de simulación de la desviación y la desviación estimada en la simulación. La desviación simulada retrata la deformación producida por la barra de cuerpos agrupados, las desviaciones estimadas por los acelerómetros y las galgas que muestran el desplazamiento estimado después de los cálculos con los algoritmos explicados anteriormente. En el gráfico, se puede observar que la estimación realizada por el método de los acelerómetros se comporta como la barra flexible simulada, pero existe un error que se acumula a lo largo del tiempo. Éste error se debe al proceso de integración sobre la velocidad angular calculada para obtener la pendiente. El error se podría reducir notablemente mediante la composición de la señal con la de un giróscopo o algún otro sensor. La desviación estimada por las galgas extensométricas tiene un comportamiento bastante similar al comportamiento del manipulador flexible. La estimación trabaja sin integración de la señal y el hecho de que la señal de las galgas en la simulación no tiene ruido hacen que el error en estado estacionario sea nulo. 4.3.2 Resultados en el Prototipo Experimental En el prototipo experimental el desplazamiento del extremo del manipulador se estima midiendo los tres puentes de galgas y calculándolo mediante el algoritmo previo. 0 0.5 1 1.5 2 2.5 3 -0.1 -0.05 0 0.05 0.1 Time (s) Slope (m) Slope Comparison between Simulation and Estimated Models 0 0.5 1 1.5 2 2.5 3-0.1 -0.05 0 0.05 0.1 Estimated Slope by Accelerometers Simulated Slope Estimated Slope by Gages 65 Para reducir el ruido, las mediciones de deformación se retrasan debido a la utilización de filtros de Kalman. Las tres deformaciones se miden a través de un medio puente de Wheatstone, la Figura 5.8 muestra la deformación medida por cada puente después de una deformación aleatoria. La mayor deformación se obtiene por el puente número uno porque está más cercano al punto de mayor tensión, donde el cilindro actúa sobre el manipulador. Figura 5.8 Gráfica de los valores medidos a la salida de cada puente de Wheatstone. 4.3.3 Comparación Simulación y Prototipo En esta sección, se muestra una comparación a la respuesta escalón de 20 grados, Figura 5.9. El modelo de Simulink se ha rediseñado con los mismos parámetros que la viga utilizada en la prueba experimental. El comportamiento del sistema rediseñado difiere del anterior ya que se acortó la longitud del manipulador, por lo que es más rígida y la frecuencia de funcionamiento se incrementa. Como muestra la figura, se obtiene una frecuencia similar en ambos sistemas. 0 10 20 30 40 50 60 -12 -10 -8 -6 -4 -2 0 2 4x 10 -3 Time (s) Strain (um/m) Filtered Ouput of the three Wheatstone Bridges Bridge 1 Bridge 2 Bridge 3 66 Figura 5.9 Comparación de la deformación estimada por las galgas en la simulación (Gráfica superior) y en el prototipo (Gráfica inferior). La desviación del extremo final muestra que hay una gran deformación al principio del movimiento, seguido de una oscilación cuando el manipulador flexible alcanza la posición deseada en el eje neutral de la barra. 4.4 Resultados del Sistema Controlado 4.4.1 Inclinación sin Controlador El sistema de bucle cerrado es un sistema con una retroalimentación del encoder. El error (la diferencia entre la salida y la entrada deseada) se retroalimenta para influir en el funcionamiento. El movimiento angular del manipulador en radianes se muestra en la Figura 5.10, y el comportamiento obtenido es similar, pero la respuesta experimental es más lenta. Ambas respuestas tienen una vibración debido a la carga parasitaria producida por el manipulador flexible. 0 0.5 1 1.5 2 2.5 3 -0.02 -0.01 0 0.01 0.02 Time (s) Deflection (m) Slope Deflection Simulink 0 0.5 1 1.5 2 2.5 3 -0.3 -0.2 -0.1 0 Time (s) Deflection (m) Slope Deflection Experimental 67 Figura 5.10 Inclinación obtenida sin controlador en simulación (Gráfica superior), y en el prototipo (Gráfica inferior). 4.4.2 Inclinación con Controlador Con el fin de probar el controlador diseñado anteriormente, el controlador se introduce en el sistema manipulador flexible. El objetivo del controlador es proporcionar un comportamiento suave para llegar a la posición deseada Esto se lleva a cabo mediante la modificación de la señal de error producida por las señales de realimentación de referencia, como se muestra en la Figura 5.11. El controlador proporciona la señal para operar la válvula hidráulica de carrete. 0 0.5 1 1.5 2 2.5 3 0 0.1 0.2 0.3 0.4 Time (s) Theta (Rad) Tilt in Radians Simulink 0 0.5 1 1.5 2 2.5 3 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 Time (s) Theta (Rad) Tilt in Radians Experimental 0 2 4 6 8 10 12 14 16 18 20 -2 -1 0 1 2 3 Time (s) Error and Control Signal Experimental Error Control Signal 68 Figura 5.11 Gráficas de entrada y salida del controlador calculado y pendiente obtenida utilizando el controlador. El comportamiento de la respuesta es lento, como muestra la Figura 5.11, de modo que el movimiento suave no hace flexionar el manipulador. Por consiguiente, la salida obtenida es sólo el ruido producido por la cepa galgas mediciones. Por lo tanto, se cambian algunos parámetros para mostrar la reacción de pendiente producida por la entrada escalón de 20 grados. En la Figura 5.12, se puede observar el aumento de la velocidad mediante el aumento del parámetro K del controlador de PT. Con el aumento de la ganancia se obtiene una respuesta rápida, pero se eleva la sobreoscilación. Figura 5.12 Respuesta con controlador de parámetros K=1 y T=0.33. En la Figura 5.13, se expone la estimación de la deformación en el extremo obtenida con el controlador anterior. 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 0.05 0.1 0.15 0.2 0.25 Time (s) Theta (Radians) Step response in Radians K=0.2, T=0.33 0 0.5 1 1.5 2 2.5 3 3.5 4 0 0.1 0.2 0.3 Time Theta (Radians) Step response to controlled K=1 Experimental 69 Figura 5.13 Estimación de la desviación con controlador (K=1, T=0.33). Hay todavía una gran deformación en el arranque, pero menor que el sistema sin el regulador. También la vibración después del movimiento tiene menos amplitud que la obtenida con el sistema no controlado. Por lo tanto, el parámetro T se sintoniza para reducir la sobreoscilación (T=0.2). A continuación, la respuesta a un escalón aún más rápido y sin sobreoscilación. Figura 5.14 Respuesta con controlador de parámetros K=1 y T=0.2. 0 0.5 1 1.5 2 2.5 3 3.5 4 -1 0 1 2 x 10 -4 Time (s) Deflection (m) Slope Deflection from Experimental Test K=1 01234567891011 -0.1 0 0.1 0.2 0.3 0.4 Time (s) Theta (Radians) Step response to controlled T=0.2 70 5 Conclusiones En este proyecto, un manipulador flexible se ha diseñado, simulado, estimado y controlado. También, fue construido un prototipo del manipulador para los trabajos experimentales. Durante el proceso de diseño, se llevó a cabo un análisis paramétrico para investigar cómo la forma y masa del extremo modifican el comportamiento del manipulador flexible. La longitud de la viga y la masa afectan inversamente a la frecuencia. El modelo servo hidráulico tuvo que ser desarrollado al mismo tiempo, debido a que la frecuencia natural del manipulador flexible y la del sistema actuador tuvieron que diferir de algunos órdenes de magnitud para evitar problemas de resonancia. Además, se ha desarrollado un sistema servo hidráulico adaptado con una realimentación de la posición, que se ocupa de manejar y evitar la resonancia parásita. Esa resonancia parásita es producida por la viga flexible. En teoría, la frecuencia del sistema parasito debería diseñarse lo más alta posible, pero en cambio, la amortiguación global se incrementó, para evitar respuestas indeseadas. Por lo tanto, sobre las pruebas de simulación, el método cuerpos agrupados presenta un buen comportamiento como un manipulador flexible, tal y como se esperaba durante el diseño. El modelo simulado funciona con la frecuencia y la deformación similar a la calculada en el estudio paramétrico. Gracias al Toolbox SimMechanic MATLAB, el sistema hidráulico y el sistema de medición se pudieron adjuntar a la viga flexible en un mismo modelo. Durante los trabajos de simulación, se han testado dos sensores y algoritmos de estimación diferentes. El hecho de que las bandas extensométricas no necesitan ningún proceso de integración o filtrado, hace el método más fiable, y además la estimación de la forma del manipulador mediante la ecuación hiperbólica a lo largo de la viga se comporta como una viga flexible. El método de medición acelerómetros en base a la fuerza centrípeta es válida sólo en trayectorias circulares. Por otra parte, la señal de medición es muy ruidosa y necesita ser filtrada, aunque de éste modo, la salida se retrasa. Además, añadiendo un cálculo de cambio de dirección en la velocidad angular se presenta una discontinuidad, por tanto no es un método tan robusto como el de las galgas. Después de hacer varias configuraciones experimentales de hardware y de software, y de rediseño del modelo de barra flexible en Simulink, finalmente, se puso a prueba el método de estimación mediante bandas extensométricas. Los resultados muestran un comportamiento similar entre los sistemas simulados y experimentales, se obtuvo una frecuencia similar en ambos sistemas. Las señales obtenidas de las galgas extensométricas eran muy ruidosas y tuvieron que ser filtradas bastante bien debido a que afectaba a la estimación de la posición, pero a pesar de eso, el comportamiento de la 71 barra se estimó correctamente. Desafortunadamente, debido al programa de utilización del laboratorio el método acelerómetros pudo no ser probado. Como conclusión, la posición final de la estimación de una ecuación hiperbólica con galgas extensométricas es una manera apropiada para estimar la forma de un manipulador flexible, así como los el método de cuerpos agrupados para la simulación de dicha flexibilidad. Esta tecnología se puede aplicar a diferentes formas y diferentes tipos de manipuladores. Por otra parte, se trata de un sistema mecatrónico de bajo coste debido a la ligereza del manipulador y a la simplicidad de los sensores necesarios. Un manipulador ligero siempre necesita un actuador de baja potencia, esto le hace ser más ágil y preciso que cualquier otro tipo de manipulador. Como trabajos futuros, para un mayor desarrollo de ésta tecnología, se puede incluir el método de estimación de la deformación en el control de la posición, añadiendo una realimentación y calculando el error con la posición del extremo deseada. 72 6 Referencias [wikipedia Beronoulli 2013] http://en.wikipedia.org/wiki/Euler%E2%80%93Bernoulli_beam_theory, Marzo 2013. [vlab 2013] http://iitg.vlab.co.in/?sub=62&brch=175&sim=1080&cnt=1, Marzo 2013. [Meirovitch, 1967] Analytical methods in vibrations, by Leonard Meirovitch, Macmillan, 1967. [Wachel & Bates 1976] Wachel, J. C., and Bates, C. L., "Techniques for controlling piping vibration failures", ASME Paper, 76-Pet-18, 1976. [Schiavo, Vigan`o & Ferretti 2006] Object-oriented modelling of flexible beams Francesco Schiavo , Luca Vigan`o and Gianni Ferretti. Springer Science, Business Media B.V. 2006. [wikipedia, Buckling 2013] http://en.wikipedia.org/wiki/Buckling, Abril 2013. [Chudnovsky, Mukherjee, Wendlandt & Kennedy, 2006] Modeling Flexible Bodies in SimMechanics, Victor Chudnovsky, Arnav Mukherjee, Jeff Wendlandt and Dallas Kennedy, August 30, 2006. [Shye & Richardson, 1987] Mass, Stiffness, and Damping Matrix Estimates from Structural Measurements, by Ken Shye and Mark Richardson, Structural Measurement Systems, Inc. San Jose, California, April 6, 1987. [Rydberg 2008] Hydraulic servo systems Karl-Erik Rydberg, Linköpings university (2008). [Mattila 2013] Modelling of a Servo Hydraulic System notes, Jouni Mattila. [Formula Book, 2008] Formula Book for Hydraulics and Pneumatics, Department of Management and Engineering Linköping University (2008). [Jelai & Kroll 2004] Hydraulic Servo-systems, Modelling, Identification and Control, Mohieddine Jelai and Andreas Kroll. Springer (2004). [Korayem, Shafei & Absalan 2012] Estimate a Flexible Link's Shape by the Use of Strain Gauge Sensors, M.H. Korayem, A.M. Shafei and F. Absalan, Mechanical Engineer Department , Iran University of Science and Technology, Narmak Tehran (Iran), August 2012. 73 [Naghshineh, Ameri, Zereshki, Krishnan & Abdoli] Human Motion capture using TriAxial ccelerometers, Sam Naghshineh, Golafsoun Ameri, Mazdak Zereshki and Dr. S. Krishnan, Dr. M. Abdoli-Eramaki. [Abedi, Nadooshan & Salehi 2008] Dynamic Modeling of Tow Flexible Link Manipulators, E. Abedi, A. Ahmadi Nadooshan, and S. Salehi, World Academy of Science, Engineering and Technology 22, 2008 [Kennedy 2009] Force Control of a Hydraulic Servo System, A Thesis presented to the Faculty of the Graduate School, University of Missouri, by Joseph L. Kennedy, Dr. Roger Fales, Thesis Supervisor, May 2009. [Quijano & Passino 2002] A Tutorial Introduction to Control Systems Development and Implementation with dSPACE, Nicanor Quijano and Kevin Passino, Dpt. of Electrical Engineering, The Ohio State University, 2002. Otras Fuentes Consultadas: [Tokhi & Azad 2008] M.O. Tokhi and A.K.M. Azad. Flexible Robot Manipulators Modelling, simulation and control. By The Institution of Engineering and Technology, London, United Kingdom. [Lurie 1952] Lurie, H., "Lateral vibration as related to structural stability", Journal of Applied Mechanics, Vol. 19, pp. 195-204, 1952. [Stokey 2009] Harri's Shock and Vibration Handbook, Chapter 7 Vibration of systems having a distributed mass and elasticity, by William F. Stokey, Mc Graww Hill Professional 2009. [Yan & Zhang 2008] Numerical analysis and control for cantilever flexible beams using PZT patches, Shi Yan *a, Hao Zhang, School of Civil Engineering, Shenyang Jianzhu University, No.9 Hunnan East Road, Hunnan New District, Shenyang, Liaoning, China110168. [Jenkins 2001] Mechanics of Materials Laboratory, Department of Mechanical Engineering University of Washingto. Seattle, Washington. By Michael G. Jenkins, Associate Professor, Mechanical Engineerin, 2001. [Ueno, Ito, Ma & Ikeo 2005] Design of Robust Position/Pressure Controller for Cylinder using Hydraulic Transformer, Tomohiro UENO, Kazuhisa ITO, Weidong MA and Shigeru IKEO. Division of Science and Technology, Graduate school of SOPHIA University. Department of Mechanical Engineering, Faculty of Science and Technology (2005). 80 Mass=50 Kg (32%165) 20 30 40 50 60 70 80 90 100 110 120 2 4 6 8 Mass [Kg] 20 30 40 50 60 70 80 90 100 110 120 0 5 10 15 x 10 7 Frequency [Hz] Tmax [MPa] 81 7.2 Anexo II: Test Bench, Previous Works 7.2.1 Cylinder test bench Design requirements of the cylinder test bench: Natural frequency as low as possible Minimum torque of the joint 1100 Nm Boom angle range –20…100 ° Required joint speed 45 deg/s Joint parameters and the cylinder were selected based on: Max force on the cylinder is high to obtain low natural frequency Higher force (and lower natural freq.) requires bigger piston rod (buckling) Smaller cylinder diameter gives smaller min natural frequency Stroke is minimized; shorter stroke means shorter torque arm and higher force When cylinder is at maximum stroke, there is enough space between joint axle and the cylinder The total length of the cylinder is enough to fit load cell at the end of the piston rod Joint parameters: L1 = 170 mm L2 = 525 mm α αα α1 = 2.5 ° °° ° α αα α2 = 42.5 ° °° ° Boom angle range: -20…100 ° °° ° from horizontal Cylinder: 25/32 – 300 mm Picture 1. Joint parameters Calculation results 82 Minimum natural frequency of the joint: 8.8 Hz (Figure 1.) Maximum force on the cylinder: 11.3 kN (Figure 2.) Minimum torque arm: 81.1 mm (Figure 3.) Required flow for 45 deg/s speed: 6.4 l/min (Figure 4.) Minimum pressure to achieve max force: 177 bar Maximum force of the cylinder (210bar): 13.5 kN Torque at cylinder max force and min torque arm: 1100 Nm (Figure 5.) Cylinder minimum length: 562 mm Figure 1. Natural frequency -20 0 20 40 60 80 9 10 11 12 13 14 15 16 Natural frequency (Hz) Joint angle (deg) Natural frequency 83 Figure 2. Cylinder force with max torque load at every point Figure 3. Torque arm -20 0 20 40 60 80 -2000 0 2000 4000 6000 8000 10000 Force (N) Joint angle (deg) Cylinder force with minimum torque load -20 0 20 40 60 80 90 100 110 120 130 140 150 160 170 Torque arm (mm) Joint angle (deg) Torque arm 84 Figure 4. Required flow for 45 deg/s velocity Figure 5. Joint torque with maximum cylinder force -20 0 20 40 60 80 3.5 4 4.5 5 5.5 6 Flow (liter/min) Joint angle (deg) Required flow for 45deg/s speed -20 0 20 40 60 80 1200 1400 1600 1800 2000 2200 Torque (Nm) Joint angle (deg) Joint torque with max cylinder force 85 7.2.2 Joint construction Maximum radial force at the joint: 11.9 kN Equivalent static bearing load: 5.9 kN Selected bearing: SKF NU 1006 Cylindrical roller bearing Inner diameter: 30 mm Outer diameter: 55 mm Height: 13 mm Static load rating: 17.3 kN Safety factor: 2.9 Safety factor should be higher than 1 for roller bearings. Safety factor for joint axle: 5.2 (Fe 52) Required bending strength of the beam material: 1.2*10^4 mm^3 Selected tubular beam: 100x50-5 Bending strength (x-direction): 3.2*10^4 mm^3 Weight: 10.5kg/m 7.2.3 Strength calculations Safety factor of the boom connection: 3.6 Safety factor of the weakest section between joint axle and cylinder support on the beam: 2.6 Safety factor of the cylinder lower supports: 5.4 Safety factor of the cylinder shafts: 2.0 86 87 Sizing of the cylinder Joint parameters Lengths of hypotenuses Cylinder torque arm when boom is inclined -20deg from horizontal Angle of boom to horizontal Cylinder minimum and maximum lengths Maximum force exerted on the cylinder Force on the cylinder, boom -20deg from horizontal The minimum diameter of the cylinder Mechanical efficiency of the cylinder l1170 mm := l2525 mm := α12.5deg:= α242.5deg:= l1_ l1 cos α1 ( ) := l1_ 170.162 mm = l2_ l2 cos α2 ( ) := l2_ 712.0794 mm = αB20−deg:= st l1 l2 cos α2 ( ) sin αBα2 − ( ) ⋅+       2 l2 cos α2 ( )       cos αBα2 − ( ) ⋅l1tan α1 ( ) ⋅−             2 +:= d l2_ cos 90deg acos l1_ 2 l2_ 2 −st 2 − 2 st⋅l2_ ⋅                 −         ⋅:= d 91.0405 mm = cyl_min l1 l2 cos α2 ( ) sin 20−deg α2 − ( ) ⋅+       2 l2 cos α2 ( )       cos 20−deg α2 − ( ) ⋅l1tan α1 ( ) ⋅−             2 +:= cyl_min 562.4764 mm = cyl_max l1 l2 cos α2 ( ) sin 100deg α2 − ( ) ⋅+       2 l2 cos α2 ( )       cos 100deg α2 − ( ) ⋅l1tan α1 ( ) ⋅−             2 +:= cyl_max 857.0435 mm = Tmax 1100 Nm := Fmax Tmax dcos 20−deg( ) ⋅:= Fmax 11.3539 kN = pref 210 bar := ηm 0.8 := rFmax πpref ⋅ η m ⋅ := r 14.667 mm = 88 The minimum diameter of piston rod Safety factor Modulus of elasticity (AISI 316) Length of the cylinder Proposed cylinder: Maximum allowed force of the cylinder for buckling Minimum pressure to achieve necessary force Fs Piston side of the cylinder considered Maximum force on rod side Maximum force of the cylinder Minimum torque of the joint with maximum cyl force Min torque arm from Matlab Required flow for 45deg/s speed Dc2 r ⋅:= Dc29.3341 mm = n 4 := E 193000000000 Pa := Lc cyl_max := dr 4 Fmaxn⋅Lc 2 ⋅64⋅     π 3 E⋅ := dr24.4397 mm = dR25 mm := DC32 mm := Fbmax π 3 E⋅dR 4 ⋅ n 64⋅Lc 2 ⋅ := Fbmax 12.4314 kN = p_min Fmax4⋅ DC 2 π⋅     η m ⋅ := p_min 176.5 bar = Fminus pref DC 2       2 dR 2       2 −         π⋅         ⋅ η m ⋅:= Fminus 5.2647 kN = Fplus pref DC 2       2 π⋅         ⋅ η m ⋅:= Fplus 13.5114 kN = dmin 81.1 mm := Tmin Fplus d min ⋅:= Tmin 1095.7714 Nm = 89 Minimum change in angle when cylinder moves 1mm (from Matlab) Cylinder area Required flow to achieve 45deg/s speed Minimum natural frequency of the cylinder (stroke value from Matlab) Stroke of the cylinder Stroke at min wn, from Matlab Estimated compression modulus Dead volumes 5% of cyl volumes Linear hydraulic stiffness Torque arm at min wn, from Matlab Torsional stiffness of the arm ∆αmin 0.3367deg:= v 45 deg s ∆αmin 1⋅ mm := v 133.6501 mm s = A1πDC 2       2 ⋅:= A10.0008m 2 = Q v A 1 ⋅:= Q 6.4493 liter min = A2A1πdR 2       2 ⋅−:= A20.0003m 2 = Lscyl_max cyl_min −:= Ls294.5671 mm = y 269.23 mm := Bs2.3 10 9 ⋅ Pa := V10 A1Ls ⋅ 0.05 ⋅:= V20 A2Ls ⋅ 0.05 ⋅:= Kh A1 2 Bs ⋅ y A1 ⋅V10 + A2 2 Bs ⋅ Lsy− ( ) A2 ⋅V 20 + +:= Kh24503779.1991N m = dk111.62 mm := k dk 2 K h ⋅:= k 305293.1829 Nm = I 100kg 1⋅m 2 := ωnk I := 96 Allowed stress (0.67*yield point) Fe37? Cylinder shafts Bending stress Two sides taking the radial force Width of the side wall Diameter of the shaft Bending stress in shaft Comparable stress σv27.04 N mm 2 = σall 2 3220⋅N mm 2 := σall 146.667 N mm 2 = sσall σv := s 5.424= Fr2 F smax 2 := MbFr2 10⋅ mm := Mb0.068 kNm = d 20 mm := Iπ 64 d 4 ⋅:= I 7854mm 4 = e 0.5 d ⋅:= e 10 mm = σb Mbe⋅ I := σb85.9 N mm 2 = Aτπ 4d 2 ⋅= τFr Aτ := τ37.8 N mm2 = σvσb 2 3τ 2 ⋅+:= σv108 N mm2 = σall 0.67320⋅N mm 2 := 97 Yield point of the shaft material (Fe 52) = 320N/mm^2 Allowed bending stress of the material Safety factor σall 214.4 N mm 2 = nσ all σv := n 1.99= 98 7.3 Anexo III: SimMechanics Model SimMechanics Model based on the Lumped-Parameter Method which discretizes the beam length in n beam elements. Each element is a body joint body combination: Figure 5 Also Spring and Damper parameters are added to the joint. Figure 6 Spring-Damper Figure 7 99 The normalized moment equation: l=+2#9m=+9 ) ==4n:=k9n=4E Where ω o2 = k/I and I is the moment of inertia. The damping coefficient 2ξω o is a quasi-empirical value that accounts for energy lost to visco-elastic effects. In SimMechanics, the relative angle θ n and joint velocity θn m are measured at the joint using a Joint Sensor block. Then k is multiplied by the angle, and a material damping coefficient 2ξω o is multiplied by the angular velocity (Figure 7). The two resulting moments are added and applied back to the joint using a joint actuator. Entire Model Figure 8 Global view of the entire simulation model, where there are some blocks in different colors (Figure 8). The black one is just a simulink block where the trigonometric operations are done to adapt the input (theta) to the hydraulic circuit input (stroke displacement). The red one is the Hydraulic circuit block, developed using simHydraulic components, it will be explained later. The blue blocks are the mechanical part of the system composed by simMechanics bodies and joints. The actuator block has been done to connect the cylinder to the beam, 100 which is composed by two bodies, one attached to the beam and the other to the ground according to the new test bench. This bodies are connected each other by a prismatic joint (Figure 9). Figure 9 In the flexible beam block, the flexible manipulator is simulated by the lumped bodies shown previously. The beam is connected to the actuator and the ground by single rotational joints, according to the test bench. The first body, which is located between the two joints (actuator and ground), is rigid and not lumped as the followers. All bodies have similar properties to the proposed flexible beam, as inertia, young modulus, shape, density and bending moment. Hydraulic Circuit 101 Figure 10 The hydraulic circuit is composed by a proportional valve with 4 ways and 3 positions and several manometers to measure cylinder and supply pressure (Figure 10). The input is adapted to the valve spool position. Simulation Figure 11 102 The tip end position (X, Y, Z) is shown in Figure 12 after a step input of 20 degrees: Figure 12. Tip end coordinates (X yellow, Y purple, Z blue) 103 7.4 Anexo IV: Hollow Section Profile 104 105 112 h^2 + fBeamElement.length^2 ] ; % principalMoments = [ (4/3)*material.density*fBeamElement.length*e*(h1^3+b1^3+(h1-e)^3+(b1e)^3) % (1/12) * fBeamElement.mass * ( fBeamElement.length^2 + b^2 ) % (1/12) * fBeamElement.mass *( h^2 + fBeamElement.length^2 ) ]; fBeamElement.inertia = diag( principalMoments ); principalMoments = (1/3) * fBeamElement.mass * [ h^2 + b^2 - (b2*e)^2 - (h-2*e)^2 % http://www.prt.fernunihagen.de/lehre/KURSE/PRT001/course_main/node25.html (L1^2) + b^2 h^2 + (L1^2) ] ; fBeamElement.inertia1 = diag( principalMoments ); yzyBendingMoment = ( b * h^3 - (b-2*e)*(h-2*e)^3) / 12 yzz = fBeamElement.length*h^3/12 springConstantAtTip = 3 * material.youngsModulus * yzyBendingMoment / ( fBeamElement.length^3 ); fBeamElement.EI=material.youngsModulus * yzyBendingMoment fBeamElement.matDamping =49; fBeamElement.material = material; %%%% material.youngsModulus * yzyBendingMoment/fBeamElement.length; 7.6.3 Hydraulic System Data %% Hydraulic Cylinder Model %based on the Jouni Mattila's notes %Data Rod.L=Lenght; Rod.L1=Lenght_1; Cyl.L=Stroke; %[m] Cyl.D1=Piston_Diameter; %[m] Cyl.D2=Rod_Diameter; %[m] Cyl.A1=(Cyl.D1^2)/4*3.1416; Cyl.A2=Cyl.A1-(Cyl.D2^2)/4*3.1416; Cyl.V0 = Dead_Volume; % Dead Volume Cyl.Vt=Cyl.L*Cyl.A1; Cyl.V1=Cyl.Vt/2; % Linearisation point at midstroke Cyl.V2=Cyl.Vt/2-(Cyl.D2^2)/4*3.1416*Cyl.L/2; Load.M=Rod.L*TipMass*9.8/Rod.L1*cos(3.1416/9); % [3000Kg] Cyl.B=Bulk; % [1e6] Cyl.Cp=0; Cyl.Fv=Viscous_Force; %[N] Valve.QN=24*60/1000; %[m3/s] 24 l/min Valve.Xv_max=0.01; 113 Valve.dpN=3.5e6; %[MPa] 35 bar Valve.Tau = 10e-3; % [s] Supply.pP = Supply_Pressure; % [100e5Pa] Supply.pT = Tank_Pressure; % [Pa] %% Initial conditions Init.p1 = Load.M/Cyl.A1; % [Pa] Init.p2 = 0; % [Pa] Init.x = Initial_Position_Cyl; % [m] %% Linearization point %Valve coefs for pl=0 and Xv1=0 Valve.Kc1=0; %%%%% Valve.Kc2=Valve.Kc1; %%%%% Valve.Kq1=Valve.QN/Valve.Xv_max*sqrt((Supply.pP)/2/Valve.dpN); Valve.Kq2=Valve.Kq1; Valve.Kp=0.8*Supply.pP/0.05/Valve.Xv_max; %% System Forwards; Output Velocity NUM1=Valve.Kq1/Cyl.A1; DEN1=[Cyl.Vt*Load.M/4/Cyl.B... (Load.M*(Valve.Kc1+Cyl.Cp)+Cyl.Vt*Cyl.Fv/4/Cyl.B)/(Cyl.A1)^2 1 0]; %for Output position % DEN1=[Cyl.Vt*Load.M/4/Cyl.B... % (Load.M*(Valve.Kc1+Cyl.Cp)+Cyl.Vt*Cyl.Fv/4/Cyl.B)/(Cyl.A1)^2 1 0]; Cyl_sys=tf(NUM1,DEN1); Cyl_sys=minreal(Cyl_sys); %% System Forwards; Output Load Pressure NUM2=Valve.Kq2; DEN2=[Cyl.Vt*Load.M/4/Cyl.B Load.M*Valve.Kc1+Cyl.Vt*Cyl.Fv/4/Cyl.B Valve.Kc1*Cyl.Fv+Cyl.A1^2]; Cyl_sys2=tf(NUM2,DEN2); Cyl_sys2=minreal(Cyl_sys2); %% State-space Model %X'=A*X+B*U X=[x' p1' p2'] U=Xv A=[-Cyl.Fv/Load.M Cyl.A1/Load.M -Cyl.A2/Load.M;... -Cyl.B*Cyl.A1/Cyl.V1 -Cyl.B*Valve.Kc1/Cyl.V1 0;... Cyl.B*Cyl.A2/Cyl.V2 0 Cyl.B*Valve.Kc2/Cyl.V2]; B=[0 Cyl.B*Valve.Kq1/Cyl.V1 -Cyl.B*Valve.Kq2/Cyl.V1]'; 114 C=[1 0 0]; D=0; [NUM,DEN]=ss2tf(A,B,C,D); System_sys=tf(NUM,DEN); System_sys=minreal(System_sys) % bode(System_sys); [Wn,Z,P]=damp(System_sys) %Natural frequency, zeros and poles % hold on; System_sys=feedback(System_sys,1) % bode(System_sys); % legend('Before Feedback','After Feedback'); [Wn,Z,P]=damp(System_sys) %% Kh1=Cyl.A1^2*Cyl.B/(Cyl.V1); Kh2=Cyl.A2^2*Cyl.B/(Cyl.V2); Kh=Kh1+Kh2; % Stiffness of the Cylinder [MN] M_cyl=Load.M; Wnat=sqrt(Kh/M_cyl) % Natural frequency Kqa=Valve.Kq1/Cyl.A1; % Velocity Gain a=Rod.L1*cos(-3.1416/9); J=mass*L^2/4+TipMass*L^2; W=sqrt(Kh*a^2/J) 7.6.4 Estimation functions %% Hyperbolic and Parabolic Estimation % Coefficient Calculations % coef=[A B C D E F G] X=[0 L x1 x2 x3] % % Strain_pos = [0, Lenght_1, Lenght, x1, x2, x3]; % Hyp_C = [0 0 0 0 0 0 0.07 0.045 0.028]; [Hyp_A] = Polinomic_Eq_7Coeff_Acc(Strain_pos7) % A*coef=C % [X,FVAL,EXITFLAG] = fsolve('Hyperbolic_Equations',x0) Hyp_APA_inv = inv(Hyp_A); [Hyp_A] = Polinomic_Eq_7Coeff(Strain_pos7) % A*coef=C % [X,FVAL,EXITFLAG] = fsolve('Hyperbolic_Equations',x0) Hyp_APG_inv = inv(Hyp_A); [Hyp_A] = Hyperbolic_Eq_9Coeff(Strain_pos) % A*coef=C % [X,FVAL,EXITFLAG] = fsolve('Hyperbolic_Equations',x0) Hyp_AHG_inv = inv(Hyp_A); 115 [Hyp_A] = Hyperbolic_Eq_9Coeff_Acc(Strain_pos) % A*coef=C % [X,FVAL,EXITFLAG] = fsolve('Hyperbolic_Equations',x0) Hyp_AHA_inv = inv(Hyp_A); %%%%%%%%%%%%%%%%%%%%%%%%% Hyperbolic_Eq_9Coeff .m; function [A]= ecuaciones(x) % coef=[A B C D E F G H I] X=[0 L1 L x1 x2 x3] A = [sinh(x(1)) sinh(2*x(1)) sinh(3*x(1)) sinh(4*x(1)) cosh(x(1)) cosh(2*x(1)) cosh(3*x(1)) cosh(4*x(1)) 1; %% w(0)=0 w''(0)=0 w(L1)=0 w''(L1)=0 w''(L)=0 w'''(L)=0 sinh(x(1)) 4*sinh(2*x(1)) 9*sinh(3*x(1)) 16*sinh(4*x(1)) cosh(x(1)) 4*cosh(2*x(1)) 9*cosh(3*x(1)) 16*cosh(4*x(1)) 0; sinh(x(2)) sinh(2*x(2)) sinh(3*x(2)) sinh(4*x(2)) cosh(x(2)) cosh(2*x(2)) cosh(3*x(2)) cosh(4*x(2)) 1; %% sinh(x(2)) 4*sinh(2*x(2)) 9*sinh(3*x(2)) 16*sinh(4*x(2)) cosh(x(2)) 4*cosh(2*x(2)) 9*cosh(3*x(2)) 16*cosh(4*x(2)) 0; sinh(x(3)) 4*sinh(2*x(3)) 9*sinh(3*x(3)) 16*sinh(4*x(3)) cosh(x(3)) 4*cosh(2*x(3)) 9*cosh(3*x(3)) 16*cosh(4*x(3)) 0; cosh(x(3)) 8*cosh(2*x(3)) 27*cosh(3*x(3)) 64*cosh(4*x(3)) sinh(x(3)) 8*sinh(2*x(3)) 27*sinh(3*x(3)) 64*sinh(4*x(3)) 0; -(sinh(x(4)))/25 -(4*sinh(2*x(4)))/25 -(9*sinh(3*x(4)))/25 - (16*sinh(4*x(4)))/25 -(cosh(x(4)))/25 -(4*cosh(2*x(4)))/25 - (9*cosh(3*x(4)))/25 -(16*cosh(4*x(4)))/25 0; -(sinh(x(5)))/25 -(4*sinh(2*x(5)))/25 -(9*sinh(3*x(5)))/25 - (16*sinh(4*x(5)))/25 -(cosh(x(5)))/25 -(4*cosh(2*x(5)))/25 - (9*cosh(3*x(5)))/25 -(16*cosh(4*x(5)))/25 0; -(sinh(x(6)))/25 -(4*sinh(2*x(6)))/25 -(9*sinh(3*x(6)))/25 - (16*sinh(4*x(6)))/25 -(cosh(x(6)))/25 -(4*cosh(2*x(6)))/25 - (9*cosh(3*x(6)))/25 -(16*cosh(4*x(6)))/25 0; ]; 116 7.7 Anexo VII: Test Bench Drawing 117 7.8 Anexo VIII: Finite Element Analysis with MATLAB Code modification to project boundary conditions from: MATLAB Codes for Finite Element Analysis, Solids and Structures, A. J.M. Ferreira, Universisdade do Porto, Springer 2009 %% FEM (directly from the example file provided w. the book) %--------------------------------------------------------------------- ----- % Example 8.11.1 % to find the transient response of a cantilever beam with a tip load % % Problem description % Find the transient response of a cantilever beam whose length is 1 m. % long. The beam has the cross-section of 0.02 m by 0.02 m and the mass % density is 1000 Kg/m^3. The elastic and shear modulus is 100 GPAa and % 40 GPa, respectively. Use 4 elements. % % Variable descriptions % k = element stiffness matrix % kk = system stiffness matrix % m = element mass matrix % mm = system mass matrix % index = a vector containing system dofs associated with each element % bcdof = a vector containing dofs associated with boundary conditions %--------------------------------------------------------------------- --- clear all clc nel=4; % number of elements nnel=2; % number of nodes per element ndof=2; % number of dofs per node nnode=(nnel-1)*nel+1; % total number of nodes in system sdof=nnode*ndof; % total system dofs el=100*10^9; % elastic modulus tleng=1; % total beam length leng=tleng/nel; % same size of beam elements xi=0.02^4/12; % height (or thickness) of the beam area=0.004; % cross-sectional area of the beam rho=1000; % mass density of the beam ipt=1; % option flag for mass matrix dt=0.0001; % time step size ti=0; % initial time tf=0.2; % final time nt=fix((tf-ti)/dt); % number of time steps 118 nbc=2; % number of constraints bcdof(1)=1; % transverse displ. at node 1 is constrained bcdof(2)=2; % slope at node 1 is constrained kk=zeros(sdof,sdof); % initialization of system stiffness matrix mm=zeros(sdof,sdof); % initialization of system mass matrix force=zeros(sdof,1); % initialization of force vector index=zeros(nel*ndof,1); % initialization of index vector acc=zeros(sdof,nt); % initialization of acceleartion matrix vel=zeros(sdof,nt); % initialization of velocity matrix disp=zeros(sdof,nt); % initialization of displ. matrix vel(:,1)=zeros(sdof,1); % initial zero velocity disp(:,1)=zeros(sdof,1); % initial zero displacement force(9)=100; % tip load of 100 for iel=1:nel % loop for the total number of elements index=feeldof1(iel,nnel,ndof); % extract system dofs associated with element [k,m]=febeam1(el,xi,leng,area,rho,ipt); % compute element stiffness matrix kk=feasmbl1(kk,k,index); % assemble each element matrix into system matrix mm=feasmbl1(mm,m,index); % assemble each element matrix into system matrix end % original simulation mminv=inv(mm); % invert the mass matrix % central difference scheme for time integration for it=1:nt acc(:,it)=mminv*(force-kk*disp(:,it)); for i=1:nbc ibc=bcdof(i); acc(ibc,it)=0; end vel(:,it+1)=vel(:,it)+acc(:,it)*dt; disp(:,it+1)=disp(:,it)+vel(:,it+1)*dt; end acc(:,nt+1)=mminv*(force-kk*disp(:,nt+1)); figure time=0:dt:nt*dt; plot(time,disp(9,:)) xlabel('Time(seconds)') 119 ylabel('Tip displ. (m)') title('Tip displacement according to the book') grid on %% SS analysis % The model is made up of 4 beam elements (5 nodes,10variables) % use the system massand stifness matrices % calculated above M = mm; K = kk; % Apply constraints. % e.g. if the first node is constrained, set K(1,:) = 0 K(1:2,:) = 0; % add some damping tdamp = 5; % tip damping damping = [16*tdamp 0 8*tdamp 0 4*tdamp 0 2*tdamp 0 tdamp 0]; D = diag(damping); D = 0*D; % Damping removed.. comment or uncomment.. % input influence matrix F = zeros(sdof); for i=1:2:sdof F(i,i) = 1; % only forces allowed; no torgues end % create a state space model [ A, B ] = FEM_2nd_to_SS( M, D, K, F ); clear D C = [eye(sdof,sdof) zeros(sdof,sdof)]; D = zeros(sdof,sdof); sys = ss(A,B,C,D); % sys.StateName = {'q1' 'q2' 'q1dot' 'q2dot'}; % sys.InputName = {'f1' 'f2'}; % sys.OutputName = {'endpos'}; %sys % simulation t = 0:0.1e-4:0.2; u = zeros(sdof,length(t)); u(sdof-1,:) = 100; % 100N load at the tip x0 = zeros(sdof*2,1); [Y,T,X] = lsim(sys,u,t,x0,'zoh'); figure plot(T,X(:,sdof-1)) xlabel('Time(seconds)') ylabel('Tip displ. (m)') title('Tip displacement according to CST') grid on % compare solutions figure plotyy(time,disp(end-1,:),T,X(:,sdof-1)) 120 title('Comparison') legend('kirjan','CST') grid on % export matrices save Amatrix.txt A -ascii -double -tabs save Bmatrix.txt B -ascii -double -tabs save Cmatrix.txt C -ascii -double -tabs save Dmatrix.txt D -ascii -double -tabs function [kk,ff]=feaplyc2(kk,ff,bcdof,bcval) %---------------------------------------------------------- % Purpose: % Apply constraints to matrix equation [kk]{x}={ff} % %----------------------------------------------------------- n=length(bcdof); sdof=size(kk); for i=1:n c=bcdof(i); for j=1:sdof kk(c,j)=0; end kk(c,c)=1; ff(c)=bcval(i); end function [kk]=feasmbl1(kk,k,index) %---------------------------------------------------------- % Purpose: % Assembly of element matrices into the system matrix % %----------------------------------------------------------- edof = length(index); for i=1:edof ii=index(i); for j=1:edof jj=index(j); kk(ii,jj)=kk(ii,jj)+k(i,j); end end function [k,m]=febeam1(el,xi,leng,area,rho,ipt) %------------------------------------------------------------------- % Purpose: % Stiffness and mass matrices for Hermitian beam element % nodal dof {v_1 theta_1 v_2 theta_2} % %------------------------------------------------------------------- 121 % stiffness matrix c=el*xi/(leng^3); k=c*[12 6*leng -12 6*leng;... 6*leng 4*leng^2 -6*leng 2*leng^2;... -12 -6*leng 12 -6*leng;... 6*leng 2*leng^2 -6*leng 4*leng^2]; % consistent mass matrix if ipt==1 mm=rho*area*leng/420; m=mm*[156 22*leng 54 -13*leng;... 22*leng 4*leng^2 13*leng -3*leng^2;... 54 13*leng 156 -22*leng;... -13*leng -3*leng^2 -22*leng 4*leng^2]; % lumped mass matrix elseif ipt==2 m=zeros(4,4); mass=rho*area*leng; m=diag([mass/2 0 mass/2 0]); % diagonal mass matrix else m=zeros(4,4); mass=rho*area*leng; m=mass*diag([1/2 leng^2/78 1/2 leng^2/78]); end function [index]=feeldof1(iel,nnel,ndof) %---------------------------------------------------------- % Purpose: % Compute system dofs associated with each element in one- % dimensional problem % %----------------------------------------------------------- edof = nnel*ndof; start = (iel-1)*(nnel-1)*ndof; for i=1:edof index(i)=start+i; end function [ A, B ] = FEM_2nd_to_SS( M, D, K, F ) % The input should be a linearized 2nd order system % of the form: % M*q´´ + D*q´ + K*q = F*u