Análisis mecánico de la cabeza de rotor en una aeronave con capacidad de despegue y aterrizaje vertical (VTOL)

i Equation Chapter 1 Section 1 Proyecto Fin de Máster Ingeniería Industrial Análisis mecánico de la cabeza de rotor en una aeronave con capacidad de despegue y aterrizaje vertical (VTOL) Autor: José Antonio Beardo Duque Tutores: Daniel García Vallejo y Sergio Esteban Roncero Dpto. Ingeniería Mecánica y Fabricación Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2021 iii Proyecto Fin de Máster Ingeniería Industrial Análisis mecánico de la cabeza de rotor en una aeronave con capacidad de despegue y aterrizaje vertical (VTOL) Autor: José Antonio Beardo Duque Tutores: Daniel García Vallejo Sergio Esteban Roncero Dpto. Ingeniería Mecánica y Fabricación Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2021 v Proyecto Fin de Máster: Análisis mecánico de la cabeza de rotor en una aeronave con capacidad de despegue y aterrizaje vertical (VTOL) Autor: José Antonio Beardo Duque Tutores: Daniel García Vallejo y Sergio Esteban Roncero El tribunal nombrado para juzgar el Proyecto arriba indicado, compuesto por los siguientes miembros: Presidente: Vocales: Secretario: Acuerdan otorgarle la calificación de: El Secretario del Tribunal Sevilla, 2021 A Valentín I Agradecimientos on este proyecto termina una de las etapas más duras y enriquecedoras de mi vida que sólo ha sido posible culminar gracias al apoyo de todos los que en algún momento han estado a mi lado. Estos más de 7 años en la Escuela no habrían sido posibles sin todas las personas que me han acompañado. A mis amigos, por aguantar mis insufribles momentos de incertidumbre y por levantarme el ánimo en las circunstancias difíciles. A mis compañeros, con los que he pasado interminables tardes en la biblioteca. Su ayuda ha sido inestimable y gracias a ellos he podido llegar hasta aquí. A mi familia, a mi madre, que nunca han perdido la paciencia y cuyo apoyo incondicional ha sido en gran parte la causa de mis logros. A Dani y Sergio, por su gran labor y esfuerzo a lo largo del desarrollo de este trabajo. José Beardo Duque Sevilla, 2021 C 6.2 Modelo de Findley 35 6.2.1 Cálculo del número de ciclos hasta el fallo 37 6.3 Modelo SWT 38 6.4 Modelo Fatemi-Socie 38 6.5 Modelos Uniaxiales 40 6.5.1 Método uniaxial puro 41 6.5.2 Método multiaxial simplificado 41 7 Aplicación al Modelo 43 7.1 Punto crítico 43 7.2 Análisis de la aplicación de esfuerzos 44 7.2.1 Compresión ejercida por el tornillo 44 7.2.2 Reestructuración de las zonas de aplicación de los esfuerzos 45 7.3 Otros aspectos a tener en cuenta 46 7.3.1 No-linealidad geométrica 46 7.3.2 Fuerza centrífuga 47 7.4 Fases de vuelo 47 7.4.1 Cálculo de los parámetros de vuelo 47 7.5 Procedimiento de daño acumulado 51 7.5.1 Definición de los ciclos de fatiga 51 7.6 Amortiguamiento y vibraciones 55 7.6.1 Amortiguamiento 55 7.6.2 Vibraciones y transitorios 56 7.7 Resultados 59 7.7.1 Punto crítico de la fase de despegue (PCD) 60 7.7.2 Punto crítico del resto de fases (PCR) 61 8 Conclusiones 63 9 Futuras Líneas de Investigación 65 Anexo A 67 Anexo B 69 Anexo C 75 Código del modelo de fuerzas de la co-simulación 75 Código del modelo de fuerzas para el cálculo de parámetros 80 Anexo D 85 D.1 Remote Point 85 D.2 Mesh 86 D.3 Transient 87 D.4 Solution 88 Referencias 89 IX ÍNDICE DE TABLAS Tabla 2-1. Masa y centros de masa de los elementos según el programa Solid Edge. 10 Tabla 3-1. Definición de conceptos relevantes. 11 Tabla 4-1. Propiedades básicas del material empleado. 24 Tabla 6-1. Parámetros y valores empleados. 39 Tabla 6-2. Datos necesarios para la obtención de la vida a fatiga mediante el método uniaxial puro. 41 Tabla 7-1. Valores empleados para el cálculo del apriete del tornillo. 45 Tabla 7-2. Resumen de la aplicación de fuerzas. 46 Tabla 7-3. Parámetros de vuelo iniciales. 48 Tabla 7-4. Parámetros de vuelo en la fase de despegue. 49 Tabla 7-5. Parámetros de vuelo en la fase de crucero nominal. 49 Tabla 7-6.Parámetros de vuelo en la fase de crucero lento. 50 Tabla 7-7. Parámetros de vuelo de la 1º parte de la transición. 50 Tabla 7-8.Parámetros de vuelo de la 2º parte de la transición. 51 Tabla 7-9.Parámetros de vuelo de la fase crucero ante aterrizaje. 51 Tabla 7-10. Frecuencias naturales de los dos primeros modos de vibración del elemento. 55 Tabla 7-11. Resultados del análisis de vida a fatiga para cada fase individual. 59 Tabla 7-12. Secuencia de fases empleada para el cálculo de vida a fatiga. 60 Tabla 7-13. Ángulos de Euler y ángulo del plano crítico para el punto PCD. 61 Tabla 7-14. Ángulos de Euler y ángulo del plano crítico para el punto PCR. 61 Tabla A-1. Tensiones del punto en crítico en cada fase estudiada de forma aislada. 67 Tabla C-1. Parámetros empleados en TCMEP. 76 XI ÍNDICE DE FIGURAS Figura 1-1. Modelo 3D del vehículo en desarrollo [1]. 1 Figura 1-2. Modos de vuelo principales de una aeronave tilt rotor [3]. 2 Figura 1-3. XV-15 (izq) y V-22 (dcha). 3 Figura 1-4. TR-60 (izq) y AKITSU (dcha). 3 Figura 1-5. ROTORwing. 4 Figura 1-6. Diagrama de flujo del procedimiento expuesto. 6 Figura 2-1. Rotor MEVP22. 7 Figura 2-2. Representación gráfica de los movimientos principales de las palas. 8 Figura 2-3. Objeto 3D (extrusión) a partir del boceto sobre el plano de referencia. 8 Figura 2-4. Relación de alineación axial del elemento ‘cabeza’ con la pieza ‘mástil fino’ (izq) y resultado del establecimiento de todas las relaciones de la pieza ‘cabeza’ (dcha). 9 Figura 2-5. Conjunto completo modelado en Solid Edge (izq) y detalle de las barras de control (dcha). 9 Figura 2-6. Dimensiones de la pala empleada [mm] (izq) y juegos de palas incluidas en el sistema (dcha). 10 Figura 3-1. Esquema simplificado del flujo de información de la co-simulación. 11 Figura 3-2.Ventana de importación en Adams View. 12 Figura 3-3. Herramienta ‘Marker’ (resaltada en amarillo) empleada para la definición de los centros de masa de los elementos. 12 Figura 3-4. Ventana de modificación de elemento. 13 Figura 3-5. Herramienta ‘Joints’. 13 Figura 3-6. Restricción fija entre la pala y la pieza ‘Unión’. 14 Figura 3-7. Movimientos aplicados en el sistema. 14 Figura 3-8. Configuración de FZ en la restricción fija existente entre la pala y su anclaje. 15 Figura 3-9. Ventana de configuración de la simulación. 15 Figura 3-10. Herramienta ‘Six component force’. 16 Figura 3-11. Variable ‘input’ (izqda.) y variable ‘output’ (dcha). 16 Figura 3-12. Archivos generdos en la exportación de Adams. 18 Figura 3-13. Bloque mecánico en Simulink. 18 Figura 3-14. Modelo de fuerzas (izqda.) y modelo mecánico (dcha). 19 Figura 3-15. Punto de aplicación para las condiciones de 1 fuerza por pala. 19 Figura 3-16.Disposición de las fuerzas en Adams para la simulación de 17 puntos de aplicación por pala (izq) y visualización de los puntos de aplicación sobre la pala (dcha). 20 Figura 3-17. Esfuerzos obtenidos para los distintos números de puntos de aplicación por pala. 20 Figura 4-1. Módulo ‘Transient Structural’ de Ansys Workbench. 23 Figura 4-2. Representación gráfica de la resistencia a fatiga de algunos componentes de la familia 7000 de aluminio (obtenido del programa CESEdupack). 24 Figura 4-3. Geometría de la pieza analizada. 25 Figura 4-4. Módulos de la ventana ‘Mechanical’. 25 Figura 5-1.’Remote Point’ (círculo amarillo) y caras de aplicación (en rojo). 27 Figura 5-2. Sección de la pieza mallada utilizando la función ‘Hexa Dominant’. 28 Figura 5-3. Mallado hexaédrico y resultado de la aplicación de la herramienta ‘Sizing’ en la zona sombreada. 29 Figura 5-4. Caras de aplicación de las restricciones cilíndricas (flechas amarillas) y redondeo de 1 mm sobre las aristas (flechas rojas). 29 Figura 5-5. Variación de las tensiones máximas en función del número de nodos empleados en la malla hexaédrica. 30 Figura 5-6. Variación de las tensiones máximas en función del número de nodos empleados en la malla tetraédrica. 30 Figura 5-7. Malla tetraédrica empleada en el estudio. 31 Figura 5-8. Módulo ‘Solution’ en Ansys Workbench. 31 Figura 5-9. Variación de la velocidad de avance máxima con el ángulo de paso (en grados). 32 Figura 6-1. Representación genérica de los círculos de Mohr [9]. 36 Figura 6-2. Posición del punto crítico obtenido, señalado con flecha (izq), y resultado del cálculo del parámetro de daño (dcha). 37 Figura 6-3.Curva S-N empleada. 38 Figura 6-4. Daño y número de ciclos hasta el fallo obtenidos para cada plano del punto crítico mediante el criterio de Fatemi-Socie. 40 Figura 7-1. Representación gráfica del parámetro de daño escogido. 43 Figura 7-2. Tornillo M2 de sujeción de la pala. 44 Figura 7-3. Zona de aplicación de la fuerza de apriete del tornillo. 45 Figura 7-4. Orificio de los tornillos (a), caras superiores (b) y caras inferiores (c). 46 Figura 7-5. Diferencial de empuje a lo largo de la pala para la fase de crucero nominal. 50 Figura 7-6. Secuencia de los ángulos de Euler en Ansys [16]. 52 Figura 7-7. Diagrama de picos y valles para la plicación del procedimiento de Rainflow [4]. 54 Figura 7-8. Representación de los ciclos de carga en la curma tensión-deformación [4]. 54 Figura 7-9. Visualización de los dos primeros modos de vibración del elemento. 55 Figura 7-10. Pestaña ‘Damping Controls’ en Ansys. 56 Figura 7-11. Señal original. 56 Figura 7-12. Señal en el dominio de la frecuencia. 57 Figura 7-13. Resultado de la aplicación de la ventana de Hanning. 58 Figura 7-14. Compración de la señal original y la antitransformada tras la aplicación de la ventana de Hanning y la supresión de las frecuencias superiores a la mitad de la frecuencia máxima de muestreo. 58 Figura 7-15. Punto crítico de la fase de despegue. 60 Figura 7-16. Punto crítico del resto de fases. 60 Figura 7-17. Diagramas Rainflow empleados en el punto crítico de la fase de despegue. 61 Figura 8-1.Diagrama de Rainflow incluyendo una fase adicional tras la fase de crucero lento. 63 Figura D-1. Árbol de objetos emlpeado en Ansys Workbench. 85 XIII Figura D-2. Definición del ‘Remote Point’. 86 Figura D-3. Modificaciones en la pestaña ‘Sizing’ del apartado ‘Mesh’. 86 Figura D-4. Ventana de configuración (izq) y zonas de aplicación (dcha) de la herramienta ‘Face Sizing’. 86 Figura D-5. Pestaña ‘Statistics’. 86 Figura D-6. Ventana ‘Step Controls’. 87 Figura D-7. Ventana de definición de ‘Cylindrical Support’. 87 Figura D-8. Definición como vector (izq) y en forma de tabla (dcha) de las fuerzas/momentos. 87 Figura D-9. Selección del Remote Point como punto de aplicación de esfuerzos. 88 Figura D-10. Ventana de configuración del resultado ‘Stress’. 88 Figura D-11. Definición del resultado ‘User Defined result’. 88 1 1 INTRODUCCIÓN n este documento se reflejan los resultados obtenidos del estudio y análisis de la vida a fatiga multiaxial de la zona de anclaje de la pala de un rotor de pequeña envergadura de tipo paso variable enmarcado en las investigaciones y análisis realizados para el desarrollo del proyecto EMERGENTIa. Para su realización se ha diseñado un rotor tomando como base un modelo real y se han estudiado, mediante el empleo de simulaciones y procedimientos analíticos, los diferentes factores que intervienen en el estado tensional necesario para el cálculo de la vida a fatiga, obteniéndose como resultado un procedimiento inicial para su determinación. 1.1 Proyecto EMERGENTIa EMERGENTIa (DevElopment of an unManned convERtible aircraft for rapid and efficient deployment in emerGENcy situaTIons) es un proyecto multidisciplinar e internacional que surge del auge y la necesidad de vehículos no tripulados de intervención rápida para tareas de búsqueda y rescate. Actualmente se emplean vehículos de este tipo en multitud de escenarios y situaciones catastróficas, sin embargo, las tecnologías más empleadas se basan en aeronaves de ala fija o rotatoria mientras que el objetivo de este proyecto es el diseño y fabricación de un vehículo que incorpore ambas. Los dispositivos de ala fija poseen características que los hacen ideales para vuelos de larga distancia con cargas relativamente pesadas, siendo su mayor desventaja la necesidad de infraestructuras para el despegue y el aterrizaje (catapultas, pistas…). Por otra parte, las aeronaves de ala rotatoria como helicópteros o quadcópteros reducen esta necesidad, pero sus capacidades están limitadas en lo referido al vuelo horizontal y carga. La propuesta del proyecto EMERGENTIa (Figura 1-1) se basa en un tercer tipo de aeronave de ala fija que además posee la capacidad de despegue y aterrizaje horizontal (VTOL según sus siglas en inglés). Para incorporar estas capacidades se emplea la tecnología ‘tilt rotor’ (rotor basculante) que otorga la posibilidad de variar la orientación del eje de giro de las palas en función de la fase de vuelo, siendo su orientación vertical durante el despegue y el aterrizaje, y horizontal para el vuelo a gran velocidad. Esta tecnología consigue aunar las ventajas de los dos tipos de aeronaves principales mencionados y dota a la aeronave de una gran versatilidad, sin embargo, su uso es casi inexistente debido al riesgo inherente a este tipo de tecnología, por lo que se requiere de un análisis exhaustivo y multidisciplinar que asegure su idoneidad. Figura 1-1. Modelo 3D del vehículo en desarrollo [1]. Además de esta característica y a partir de estudios llevados a cabo en conjunción con compañías del sector de E Introducción 2 las emergencias, se definen los siguientes requisitos: • Tamaño. Las características dimensionales de la aeronave no deben sobrepasar el espacio disponible en los vehículos de emergencia para que su transporte se desarrolle de manera rápida. • Versatilidad. Las tareas desempañadas por el vehículo deben abarcar un amplio rango de operaciones. • Simplicidad de operación. La aeronave debe poder ser manipulada y operada sin necesidad de un entrenamiento exhaustivo. • Cumplimiento de la normativa y requerimientos de aeronavegabilidad. Información más detallada sobre el proyecto, así como la descripción de sus componentes y las principales áreas de investigación se pueden encontrar en la propuesta científica desarrollada por el Departamento de Ingeniería Aeroespacial de la Universidad de Sevilla [2]. 1.2 Estado del arte 1.2.1 Aeronaves tilt rotor no tripuladas Los vehículos ‘tilt-rotor’ se caracterizan por el uso de alas y rotores para la generación de la sustentación necesaria a lo largo de las fases de vuelo. Por lo general se establecen tres modos de vuelo principales: vuelo vertical o de helicóptero, transición y vuelo horizontal o de avión [3]. Figura 1-2. Modos de vuelo principales de una aeronave tilt rotor [3]. La ventaja principal de este tipo de vehículos reside en la capacidad de despegar y aterrizar sin necesidad de movimiento horizontal, lo que les permite realizar estas fases en un espacio reducido. Durante estos procedimientos el eje de giro del rotor se encuentra en posición vertical y son las propias palas las que generan el empuje necesario en esta dirección. En la fase de vuelo horizontal el eje de giro se coloca en la dirección del avance proporcionando la velocidad necesaria para que sean las alas fijas las que generen la sustentación necesaria, lo que permite aumentar la eficiencia del vuelo y alcanzar velocidades muy superiores a las de aeronaves como los helicópteros. La concepción de este tipo de aeronaves se originó a principios de la década de 1930 pero no fue hasta los años 70 cuando se demostró su viabilidad con la llegada del Bell XV-15 y más tarde con el diseño y fabricación por parte de Bell y Boeing del V-22 (‘Osprey’). 3 Figura 1-3. XV-15 (izq) y V-22 (dcha). Con la evolución de técnicas de comunicación y control también se comenzaron a desarrollar aeronaves de rotores basculantes no tripuladas. Estos vehículos se enfrentan tanto a los retos propios del control en cada fase del vuelo como a las exigencias estructurales derivadas de la incorporación de esta tecnología. En el área de control, además de las inestibalididades y peculiaridades de las fases de despegue vertical y avance horizontal, se introduce una fase de transición con sus propias dinámicas de vuelo que deben ser analizadas y estudiadas para asegurar la viabilidad de este tipo de aeronaves. En este campo se han realizado diversas investigaciones dando lugar a diseños variados, como por ejemplo el modelo TR-60 del Korean Aerospace Research Institute o el vehículo de 4 rotores AKITSU de la agencia japonesa de exploración aeroespacial. Figura 1-4. TR-60 (izq) y AKITSU (dcha). El modelo AKITSU posee una característica diferenciadora del modelo desarrollado en el proyecto EMERGENTIa al no ser únicamente los rotores los que basculan sino toda la superficie alar, lo que lo convierte en un vehículo de tipo ‘tilt wing’. A este efecto podría considerarse una última categoría en desarrollo denominada ‘Rotor wing’ en la que los rotores se mantienen en una posición horizontal y fija relativa a las alas y son éstas las que durante el despegue y el aterrizaje giran para proporcionar el empuje necesario. Un ejemplo de esta tecnología se puede apreciar en la Figura 1-5 con el modelo ROTORwing de DZYNE Technolgies. Modelado 10 Por otra parte, se incluyen en el sistema 3 juegos de palas que permiten crear un rotor de 7, 8 y 9 pulgadas de diámetro (Figura 2-6). Sin embargo, y con objeto de simplificar el modelado, se decide diseñar una pala simétrica (sin torsión a lo largo de su envergadura) de dimensiones simiares a la pala de tamaño intermedio utilizando como referencia una pala simétrica disponible en línea en GrabCad (Figura 2-6). Figura 2-6. Dimensiones de la pala empleada [mm] (izq) y juegos de palas incluidas en el sistema (dcha). Por último, de este programa también se pueden obtener las posiciones de los centros de gravedad, masas 1 y matrices de inercia (aunque estas últimas son nulas debido a la baja masa de los componentes). Estos datos serán necesarios en posteriores fases del análisis. Tabla 2-1. Masa y centros de masa de los elementos según el programa Solid Edge. Elemento Masa [g] CM-X [mm] CM-Y [mm] CM-Z [mm] Barra1 ~ 0 -7.85 2.00 91.10 Barra2 ~ 0 7.85 2.00 91.10 Buje 1 0.00 0.00 85.12 Cabeza 1 0.00 0.00 97.09 Mástil fino 1 0.00 0.00 50.00 Mástil grueso ~ 0 0.00 0.00 71.47 Unión1 2 -0.47 10.82 85.12 Unión2 2 0.47 -10.82 85.12 Pala1 3 -1.36 74.48 85.11 Pala2 3 1.36 -74.48 85.11 Observando la Figura 2-1 y la Figura 2-5 se puede identificar un elemento que no aparece en el modelado del rótor. Como simplificación del modelo y debido a la complejidad que supone el establecimiento de restricciones de contacto en fases posteriores del análisis, no se va a tener en cuenta para el estudio el tornillo que sirve de unión entre la pala y su anclaje. En los siguientes apartados se comentará cómo se ha diseñado está unión para que, aunque no exista este elemento, el conjunto completo se asemeje lo más fielmente posible a la realidad. 1 Se ha supuesto como material para el cálculo de todos los componentes, a excepción de las palas, el Al 1060 con una densidad de 2712 kg/m3. Para las palas se ha supuesto PVC con una densidad de 1385 kg/m3. Esta selección de materiales se ha realizado sólo para calcular la masa. Pueden existir pequeñas discrepancias entre la masa obtenida y la masa real de los componentes debido a la suposición de los materiales y a ciertos errores a la hora de calcular el volumen de zonas determinadas en el modelo 3D. 11 3 CO-SIMULACIÓN ara calcular la vida a fatiga de un elemento en cualquier aplicación es necesario obtener los esfuerzos que soporta y a partir de ellos las tensiones que se producen, responsables del fallo por fatiga. En este apartado se discute el procedimiento seguido para el cálculo de los esfuerzos. Teniendo en cuenta que en este caso están involucradas fuerzas aerodinámicas, se hace necesario disponer de un modelo que genere dichas fuerzas a partir de los parámetros de vuelo (velocidad de giro, ángulo de paso…). A este fin se decide realizar una co-simulación entre los programas Adams View y Matlab cuyo diagrama de flujo simplificado se puede observar en la Figura 3-1. Figura 3-1. Esquema simplificado del flujo de información de la co-simulación. Durante el documento se hace referencia a diversos conceptos relacionados con los parámetros de vuelo. Con la intención de mejorar la comprensión de estos conceptos, la Tabla 3-1 recoje aquellos más empleados junto con su definición. Tabla 3-1. Definición de conceptos relevantes. Concepto Definición Velocidad de giro/angular Velocidad de giro del eje del rotor Velocidad de avance Velocidad relativa del sistema respecto al fluido circundante Ángulo de paso Ángulo de las palas del rótor (Véase Figura 2-2) T (tracción) Empuje producido por el sistema Ciclo de vuelo Conjunto de fases de un vuelo desde el despegue hasta el aterrizaje Ciclo a fatiga Ciclo de carga que influye en la vida a fatiga del elemento P Co-simulación 12 3.1 Adams View Este programa dispone del modelo 3D del conjunto completo y en él se introducen la gran mayoría de parámetros de vuelo. Estos datos son después trasladados al programa Matlab para el cálculo de las fuerzas aerodinámicas y los resultados obtenidos son devueltos de nuevo al programa Adams donde se calculan los esfuerzos en las diversas partes del sistema. 3.1.1 Importación del conjunto La extensión (formato) por defecto de los elementos 3D diseñados en Solid Edge es la propia de este programa (’.par’), sin embargo, este formato no es compatible con el software del programa Adams, por lo que es necesario convertir cada componente en un elemento ‘parasolid’ para su correcta importación. Existen otros muchos formatos compatibles y no sería extraño encontrar casos en los que la extensión adecuada para esta importación fuera otra. En este caso se ha probado con distintos formatos siendo el comentado el que mejor mantiene las propiedades dimensionales de cada elemento. La importación de los elementos se realiza de manera individual a partir de la herramienta ‘Import’ del programa indicando el fichero correspondiente (Figura 3-2). Es recomendable crear para cada pieza un elemento denominado ‘Part’ ya que simplifica los procesos posteriores. Dentro de la ventana de importación se puede indicar el nombre que se asignará a cada elemento, así como su localización y orientación, sin embargo, la posición por defecto es aquella en la que se había guardado cada fichero, que en este caso es la que corresponde a su posición en el conjunto por lo que no es necesario modificar este apartado. Figura 3-2.Ventana de importación en Adams View. Una vez se han importado todos los elementos, se definen en el programa las masas y los centros de masa. Utilizando la opción ‘Marker’ dentro del módulo ‘Construction’ (Figura 3-3) e introduciendo las cooordenadas calculadas en el apartado anterior (Tabla 2-1), se genera un punto en el espacio que corresponderá al centro de masa de la pieza. Completada esta operación para cada elemento, se modifican las propiedades introduciendo su masa y asignando como centro de masa el ‘Marker’ creado (Figura 3-4). Figura 3-3. Herramienta ‘Marker’ (resaltada en amarillo) empleada para la definición de los centros de masa de los elementos. 13 Figura 3-4. Ventana de modificación de elemento. 3.1.2 Restricciones Llegados a este punto cada pieza está completamente definida y colocada en su posición, pero como un elemento aislado. Para vincular las distintas partes del sistema es necesario especificar las relaciones que existen entre ellas. Las restricciones se establecen mediante la herramienta ‘Joints’ que ofrece 9 tipos diferentes de relaciones entre elementos (siguiendo el orden de la Figura 3-5, de izquierda a derecha y de arriba abajo): ‘Fixed Joint’ (unión fija), ‘Revolute Joint’ (restricción de revolución), ‘Translational Joint’ (par prismático), ‘Cylindrical Joint’ (par cilíndrico), ‘Spherical Joint’ (par esférico), Constant Velocity Joint’ (restricción de velocidad constante), ‘Hooke Joint’ (junta cardan), ‘Screw Joint’ (restricción de rosca) y ‘Planar Joint’ (par plano). Para este caso se emplearán exclusivamente 3: unión fija, par de revolución y par cilíndrico. Para vincular dos cuerpos es necesario seleccionar el tipo de restricción, los elementos involucrados y el punto en el que se colocará dicha restricción. Al igual que en el caso anterior, este punto se genera mediante la opción ‘Marker’. Un ejemplo de estas restricciones es la que ocurre entre la pala y la pieza ‘Unión’. Como se ha comentado anteriormente, no se ha incorporado al modelo el tornillo que fija ambas partes, pero se ha establecido una unión fija entre ellas mediante la herramienta ‘Fixed Joint’ en un punto que se encuentra en el plano medio de la pala y en el centro del cilindro ficticio (tornillo) que uniría ambas partes (Figura 3-6). Aunque se podría suponer que una unión en un punto no es representativa de la vinculación provocada por la colocación de un tornillo, se explicará más adelante cómo se interpretan las fuerzas resultantes en este punto para que la pérdida de información sea mínima. Figura 3-5. Herramienta ‘Joints’. Co-simulación 14 Figura 3-6. Restricción fija entre la pala y la pieza ‘Unión’. 3.1.3 Mediciones y movimientos Dentro de la co-simulación el programa creado será interpretado por el modelo de fuerzas aerodinámicas de Matlab como un modelo real del que se extraen las mediciones necesarias de los parámetros de vuelo, por lo tanto, para transmitir esta información desde el modelo es necesario generar y ‘medir’ las variables involucradas. Los movimientos se aplican a las restricciones creadas en el apartado anterior mediante el módulo ‘Joint Motions’ dentro del apartado ‘Motions’. Para cada movimiento es necesario definir la unión en la que están aplicadas, el tipo y la función que definirá el valor de la variable a lo largo del tiempo. En este caso se han empleado dos movimientos (Figura 3-7): el primero tiene lugar en la restricción de revolución creada en el mástil grueso y es la encargada de simular la velocidad de giro del sistema; el segundo movimiento es el encargado de variar el ángulo de paso de las palas. Para representar la situación real, este último desplazamiento se aplica en el mástil fino en dirección vertical y, como se ha explicado anteriormente, a través de las barras de control se transmite hasta las palas, lo que hace variar su ángulo de paso. Figura 3-7. Movimientos aplicados en el sistema. Las mediciones se realizan sobre las propias restricciones y con ellas se obtiene tanto el valor de los 15 desplazamientos y velocidades de giro generadas a partir de los movimientos, como los esfuerzos resultantes en la unión entre la pala y su anclaje (pieza ‘Unión’), que serán necesarios para el cálculo de las tensiones en este elemento. Como ambas palas están sometidas a los mismos esfuerzos, la medición de éstos se realiza en un único anclaje. Para generar estas mediciones se emplea el módulo ‘Measures’ dentro del apartado ‘Design Exploration’ y se selecciona el elemento correspondiente. Un ejemplo de estas mediciones es la realizada sobre la restricción fija comentada previamente (Figura 3-6) sobre la que se miden los esfuerzos generados, entre ellos la fuerza en dirección vertical (FZ) cuya configuración se espone en la Figura 3-8, donde la restricción denominada ‘JOINT 15’ representa la unión fija aplicada entre ambos elementos. Figura 3-8. Configuración de FZ en la restricción fija existente entre la pala y su anclaje. Una vez terminado este proceso se puede comprobar que las restricciones, movimientos y mediciones sean los deseados realizando una simulación a través del módulo ‘Simulate’ dentro del apartado ‘Simulation’ (Figura 3-9) y observando posteriormente los resultados en el módulo ‘Postprocessor’ del apartado ‘Results’. Figura 3-9. Ventana de configuración de la simulación. 3.1.4 Fuerzas Las fuerzas serán calculadas durante la simulación en el programa Matlab, obteniéndose como resultado los valores de las fuerzas y momentos aplicados de manera continua a lo largo de la superficie de las palas. Sin embargo, el programa Adams no permite la aplicación de fuerzas distribuidas a lo largo de la pala, por lo que es necesario establecer una simplificación. Ésta consiste en aplicar todas las fuerzas sobre una línea ficticia colocada a una distancia del 25% de la cuerda de la pala de manera discreta. El número de puntos de aplicación será explicado y comentado con detalle con posterioridad. Co-simulación 16 La colocación de estas fuerzas en el programa se realiza mediante la definición tanto de los puntos de aplicación como de la dirección de estos parámetros. Para ello se procede a efectuar una división equidistante de la línea ficticia comentada a partir de la que se calculan las coordenadas para ambas palas de los puntos de aplicación. En el programa se introduce un ‘Marker’ por cada punto resultante y se coloca sobre él un vector de fuerzas de 6 coordenadas mediante la herramienta ‘Six component Force’ del apartado ‘Forces’ (Figura 3-10). Los esfuerzos aplicados consisten en una fuerza en dirección vertical y un momento en el sentido contrario al del giro del sistema, lo que conlleva que el resto de fuerzas y momentos tendrán un valor nulo. Para introducir las fuerzas obtenidas se emplea el comando ‘VARVAL’ que extrae la magnitud del valor introducido desde el modelo de fuerzas aerodinámicas en Matlab. Este procedimiento se realiza mediante el empleo de variables denominadas ‘State variables’, cuyo significado se explica en el siguiente apartado. Figura 3-10. Herramienta ‘Six component force’. 3.1.5 Variables de comunicación y exportación Todas las variables de medición, movimientos y fuerzas descritas hasta ahora tienen una identidad local, es decir, sólo pueden ser consultadas desde el propio programa Adams. Sin embargo, es necesario establecer una comunicación entre los programas implicados en la co-simulación. Para ello se crean unas variables de comunicación denominadas ‘State variables’ que permiten el tránsito de la información entre Adams y Matlab. Estas variables se crean dentro del módulo ‘System Elements’ y pueden ser de dos tipos en función de la dirección de la comunicación. Los valores de las fuerzas aplicadas son un ‘input’ del programa mientras que las mediciones de las reacciones y los movimientos son un ‘output’. Para introducir una variable ‘input’ sólo es necesario asignarle un nombre y colocar un ‘0’ en el apartado de valor a lo largo del tiempo (‘F(time)’). Para las variables ‘output’ se sigue el mismo procedimiento, pero en este caso se coloca en el apartado ‘F(time)’ la variable procedente de Matlab de la que extrae el valor. Para ambos tipos es necesario definir la variable como ‘Run-Time Expression’. La ventana de configuración de estas variables se puede observar en la Figura 3-11. Figura 3-11. Variable ‘input’ (izqda.) y variable ‘output’ (dcha). 17 Una vez completados todos los procesos es necesario crear el módulo que será empleado durante la cosimulación. Para ello se utiliza la herramienta ‘Plant Export’ dentro del apartado ‘Plugins’ y se indican las variables ‘input’ y ‘output’ definidas. También se puede seleccionar el tipo de análisis que se va a llevar a cabo y el ‘solver’ a emplear. En este caso se ha escogido ‘FORTRAN’ y el análisis no lineal, que está seleccionado por defecto. 3.2 Modelo de fuerzas (Matlab) Este programa se encarga de generar las fuerzas y momentos en las palas a partir de la velocidad de giro, el ángulo de paso de las palas y la velocidad axial (de avance) del sistema, empleando un modelo de fuerzas basado en la Teoría de Cantidad de Movimiento (TCM) y en la Teoría del Elemento de Pala (TEP) proporcionado por el Departamento de Ingeniería Aeroespacial de la Escuela. En la TEP se asume cada sección de pala como un perfil bidimensional aislado en el que se calcula, basándose en la TCM, el coeficiente de tracción (dCT) y el coeficiente de par (dCQ) entre otros parámetros como el coeficiente de potencia (dCP). La integración de estos coeficientes en las divisiones de la pala comentadas anteriormente da como resultado la tracción (T) y el par (Q) de cada punto de aplicación considerado en la cosimulación. Las expresiones (3-1), (3-2) y (3-3) muestran algunas de las relaciones entre los diferentes parámetros empleadas para el cálculo de los esfuerzos, donde ω representa la velocidad de giro, R es la mitad de la envergadura total del rotor, Cl es el coeficiente de sustentación, Cd representa el coeficiente de resistencia del perfil, ϕ es el ángulo de entrada de corriente, σ representa la solidez de la pala, x es la distancia adimensionalizada desde el eje de giro hasta el punto donde se calculan los esfuerzos y ρ y A son la densidad del fluido circundante y el área barrida por las palas durante el giro respectivamente. 𝑑𝐶𝑇=𝑑𝑇 𝜌𝐴(𝜔𝑅)2=1 2𝜎𝐶𝑙𝑥2𝑑𝑥 (3-1) 𝑑𝐶𝑄= 𝑑𝐶𝑝=1 2𝜎(∅𝐶𝑙+𝐶𝑑)𝑥3𝑑𝑥 (3-2) 𝐶𝑄= 𝑄 𝜌𝐴(𝜔𝑅)2𝑅 (3-3) A efectos de este trabajo y, aunque se han expuesto parte de las expresiones básicas utilizadas, el modelo de fuerzas será considerado como una caja negra, ya que los cálculos internos y variables contempladas en él superan el alcance de este proyecto. El código empleado durante la co-simulación puede consultarse en el apartado C.1 del Anexo C. Por último, es necesario introducir en el programa de forma manual ciertas variables como el peso, el número de palas y la velocidad axial y crear las variables que serán importadas desde y exportadas hacia el programa Adams. Como resultado de la simulación se pueden obtener un gran número de parámetros, aunque los únicos relevantes para este análisis son las fuerzas y momentos comentados. 3.3 Co-simulación Tras generar el modelo en Adams y preparar el modelo de fuerzas se realiza la co-simulación. Para ello se emplea el módulo Simulink del programa Matlab. Tanto el modelo mecánico (Adams) como el modelo de fuerzas (Matlab) conformarán dos bloques interconectados dentro de este programa. Antes de comenzar con el procedimiento de simulación es necesario cargar los archivos generados Co-simulación 18 previamente en el programa Adams seleccionando el archivo con la extensión ‘.m’ (Figura 3-12). Figura 3-12. Archivos generdos en la exportación de Adams. Una vez ejecutado se inicializan las variables del sistema. Para generar el bloque de Simulink se emplea el comando ‘adams.sys’ con el que se abre una ventana que contiene, entre otros, el bloque que se va a utilizar en este análisis, denominado ‘adams_sub’ (Figura 3-13). Figura 3-13. Bloque mecánico en Simulink. Existen diferentes opciones dentro del bloque mecánico. En este caso, se opta por realizar una simulación en modo ‘batch’ (más rápido que ‘interactive’) y se establece el intervalo de comunicación entre los bloques (‘step’) en 0.02 segundos. Sin embargo, la simulación no se reliza en esta ventana por lo que el bloque mencionado se copia en una ventana nueva donde se añade un bloque de la biblioteca de Simulink denominado ‘Matlab function’, en el que se incorpora el modelo de fuerzas. Si todo el proceso se ha realizado de manera correcta se dispondría en este momento de dos bloques aislados. Para conectarlos se debe unir cada salida con su correspondiente entrada (si la tuviera). Todas las variables de salida del modelo de fuerzas son entradas del modelo mecánico y están ordenadas en cada bloque de la misma forma por lo que su conexión es relativamente sencilla. Por otra parte, todas las salidas del modelo mecánico están conectadas a bloques de tipo ‘Scope’ para representar y almacenar los datos de salida de la simulación empleando la opción ‘To Workspace’. Además, las variables de ángulo de paso y velocidad de giro que son salidas de este sistema deben ser conectadas al modelo de fuerzas. La disposición resultante en Simulink se muestra en la Figura 3-14. Por último, antes de ejecutar la simulación, en la opción de configuración de la ventana de Simulink se establecen todos los parámetros por defecto a excepción del tiempo de simulación, que se definirá en función 19 de las necesidades de cada co-simulación desarrollada a lo largo del estudio. Figura 3-14. Modelo de fuerzas (izqda.) y modelo mecánico (dcha) 2 . 3.4 Resultados Como resultado de la co-simulación se obtienen las tres componentes de las fuerzas y momentos en el punto de anclaje entre la pala y la pieza denominada ‘Unión’. Estos resultados variarían en función de las condiciones analizadas pero, debido la discretización de las fuerzas realizada y para establecer la fiabilidad de los datos calculados mediante este procedimiento, se deciden realizar diferentes simulaciones para distinto número de puntos de aplicación de esfuerzos por pala. Se simula desde 1 punto de aplicación por pala, situado en la línea de 25 % comentada y a un 75 % de la envergadura de la pala (Figura 3-15), hasta 33 puntos, pasando por 9 y 17 (Figura 3-16) puntos de aplicación (equidistantes). De esta manera se obtienen los resultados expuestos en la Figura 3-17, donde se comparan los esfuerzos resultantes en función del número de puntos de aplicación de fuerzas por pala para una variación lineal del ángulo de paso (ϑ) manteniendo una velocidad de giro constante de 4000 RPM y una velocidad de avance nula (vuelo a punto fijo). Figura 3-15. Punto de aplicación para las condiciones de 1 fuerza por pala. 2 Aunque se ha comentado que todas las salidas del modelo de fuerzas son entradas del modelo mecánico, también se ha supuesto como salida la velocidad axial del sistema para poder ser registrada. 27 5 ESTUDIOS PREVIOS omo si se tratara de pruebas en un laboratorio, el análisis teórico de un procedimiento como el descrito en este estudio también requiere de pruebas e iteraciones para perfeccionar los modelos empleados. En este capítulo se describen los análisis y pruebas realizadas hasta alcanzar el procedimiento final para las simulaciones. Aunque los resultados de algunas de estas pruebas no conduzcan a una mejora sustancial del análisis, son relevantes a la hora de comprender y evaluar los métodos utilizados. 5.1 Modelo de elementos finitos Los programas de elementos finitos proporcionan una amplia variedad de herramientas que permiten la realización de un análisis. En este apartado se exponen las empleadas dentro de la ventana ‘Mechanical’ del programa Ansys y se discute su idoneidad para el estudio. 5.1.1 Remote Point Como se comenta en el Capítulo 3, los esfuerzos obtenidos en la zona de la cogida de pala son puntuales. Esta suposición simplifica la co-simulación pero no refleja el escenario real. Para extrapolar estas cargas puntuales a las diferentes zonas de aplicación de la pieza se emplea lo que en Ansys se denomina como ‘Remote Point’. Esta herramienta, que se introduce en la sección ‘Model’, permite la definición en el espacio del punto donde son calculados los esfuerzos en la co-simulación y las caras de aplicación de éstos. Como consecuencia, las fuerzas se aplican de manera proporcional (en función de la superficie) en todas las zonas y además se añaden los momentos generados por estas cargas como consecuencia de la distancia entre las caras y el punto de aplicación de las fuerzas. En este estadio del análisis esta herramienta supone una ventaja, ya que permite automatizar un procedimiento que de otra forma se debería realizar de manera manual. Para este caso se seleccionan como caras de aplicación del ‘Remote Point’ tanto los orificios cilíndricos destinados a alojar el tornillo de unión entre ambas partes como las caras en contacto con las palas, por lo que en total las cargas son aplicadas en 6 caras. Éstas están expuestas señaladas en rojo en la Figura 5-1. La introducción de los valores de las cargas se realiza mediante la herramienta ‘Remote Force’ en el grupo ‘Transient’, donde es necesario indicar el ‘Remote Point’ creado con anterioridad. Esta configuración se emplea en las primeras pruebas del análisis, aunque será objeto de modificaciones a lo largo de este. Figura 5-1.’Remote Point’ (círculo amarillo) y caras de aplicación (en rojo). C Estudios Previos 28 5.1.2 Mallado y Convergencia de la malla La base de un análisis con elementos finitos es la malla (división en elementos) empleada para realizar el estudio de la pieza. La elección tanto del número de divisiones como del tipo de mallado (hexaédrico o tetraédrico) puede suponer una variación considerable en los resultados. Las limitaciones numéricas de la versión académica del programa implican que el número máximo de elementos o nodos no puede superar los 32000. Aunque la cifra puede parecer elevada, esta limitación podría influir en la convergencia de la malla, es decir, puede que, aunque el número de elementos/nodos alcance esta cifra, los resultados varíen y sean más precisos para un número de elementos/nodos mayor. Para intentar solventar este problema, y teniendo en cuenta las recomendaciones del propio programa, se decide emplear un mallado hexaédrico que proporciona una mayor eficacia y necesita de un número menor de elementos/nodos para realizar un análisis fiable. ‘Hex meshing provides computational efficiency where less nodes and elements are required to achieve high solution accuracy.’- ANSYS e-learning Este método es implementado en el programa Ansys mediante la herramienta ‘Method’ de la sección ‘Mesh’. Es necesario especificar la geometría a la que se aplicará la malla, que en este caso es toda la pieza, y el tipo de método a emplear. En el apartado ‘Mapped Mesh Type’ se selecciona ‘Hexa’, que es el método de mallado hexaédrico. Por otra parte, en el apartado ‘Free Mesh Type’ se debe seleccionar el tipo de mallado que se desea lograr en las zonas en las que no es posible el mallado hexaédrico. Esto puede suceder en zonas en las que una división en elementos hexaédricos no se adapte de manera correcta a la geometría de la pieza. En este caso se ha seleccionado ‘Hexa dominant’ que, aunque introduce tetraedros en las zonas conflictivas, hace prevalecer el uso de hexaedros. Este hecho se puede apreciar en la Figura 5-2, donde se representa una sección transversal de la pieza mallada en la que se observa como ciertas zonas del interior del volumen han sido malladas con tetraedros para poder completar esta operación. Figura 5-2. Sección de la pieza mallada utilizando la función ‘Hexa Dominant’. Otro motivo de fallo es el empleo de elementos de gran tamaño, lo que dificulta la adaptación a la geometría. Por ello es recomendable el empleo de la herramienta ‘Sizing’ para disminuir/aumentar el tamaño de los elementos en las zonas donde este fenómeno sea apreciable (Figura 5-3). Además, también se emplea esta utilidad para realizar la convergencia de malla, lo que permite disminuir el tamaño de los elementos en las zonas donde la exigencia es mayor y aumentar el tamaño en las partes de menor importancia para el análisis. En este caso a estudio, y como cabría esperear, la zona donde la magnitud de las tensiones es mayor es aquella que se encuentra en contacto con la pala y donde se aplican las cargas, por lo que es ésta la que centrará el estudio de convergencia. 29 Figura 5-3. Mallado hexaédrico y resultado de la aplicación de la herramienta ‘Sizing’ en la zona sombreada. Como se ha comentado, se simulan distintos modelos de mallado modificando el número de elementos en cada caso, pero aplicando las mismas cargas que corresponden a las de un vuelo con velocidad angular constante y variación del ángulo de paso a lo largo del tiempo. Antes de entrar a valorar los resultados es necesario indicar el resto de condiciones de contorno que se han aplicado. En primer lugar, se han establecido restricciones cilíndricas en los orificios de ancleje al rotor y a la barra de control, señaladas en la Figura 5-4 (flechas amarillas), dejando libre el grado de libertad que permite el giro en ambos mediante la herramienta ‘Cylindrical Support’ en el apartado ‘Transient’ (si sólo existiera uno de los orificios estas condiciones en desplazamientos permitirían el movimiento del elemento como sólido rígido, sin embargo, las restricciones en ambas caras impiden esta circunstancia). Por otra parte, se ha comentado en la sección anterior el modo de aplicación de las cargas centrándose en las fuerzas. Para introducir los momentos se emplea el mismo procedimiento, pero se usa la opción ‘Remote Moment’ en el módulo ‘Transient’. Por útlimo, en el apartado ‘Analysis Settings’ dentro de ‘Transient’ se indica el intervalo de la simulación y el intervalo de tiempo entre los cálculos que se ha fijado en 0.02 s. Las características geométricas del elemento ‘unión’ hacen inadecuada su implementación en el análisis de manera directa debido a la existencia de una arista a 90º que generaría un concentrador de tensiones infinito (aunque la solución proporcionada por el programa fuera discreta). Para solucionar este problema se aplica a esta zona un radio de acuerdo de 1 mm (señalada en la Figura 5-4 con una flecha roja). Figura 5-4. Caras de aplicación de las restricciones cilíndricas (flechas amarillas) y redondeo de 1 mm sobre las aristas (flechas rojas). El resultado del estudio de convergencia de la malla se expone en la Figura 5-5, donde se indica el valor de la tensión máxima para diferentes divisiones en elementos utilizando la misma velocidad de giro y avance que en el capítulo anterior (4000 RPM y vuelo a punto fijo). Estudios Previos 30 Figura 5-5. Variación de las tensiones máximas en función del número de nodos empleados en la malla hexaédrica. Se observa una gran disparidad en los resultados, lo que puede indicar que el método utilizado para realizar la malla no es correcto. La naturaleza casi cilíndrica de la pieza en algunas zonas hace pensar que la adaptación de los elementos hexaédricos no es la correcta y provoca estas diferencias. Por otra parte, y como se comentará posteriormente, en el análisis a fatiga es necesario establecer aquellos puntos críticos en los que la exigencia (tensiones) es mayor y empleando estos elementos se produce una gran variación en la posición y valor de la tensión de estos puntos. Por lo tanto, aunque los métodos de mallado con hexaedros suelen proporcionar mejores resultados en geometrías más simples, la geometría empleada y la limitación en el número de elementos dificultan su utilización en este caso. La única opción restante es el empleo de elementos tetraédricos. Estos elementos se adaptan mejor a la geometría y con ellos los resultados tienden a la convergencia (Figura 5-6). Estos resultados son suficientes para considerar que su uso es adecuado en el análisis, aunque sería necesario en cualquier caso un estudio en mayor profundidad sobre el comportamiento de la malla que no se puede realizar en este caso debido a las limitaciones numéricas comentadas. Figura 5-6. Variación de las tensiones máximas en función del número de nodos empleados en la malla tetraédrica. La malla empleada en el resto del análisis está formada por tanto por 20024 elementos tetraédricos y 31620 nodos (Figura 5-7). 31 Figura 5-7. Malla tetraédrica empleada en el estudio. 5.1.3 Resultados Para la obtención de los resultados de la simulación se emplea el módulo ‘Solution’ que se emplea de manera similar en Ansys WorkbenchK al módulo ‘Postprocessor’ en el tradicional Ansys APDL. Dentro de este módulo se puede escoger entre una gran variedad de resultados (tensiones, deformciones, energía…) siendo los empleados en este proyecto los correspondientes a las tensiones normales (incluyendo las principales) y las tensiones tangenciales. Además, para que las magnitudes resultantes de las tensiones no sean calculadas en los nodos sino en los elementos, se sustituye en cada resultado la opción ‘Averaged’ por ‘Elemental Mean’. Por último, otros resultados de interés necesarios en posteriores fases del análisis son los denominados como ‘Vector Principal Stress’, que permite la obtención de los ángulos de Euler, y ‘User Defined Result’, con el que se puede operar con distintos resultados (suma, resta, multiplicación…) para conseguir aquel que se ajuste de manera más adecuada a la solución requerida. Figura 5-8. Módulo ‘Solution’ en Ansys Workbench. En el Anexo D se expone una descripción más detallada de la configuración realizada en Ansys para la realización del análisis. 5.2 Influencia de los parámetros de vuelo (velocidad de giro y velocidad axial) Hasta ahora se ha empleado una velocidad de giro de 4000 RPM y una velocidad axial nula por ser estos parámetros habituales y, aunque conservadores, apropiados para la realización de las pruebas anteriores. Sin embargo, para poder sacar el máximo partido a la aeronave se requiere estudiar tanto la velocidad máxima a la que puede girar el rótor como la velocidad de avance máxima para cada caso. Los resultados obtenidos no Estudios Previos 32 serán empleados en todas las fases del análisis, pero con esta sección se pretende describir el procedimiento de cálculo de estos en función de los diferentes parámetros. 5.2.1 Velocidad de giro La velocidad de giro máxima está limitada (suponiendo que la potencia del motor es suficiente) por la velocidad lineal en la punta de las palas. Si esta velocidad se acerca a la velocidad del sonido se pueden producir fenómenos indeseados que reducen la eficiencia del sistema. Por lo tanto, se calcula la velocidad de giro máxima como aquella para la que la velocidad en punta de pala no supere 0.6 Mach. Considerando 1 atmósfera y 15 °C y sabiendo que el radio de giro del sistema a estudio es de 0.136 m esta velocidad lineal en punta de pala equivale a 14340 RPM. Las tensiones producidas en el sistema con esta velocidad de giro son muy superiores a las calculadas anteriormente, ratificando así lo comentado en los apartados previos; la fuerza que ejerce la pala sobre el elemento a estudio como consecuencia de su giro es la de mayor valor en el estudio y la que más influencia tiene sobre la magnitud de las tensiones. 5.2.2 Velocidad de avance (Vc) Una velocidad de avance nula (vuelo a punto fijo) no refleja de manera adecuada las condiciones que tienen lugar en la aeronave a estudio, donde durante la mayor parte del vuelo la velocidad axial será muy superior a 0. Esta velocidad depende tanto de la velocidad de giro del rotor como del ángulo de paso de las palas en cada instante y su máximo se da cuando la potencia necesaria para mover el sistema es negativa. Sin embargo, con el modelo aerodinámico que se emplea en este estudio, antes de alcanzar una potencia negativa se obtiene un empuje negativo (autorrotación). Este fenómeno se produce cuando es el propio aire el que impulsa el giro de las palas debido a su velocidad y será considerado el límite para el cálculo de la velocidad máxima. En la Figura 5-9 se representa la velocidad axial en función del ángulo de paso para la que se consigue un empuje negativo empleando la velocidad de giro máxima calculada en el apartado anterior. Figura 5-9. Variación de la velocidad de avance máxima con el ángulo de paso (en grados). Si se calculan las tensiones utilizando estas velocidades de avances se obtiene un resultado que podría considerarse poco intuitivo. Cuánto mayor es la velocidad de avance, menores son las tensiones producidas en 33 el elemento a estudio ya que el empuje se reduce hasta alcalzar un empuje negativo para velocidades superiores. Este efecto está generado por el fenómeno de autorrotación que será comentado más adelante. 35 6 FATIGA l grueso del análisis realizado en este estudio se centra en la investigación del método más adecuado y su aplicación para la determinación de la vida a fatiga del elemento expuesto. En este capítulo se discuten las distintas características que definen el problema y las diferentes aproximaciones a este mediante los métodos de fatiga multiaxial considerados. 6.1 Contexto Los esfuerzos producidos durante el vuelo de una aeronave de este tipo se producen en todas las direcciones y a priori no cuentan con un patron de aplicación definido. Por lo tanto, es coherente pensar que los métodos más apropiados son aquellos que permiten el análisis a fatiga multiaxial no proporcional, es decir, que son capaces de determinar la vida a fatiga ante cargas que no están en fase y que provocan tensiones en más de una dirección. Actualmente existen varios modelos de este tipo que son ampliamente utilizados en distintos ámbitos. La mayoría de los descritos en este análisis se basan en el estudio del plano crítico, es decir, no se centran en una única dirección como los modelos uniaxiales, sino que, a partir de las tensiones (principales, normales y tangenciales) determinan la dirección más desfavorable empleando un parámetro de daño que se calcula para todas las direcciones y que informa sobre la severidad (entendida como la variable que rige el número de ciclos hasta el fallo) en cada una de ellas. Esta dirección a su vez es calculada en un punto crítico. Otra forma de analizar el problema sería mediante la simplificación de las tensiones (empleo de tensiones equivalentes) y el uso de modelos uniaxiales. En algunos casos estos métodos más sencillos pueden ofrecer resultados fiables simplificando el procedimiento, por lo que también serán objeto de estudio y comparación con los modelos multiaxiales en este capítulo. Por último, en la resistencia a fatiga influyen una gran variedad de parámetros tales como la temperatura, el acabado superficial, el método de fabricación, etc…Sin embargo, debido a que se desconocen en el caso a estudio y a que el objetivo principal es encontrar un método adecuado, su influencia no se va a contemplar, aunque es importante destacar que en análisis posteriores deberían ser tenidos en cuenta. A continuación, se exponen los diferentes modelos de fatiga estudiados para la determinación de la vida del componente, analizando las ventajas e inconvenientes de cada uno de ellos y su grado de utilidad para este estudio. 6.2 Modelo de Findley Este criterio de plano crítico emplea un parámetro de daño ‘f’ que está relacionado con la resistencia a fatiga y compuesto por tensiones tangenciales, responsables del inicio de la grieta, y normales, que influyen en la propagación de ésta, y por lo tanto se trata de un modelo basado en tensiones. Estos criterios se diferencian de los basados en deformaciones en que predicen de forma más precisa la vida a fatiga para los casos en los que la influencia de la deformación plástica no es significativa [8], como es posible que ocurra en este caso. El parámetro de daño se define en la siguiente ecuación donde ∆τ representa la variación de tensiones tangenciales en el plano, 𝜎 es la tensión normal al plano y k es una constante del material relacionada con la sensibilidad del material a las tensiones normales. El valor de esta constante (0.3 para el aluminio) y la definición del criterio se han obtenido de la herramienta virtual eFatigue. ∆τ 2+𝑘𝜎=𝑓 (6-1) E 43 7 APLICACIÓN AL MODELO iseñado por completo el procedimiento de obtención de tensiones y seleccionado el criterio de fatiga a emplear, la siguiente fase del análisis consiste en aplicar al modelo las cargas características de un ciclo de vuelo completo y obtener el número de estos ciclos que ocurren antes de que se produzca el fallo. En este capítulo también se estudian y ajustan algunos de los parámetros y procedimientos que influyen en la obtención de la vida a fatiga. 7.1 Punto crítico Previamente se ha definido el punto crítico como aquel para la que la tensión de von Mises es máxima, sin embargo, el signo de esta tensión es siempre positivo y puede que se esté escogiendo como punto crítico un punto en el que las tensiones sean altas pero de compresión, y por consiguiente es probable que se pierda información muy importante al ser las tensiones de tracción las que provocan realmente el daño a fatiga. Otra opción consiste en utilizar como parámetro para determinar el punto crítico la tensión máxima principal. Aunque es una opción que se podría considerar factible, del criterio a fatiga se extrae que no es sólo este parámetro el que influye en el daño, sino que también hay que tener en cuenta las tensiones/deformaciones tangenciales. Por lo tanto, el parámetro utilizado para hallar el punto crítico debe reflejar, en la medida de lo posible, el daño de manera similar al método de Fatemi-Socie. La opción propuesta consiste en emplear no un parámetro aislado, sino una ecuación de la forma: 𝑃𝑎𝑟á𝑚𝑒𝑡𝑟𝑜 𝑐𝑟í𝑡𝑖𝑐𝑜=(𝑆1−𝑆3)×(𝑆1+𝑆3) (7-1) Donde S1 es la tensión principal máxima y S3 la tensión principal mínima. Este parámetro refleja, de manera proporcional, aquel punto cuyo plano crítico forma el círculo de Mohr para las teniones 1 y 3 crea el rectángulo de mayor tamaño en el cuadrante superior derecho (para que el parámetro empleado fuera el áera encerrada por este rectángulo sería necesario dividir por 4). Una representación gráfica de lo comentado se espone en la Figura 7-1. Figura 7-1. Representación gráfica del parámetro de daño escogido. Para implementar este parámetro en el programa Ansys se utiliza la opción ‘User Defined Result’ dentro del apartado ‘Solution’ y se introduce la ecuación comentada cuya notación en el programa coincide con la expuesta. D Aplicación al Modelo 44 7.2 Análisis de la aplicación de esfuerzos En las fases iniciales del modelado de las zonas de aplicación se ha empleado únicamente la herramienta ‘Remote Point’ para incorporar las fuerzas y momentos sobre las caras seleccionadas. Sin embargo, con este método todas las caras soportan todos los esfuerzos y no se tiene en cuenta, por ejemplo, que las caras en contacto con las palas no resisten realmente las fuerzas en la dirección longitudinal de ésta, entre otras suposiciones. Además, tampoco se ha considerado previamente la fuerza de compresión que ejerce el tornillo de anclaje entre la pala y la pieza a estudio. 7.2.1 Compresión ejercida por el tornillo Para calcular la fuerza ejercida por el apriete del tornillo sobre el sistema se emplea como referencia la norma UNE 17108-1981 (Tornillos y tuercas de acero). En concreto, al tratarse de un tornillo M2 de cabeza cilíndrica (Figura 7-2), se utiliza la tabla II de la norma expuesta, que describe además los parámetros para el material de peor calidad, lo que se asemeja más al tipo de tornillo empleado en esta aplicación. Figura 7-2. Tornillo M2 de sujeción de la pala. De esta tabla se obtiene el momento de apriete (M) en Nm de cada tipo de tornillo, asumiendo un esfuerzo máximo del 80 % sobre el límite elástico, que se emplea para obtener la carga de tracción sobre el tornillo (Ft) en decanewtons a través de la ecuación (7-2), donde P representa el paso de la rosca en milímetros, μ es el coeficiente de rozamiento entre las superficies de contacto, d2 es el diámetro medio de la rosca del tornillo, De es el la distancia entre caras de la cabeza del tornillo o de la tuerca y Di representa el diámetro mínimo de la superficie de contacto de la cabeza del tornillo (diámetro nominal) o de la tuerca. 𝑀= 𝐹𝑡 1000[0.161 𝑃+ 0.583 𝜇 𝑑2+0.25 𝜇 (𝐷𝑒+𝐷𝑖)] (7-2) Como se desconoce el coeficiente de rozamiento entre las superficies, se supone uno intermedio de valor μ=0.1. En la Tabla 7-1 se reflejan los valores de cada variable empleados para el cálculo de la fuerza apriete. 45 Tabla 7-1. Valores empleados para el cálculo del apriete del tornillo. Variable Valor M 0.0082 Nm μ 0.1 P 0.4 mm D 2 mm d2 1.74 mm De (según norma ISO DIN 13) 4 mm Di 2 mm La fuerza que se obtiene como resultado equivale a la tracción que soporta el tornillo una vez aplicado el par de apriete, cuyo valor es de Ft = 259.6 N. Esta fuerza es muy elevada, sobre todo teniendo en cuenta que por simplicidad no se va a diseñar la zona de aplicación de esta fuerza en las roscas interiores de los orificios, sino que las superficies de aplicación se encuentar en la cabeza del tornillo y de la tuerca. Considerando que este método para la obtención de la fuerza está dirigido a aplicaciones industriales y que por consiguiente el resultado puede ser excesivo para esta aplicación, se emplea para el resto del análisis un 10% del valor obtenido. Puede parecer un porcentaje pequeño, pero en la mayoría de encastres de este tipo el grado de libertad de la pala que corresponde al giro relativo sobre el tornillo (arrastre) no se suele restringir de manera completa, debido a que las propias fuerzas de inercia colocan la pala en la posición adecuada. Figura 7-3. Zona de aplicación 4 de la fuerza de apriete del tornillo. 7.2.2 Reestructuración de las zonas de aplicación de los esfuerzos Para aplicar los esfuerzos más coherentemente se dividen las caras usadas previamente con la herramienta ‘Remote Point’ y las caras de aplicación del apriete del tornillo en 4 zonas de aplicación principales: orificios de los tornillos (Figura 7-4a), apoyos tuerca/tornillo (Figura 7-3), caras en contacto con la pala superiores (según el eje Z) (Figura 7-4b) y caras en contacto con la pala inferiores (Figura 7-4c). En la Tabla 7-2 se resumen los esfuerzos aplicados en cada zona. Las fuerzas en la dirección X e Y se aplican exclusivamente en el alojamiento del tornillo y sólo para las caras que se encuentran en la dirección de cada fuerza. Las caras en contacto con la pala superiores e inferiores resisten todos los momentos mediante un ‘Remote Point’ (que ahora sólo comprende estas 4 caras) y, además, sobre ellas están aplicadas las fuerzas que se generan como reacción de la pala al apriete del tornillo. Asismismo, sobre las caras superiores se aplica también la fuerza de sustentación (Fz). 4 En la parte no visible se encuentra la otra zona de aplicación cuya geometría es idéntica a la representada Aplicación al Modelo 46 Tabla 7-2. Resumen de la aplicación de fuerzas. Esfuerzos Orificios Caras superiores Caras inferiores Apoyo tuerca/tornillo FX FY FZ MX MY MZ Aprietedel tornillo Figura 7-4. Orificio de los tornillos (a), caras superiores (b) y caras inferiores (c). 7.3 Otros aspectos a tener en cuenta Debido las altas velocidades de rotación y las consecuentes fuerzas de inercia que se dan durante el funcionamiento normal del sistema se pueden producir fenómenos que afecten al resultado final de las tensiones en la pieza a estudio. En esta sección se comenta la influencia de la no-linealidade geométrica y la rigidización centrífuga sobre el modelo analizado. 7.3.1 No-linealidad geométrica En casos en los que se producen grandes desplazamientos es importante considerar las no-linealidades geométricas. La presencia de grandes desplazamientos no implica necesariamente unas deformaciones elevadas, por ejemplo, en el caso de una viga de una envergadura suficiente y ante las cargas adecuadas puede producirse una flecha grande, aunque las deformaciones no lo sean. Las deformaciones elevadas sólo pueden ocurrir ante gradientes de desplazamiento elevados, aunque, incluso en este caso, se podrían considerar las leyes lineales si las deformaciones son bajas (menores al 2%) [14]. En piezas sometidas a rotación (como es el caso) los efectos de los grandes desplazamientos pueden ser significativos y las simplificaciones lineales pueden dejar de ser correctas. Para comprobar la influencia de estas no linealidades en la pieza a estudio se realiza un estudio similar a los empleados previamente, pero activando las no linealidades geométricas en el modelo de elementos finitos. Esto se lleva a cabo activando ‘Large deflection’ en el apartado ‘Analysis Settings’ del modelo en Ansys. La comparación del resultado con esta opción activada e inactivada no arroja diferencias ni en la posición del punto crítico ni tampoco en los desplazamientos y deformaciones, siendo el valor de estos parámetros de 0.1 mm y 0.3 % respectivamente para 7000 RPM y 11º de ángulo de paso. Aunque se determina que la influencia de este fenómeno es despreciable, la opción ‘large deflection’ se mantiene activada durante las siguientes fases del análisis ya que los parámetros de vuelo utilizados en este análisis preliminar pueden diferir de los parámetros finales. 47 7.3.2 Fuerza centrífuga Aunque en la prueba anterior las fuerzas de inercia de la pala sí se tenían en cuenta (estaban dentro de las fuerzas que se obtenían de la co-simulación) las fuerzas de inercia de la pieza a estudio no se contemplaban. Las fuerzas de inercia debidas a la rotación pueden ser importantes y además pueden provocar fenómenos no lineales como los anteriormente comentados. Para estudiar esta variable, se introduce la fuerza de inercia con el comando ‘Rotational velocity’ y se mantiene activado ‘Large Deflection’. En ‘Rotational velocity’ se introduce el valor de la velocidad de giro y el punto y eje de rotación. Hay que tener en cuenta que las restricciones impuestas en el sistema son suficientes al quedar los desplazamientos axiales (en este caso dirección Y) también restringidos a través de las condiciones establecidas con la herramienta ‘Cylindrical Support’. Aplicando los mismos parámetros de vuelo que en el apartado anterior no se aprecian diferencias en ninguna de las variables analizadas (tensiones, desplazamientos, …). En cualquier caso, y debido a la misma razón expuesta anteriormente, este análisis se repetirá cuando se disponga de los parámetros de vuelo finales. 7.4 Fases de vuelo La definición de estas fases es muy importante para el análisis, ya que rigen la magnitud de los esfuerzos que soporta el sistema y por lo tanto afectarán de manera apreciable a la vida a fatiga del elemento a estudio. El vuelo típico de una aeronave convencional puede dividirse en 3 fases principales: despegue, vuelo de crucero y aterrizaje. Sin embargo, al tratarse de una aeronave con posibilidad de despegue y aterrizaje vertical (VTOL), debe existir una fase de transición entre el despegue y el vuelo de crucero horizontal. Esta fase a su vez se divide para el análisis en dos subfases: la primera es parecida a la fase anterior pero el empuje necesario es menor ya que comienza a producirse una velocidad horizontal que provoca cierta sustentación por parte de las alas. La segunda se asemeja más a una fase de crucero en la que la mayor parte del peso es sustentado por las alas, pero la velocidad horizontal alcanzada todavía no es la de crucero. Por otra parte, la propia fase de crucero se divide en tres fases en función de la velocidad de avance: crucero nominal, fase que que se desarrolla la mayor parte del vuelo donde la velocidad de avance es relativamente alta; crucero lento, en el que se baja la velocidad para las circunstancias que así lo requieran y crucero ante aterrizaje, que se produce como preparación de la fase de aterrizaje y en la que la velocidad de avance se ve aun más disminuida. Por último, la fase de aterrizaje implica una velocidad de avance negativa que el modelo de fuerzas aerodinámicas empleado no puede generar y por tanto esta fase no será analizada en el estudio. Además, hay que tener en cuenta que la velocidad de giro (principal parámetro que aumenta la magnitud de los esfuerzos) en esta fase es reducida y por lo tanto el daño derivado de ella no será tan relevante como en el resto de fases. 7.4.1 Cálculo de los parámetros de vuelo Cada fase de vuelo se caracteriza por 3 parámetros: velocidad de giro del rotor, ángulo de paso de las palas y velocidad axial del sistema. En una primera iteración para definir estas fases se decide emplear una velocidad de giro máxima para la que la velocidad de punta de pala sea el 245 m/s, con lo que se obtiene una velocidad de 17218.14 RPM. Esta velocidad servirá como referencia para el resto de fases en las que se empleará un porcentaje de esta. Este porcentaje y los otros dos parámetros de vuelo se escogen en función de las magnitudes empleadas normalmente para este tipo de ciclos de vuelo y están reflejadas en la Tabla 7-3. Aplicación al Modelo 48 Tabla 7-3. Parámetros de vuelo iniciales. Fase Velocidad de giro Ángulo de paso Velocidad axial Despegue 100 % (1803.08 rad/s) 12º 5 m/s Transición 1ª parte Transición 2ª parte 70 % (1262.156 rad/s) 8º 5 m/s 40 % (721.232 rad/s) 6º 15 m/s Crucero lento Crucero nominal Crucero ante aterrizaje 50 % (901.54 rad/s) 6º 25 m/s 80 % (1442.464 rad/s) 12º 35 m/s 50% (901.54 rad/s) 6º 25 m/s Sin embargo, estos datos se han seleccionado sin tener en cuenta el modelo de fuerzas aerodinámicas empleado y en la mayoría de fases los parámetros de vuelo introducidos originan errores, debido en su mayor parte a la falta de concordancia entre la velocidad de avance y el resto de parámetros. Para solucionar el problema y disponer de unos datos realistas pero que puedan ser incorporados en la simulación, se emplea el propio modelo de fuerzas para obtener los parámetros de vuelo. A diferencia del modelo usado en la co-simulación, este programa utiliza como entradas el peso, la velocidad de giro y la resistencia aerodinámica, obteniéndose como resultado el ángulo de paso que genera el empuje necesario. El código de este modelo se expone en el apartado C.2 del Anexo C. La velocidad de giro se obtiene a partir de las indicaciones de la empresa especializada en aeromodelismo APC [15] según la ecuación (7-3) donde Dr representa el diámetro del rotor completo en pulgadas. 𝑅𝑃𝑀𝑀𝐴𝑋=150000×𝐷𝑟=1.4023 104 RPM (7-3) Durante un ciclo de vuelo habitual la velocidad de giro nunca alcanza este valor, sino que se suele utilizar en torno a un 80 % de ella. Este porcentaje de utilización se denomina ‘posición de palanca’ y se suele expresar en tanto por 1. Por último, será necesario definir, para cada condición de vuelo, la resistencia aerodinámica del sistema (D). El método de cálculo de este parámetro depende en gran medida del tipo de vuelo (vertical/horizontal) y será comentado para cada fase de vuelo. 7.4.1.1 Despegue La aeronave a estudio despega de manera horizontal y cuenta con dos rotores móviles, por lo que el empuje de cada rotor (T) necesario para llevar a cabo esta operación debe ser igual a la mitad del peso (W) más la mitad de la resistencia aerodinámica totales, como queda reflejado en la ecuación (7-4). Para esta fase se considera que la velocidad de ascenso es constante e igual a 5 m/s. 2𝑇=𝑊+𝐷 (7-4) El valor de la resistencia aerodinámica para esta fase se calcula a partir de la presión dinámica (q), la superficie alar en planta (S) y el coeficiente aerodinámico de resistencia (CD). La presión dinámica se calcula mediante la densidad del fluido, en este caso la del aire, y la velocidad de avance del sistema. El resto de los parámetros son datos conocidos ya que para S se ha empleado un modelo 1:2 de la aeronave real (S=0.1225 m2) y para CD se ha supuesto un valor de 0.5. Esta estimación de parámetros es preliminar y definen datos del modelo que quedan fuera del alcance de este estudio. Por último, se ha supuesto una masa total de la aeronave de 2.6 kg. Las expresiones utilizadas para el cálculo de la resistencia aerodinámica se exponen en las ecuaciones (7-5) y (7-6). 𝐷=𝑞𝑆𝐶𝐷 (7-5) 49 𝑞=1 2𝜌𝑣2 (7-6) Estos valores y ecuaciones se introducen en el modelo de fuerzas para calcular el ángulo de paso necesario, obteniéndose un valor de 9.545°. Tabla 7-4. Parámetros de vuelo en la fase de despegue. V. Avance Pos. Palanca RPM Áng. paso T D/2 5 m/s 0.8 1.1218 104 9.545° 13.18 N 0.4307 N 7.4.1.2 Crucero nominal En esta fase las velocidades de avance son más altas pero el empuje necesario es menor, ya que es la propia sustentación de la aeronave, y no el empuje del rotor, la que contrarresta el peso, por lo que éste no se tiene en cuenta para calcular T en la ecuación (7-4). Por otra parte, La resistencia aerodinámica también disminuye al ser CD en este caso de la forma: 𝐶𝐷=𝐶𝐷0+𝐾𝐶𝐿2 (7-7) Donde CD0 es el coeficiente de resistencia parasitaria, CL es el coeficiente de sustentación y K es el coeficiente de resistencia inducida. Los valores empleados de estas variables se exponen en la Tabla 7-5 y, de manera similar a lo comentado con anterioridad, se trata de valores aceptables para el modelo a escala empleado pero su definición escapa al alcance de este trabajo, por lo que, si bien sus valores pueden ser ajustados de manera más precisa, se considerarán estos datos como adecuados para el estudio aunque se traten de simplificaciones. Al igual que en el caso anterior, se introducen estos valores en el modelo de fuerzas y se obtiene que el ángulo de paso necesario para una velocidad de avance de 15 m/s es de 7.85°. Tabla 7-5. Parámetros de vuelo en la fase de crucero nominal. V-Avance Pos.Palanca RPM Áng. paso T D/2 CD0 CL K 15 m/s 0.8 1.1218 104 7.85° 0.2326 N 0.2326 N 0.02 0.5 0.04 Cuando se aumenta la velocidad de avance se produce un fenómeno anómalo a lo largo de la pala. Al haber considerado una pala simétrica (sin torsión), cuando las velocidades de avance aumentan, la tracción que produce la pala en las zonas cercanas al eje es negativa (autorrotación), al ser el propio aire el que hace girar la pala en estas zonas donde la velocidad lineal de la pala es menor. Este efecto se contrarresta en la zona más alejada de pala que genera el empuje necesario tanto para conseguir ese avance como para anular el empuje negativo, siendo la tracción resultante de toda la pala positiva. Esta es la razón por la que en muchos casos se emplean palas torsionadas, ya que se consiguen ángulos de paso mayores en las zonas donde se produciría en caso contrario la autorrotación. Este efecto se aprecia en la Figura 7-5 donde se ilustra el valor del diferencial de empuje a lo largo de la pala cuya envergadura ha sido adimensionalizada a 1. Aplicación al Modelo 50 Figura 7-5. Diferencial de empuje a lo largo de la pala para la fase de crucero nominal. 7.4.1.3 Crucero lento Para obtener los parámetros de esta fase se han usado las mismas expresiones y valores de la resistencia que en el caso anterior. Como ya se ha comentado en capítulos previos, en este tipo de aeronaves es habitual emplear el ángulo de paso para controlar la velocidad de avance manteniendo la velocidad de giro constante. Por ello, en este caso en el que la velocidad axial se ha supuesto de 10 m/s, el ángulo de paso obtenido es menor que en la fase de crucero nominal e igual a 5.134°. Tabla 7-6.Parámetros de vuelo en la fase de crucero lento. V-Avance Pos.Palanca RPM Áng. paso T D/2 CD0 CL K 10 m/s 0.8 1.1218 104 5.134° 0.1034 N 0.1034 N 0.02 0.5 0.04 Al igual que en el caso anterior y, debido a que la velocidad de avance es todavía elevada, se produce el fenómeno de autorrotación. 7.4.1.4 Transiciones Se suponen dos fases de transición: la primera a una velocidad ligeramente superior a la de despegue y con un modelado de la resistencia aerodinámica igual que el de esta fase, pero con un peso de la mitad para simular la disminución del peso provocada por la sustentación de la aeronave; la segunda es parecida a la fase de crucero nominal, con un modelado de la resistencia realizado de manera análoga, pero con un CL superior ya que en esta fase la sustentación no es todavía suficiente como para contrarrestar la misma parte del peso. La velocidad en esta segunda fase será ligeramente inferior a la de la fase de crucero y la posición de palanca para ambas fases se mantiene en 0.8 (80%). Tabla 7-7. Parámetros de vuelo de la 1º parte de la transición. V. Avance Pos. Palanca RPM Áng. paso T D/2 8 m/s 0.8 1.1218 104 8.08° 7.479 N 1.1025 N 51 Tabla 7-8.Parámetros de vuelo de la 2º parte de la transición. V-Avance Pos.Palanca RPM Áng. paso T D/2 CD0 CL K 12 m/s 0.8 1.1218 104 6.6124° 0.893 N 0.893 N 0.02 2 0.04 7.4.1.5 Crucero ante aterrizaje Por último, en preparación para el aterrizaje, la aeronave debe volar a una velocidad horizontal menor, del entorno de 6 m/s. En este caso no es posible disminuir el ángulo de paso manteniendo la velocidad de giro para conseguir esta velocidad, ya que la velocidad de giro es demasiado elevada para estas condiciones de resistencia y velocidad, por lo que se disminuye la posición de palanca y, por lo tanto, el régimen de giro. El método de cálculo de la resistencia aerodinámica es idéntico al empleado en el resto de las fases de crucero. Tabla 7-9.Parámetros de vuelo de la fase crucero ante aterrizaje. V-Avance Pos.Palanca RPM Áng. paso T D/2 CD0 CL K 6 m/s 0.6 8.414 103 4.0688° 0.0372 N 0.0372 N 0.02 0.5 0.04 7.5 Procedimiento de daño acumulado Para analizar la vida a fatiga de cualquier elemento es necesario definir y caracterizar los ciclos de carga que soporta. En los casos más simples, con una tensión uniaxial cuya variación es constante e idéntica a lo largo del tiempo, sólo es necesario calcular número de ciclos hasta el fallo para averiguar la vida a fatiga. Sin embargo, en el caso a estudio estos ciclos no son iguales ya que, por ejemplo, las cargas soportadas en cada fase difirien. En consecuencia, se hace necesario un procedimiento capaz de calcular la vida a fatiga del elemento en estas circunstancias, siendo uno de los más simples el método de daño acumulado de Palgrem- Miner. Este método consiste en el cálculo del número de ciclos hasta el fallo para cada ciclo definido y su posterior introducción en la ecuación (7-8). En esta expresión Ni es el número de ciclos hasta el fallo para el ciclo ‘i’ y ni representa el número de repeteciones que se producen de este ciclo. D refleja el daño acumulado a medida que se introduce el número de ciclos hasta el fallo de los diferentes ciclos y, cuando D alcanza el valor unidad, se produce el fallo. Con este método es fácil obtener el daño acumulado de los ciclos producidos dentro de un vuelo completo y calcular el número de vuelos hasta el fallo realizando la inversa del resultado obtenido (1/D). 𝐷=𝑛1 𝑁1+𝑛2 𝑁2+⋯𝑛𝑖 𝑁𝑖+⋯ (7-8) 7.5.1 Definición de los ciclos de fatiga Para aplicar la regla de Palgrem-Miner es necesario definir previamente los ciclos de fatiga que soporta el elemento. Estos ciclos reflejan en el caso a estudio las dos situaciones consecutivas de tensiones que se emplean para aplicar el método de Fatemi-Socie y es necesario definir tanto el plano crítico en que será aplicado como las situaciones mencionadas, para lo que se emplea el procedimiento ‘Rainflow’ y lo que se ha denominado en este análisis como procedimiento de plano crítico global. 7.5.1.1 Procedimiento de plano crítico global Hasta ahora se ha comentado el procedimiento para la obtención del plano crítico en una fase, sin embargo, en un vuelo completo tanto este plano como el punto crítico donde es calculado pueden variar. Para realizar todo Aplicación al Modelo 58 Figura 7-13. Resultado de la aplicación de la ventana de Hanning. Para evitar, en la medida de lo posible, frecuencias parásitas provocadas por el fenómeno de aliasing se eliminan las frecuencias superiores a la mitad de la frecuencia máxima de muestreo. Con el fin de comprobar la idoneidad de los procesos aplicados se realiza la anti-transformada de la señal en el dominio de la frecuencia mediante la anti-transformada de Fourier expresada en la ecuación (7-24). Como se puede apreciar en la Figura 7-14, incluso para un número elevado de frecuencias contempladas no se consigue una señal similar a la original. 𝑥𝑙=∑𝑋𝑘𝑒𝑗2𝜋𝑘𝑙 𝑁 𝑁−1 𝑘=0 (7-24) Figura 7-14. Compración de la señal original y la antitransformada tras la aplicación de la ventana de Hanning y la supresión de las frecuencias superiores a la mitad de la frecuencia máxima de muestreo. 59 Estos resultados permiten afirmar que la simplificación de la señal no es posible y por lo tanto la única opción restante es introducir estos datos directamente en la simulación. Para ello se emplea la herramienta ‘Acceleration’ dentro de Ansys y se reduce el rango de tiempo de los datos escogidos a 5 segundos para disminuir el tiempo de simulación. Una vez calculadas las tensiones y el punto crítico, se aplica el método de conteo de ciclos utilizando la tensión máxima principal en este punto para generar el diagrama de picos y valles y posteriormente se realiza el procedimiento a fatiga comentado anteriormente. Los resultados reflejan que estas vibraciones no afectan a la vida a fatiga de la pieza al obtenerse un número infinito de ciclos hasta el fallo y, teniendo en cuenta que además son datos obtenidos de un rotor de mayor envergadura, se puede concluir que la suposición inicial sobre la influencia nula de las vibraciones es correcta. 7.7 Resultados Comprobados los efectos inocuos de las vibraciones en este análisis, los únicos parámetros que influyen en las variaciones de carga son los propios de cada fase del vuelo. Aunque también se pueden generar ciclos de tensiones producidos por los transitorios entre cada fase, el lapso entre cada una de ellas (aproximadamente 2.5 segundos) y las altas frecuencias a las que se aplica el amortiguamiento no permiten obtener ningún ciclo de carga durante el cambio de fase, es decir, la variación se produce de manera cuasi estática para el sistema. Por lo tanto, para calcular los ciclos a fatiga sólo se tendrán en cuenta las fases de vuelo que, debido a lo comentado anteriormente, se pueden analizar por separado. Además, para comprobar la fiabilidad de los datos a fatiga obtenidos, se comprarán con el resultado de aplicar los dos modelos uniaxiales expuestos en 6.5. En primer lugar, se lleva a cabo la co-simulación de cada fase empleando los parámetros expuestos en el apartado 7.4.1 para calcular todas las tensiones necesarias para el análisis. Las tensiones resultantes se exponen en el Anexo A. Una vez obtenidas las tensiones y los puntos críticos de cada fase, se realiza el procedimiento de obtención de la vida a fatiga considerando un ciclo de la forma reposo-fase de vuelo-reposo. De esta manera se determina que la fase que provoca más daño es la fase de despegue, al ser la que devuelve un número de ciclos hasta el fallo menor. En la Tabla 7-11 se comparan los resultados obtenidos para cada fase mediante el procedimiento de Fatemi-Socie y los métodos uniaxiales. Como se puede apreciar, los métodos uniaxiales arrojan en todos los casos un resultado superior, lo que indica que este procedimiento predice una vida hasta el fallo que, mediante el empleo de métodos más sencillos, sería considerada segura, y por lo tanto su utilidad queda contrastada. Tabla 7-11. Resultados del análisis de vida a fatiga para cada fase individual. Fase Nº ciclos método Fatemi- Socie Nº ciclos método multiaxial simplificado Nº ciclos método uniaxial puro Despegue 1.8370 104 2.8272 105 1.4476 105 Transición 1º parte 8.1711 104 2.9693 106 1.3217 106 Transición 2º parte 5.1709 105 3.7801 107 1.7331 107 Crucero 5.0428 105 3.6555 107 1.6763 107 Crucero lento 6.9011 105 5.5166 107 2.5389 107 Crucero ante aterrizaje INF INF INF Parece claro que se debe emplear el punto crítico en el despegue (Figura 7-15) y su ángulo crítico, sin embargo, para todas las fases excepto ésta, el punto crítico se encuentra en la misma posición (Figura 7-16) por lo que, al ser el número de puntos reducido, se calculará la vida a fatiga para ambos. En los dos puntos críticos la fase más severa es el despegue, y en consecuencia se emplea en cada uno de ellos el plano crítico obtenido en esta fase. La secuencia de cargas, en este caso la secuencia de fases de vuelo, es relevante a la hora de calcular la vida a fatiga. Para este caso se ha supuesto la secuencia descrita en Tabla 7-12 que comienza y finaliza en una fase de Aplicación al Modelo 60 reposo. Tabla 7-12. Secuencia de fases empleada para el cálculo de vida a fatiga. Fase Descripción 0 Reposo 1 Despegue 2 Transición 1ª Parte 3 Transición 2ª Parte 4 Crucero 5 Crucero lento 6 Crucero ante aterrizaje 7 Reposo Figura 7-15. Punto crítico de la fase de despegue. Figura 7-16. Punto crítico del resto de fases. 7.7.1 Punto crítico de la fase de despegue (PCD) Tanto los ángulos de Euler como el ángulo del plano crítico en corrdenadas principales se obtienen de la ya 61 comentada fase de despegue, cuyos valores se expresan en la Tabla 7-13. Tabla 7-13. Ángulos de Euler y ángulo del plano crítico para el punto PCD. α β γ Θ -86,879 ° 0,26953 ° 2,5352 ° 41.5 ° Calculadas las tensiones normales y tangenciales en este punto para todas las fases, se crea el diagrama empleado para el método ‘Rainflow’ a partir de las tensiones principales máximas. Como se comentó en el apartado 7.5.1.2, se podrían utilizar otros parámetros para la generación de este diagrama, sin embargo, todos los procedimientos empleados proporcionan el mismo resultado. En la Figura 7-17 se muestra la gráfica original de variación de cargas (línea discontinua roja) y el resultado de la obtención de picos y valles (línea continua azul) utilizando la tensión máxima principal en el punto crítico y la inversa del número de ciclos hasta el fallo (daño) de la Tabla 7-11 como variables para la generación del diagrama. Figura 7-17. Diagramas Rainflow empleados en el punto crítico de la fase de despegue. Se puede observar que sólo existen en este diagrama 3 picos/valles originados por las fases de vuelo debido a las características de la secuencia empleada en el análisis. Por otra parte, el primer ciclo creado a partir del modelo ‘Rainflow’ es el que comprende las fases 3 y 4, sin embargo, la diferencia de tensiones entre estas situaciones es tan pequeña que el número de ciclos hasta el fallo obtenido por Fatemi-Socie es infinito, por lo que para el cálculo del daño acumulado sólo se tiene en cuenta un ciclo formado por las fases de reposodespegue-reposo. Esto implica que el daño de todo el vuelo analizado es el que se produce con un ciclo de este tipo, arrojando una vida a fatiga de 18800 vuelos. 7.7.2 Punto crítico del resto de fases (PCR) Para este punto es necesario calcular de nuevo el plano crítico obtenido en la fase de despegue, obteniéndose los ángulos recogidos en la Tabla 7-14. Tabla 7-14. Ángulos de Euler y ángulo del plano crítico para el punto PCR. α β γ Θ -94,419 ° 14,395 ° -176,44 ° 41.5 ° Aplicación al Modelo 62 El diagrama de ‘Rainflow’ empleado en este apartado varía ligeramente del obtenido con la tensión principal máxima del apartado anterior, que se calcula para un punto distinto, sin embargo, ambos diagramas son similares y en esta situación también se considera como único ciclo dañino el provocado por el despegue. Teniendo en cuenta que en este punto el daño provocado por esta fase es menor, no es extraño que la vida obtenida a fatiga obtenida en este caso sea mayor, con un número de ciclos de vuelo hasta el fallo de 36500. 63 8 CONCLUSIONES os resultados obtenidos de la realización del análisis completo sobre la pieza estudiada arrojan un número de ciclos hasta el fallo que se encuentra dentro del rango esperarado para un sistema de este tipo en las circunstancias estudiadas, sin embargo, el objeto de este proyecto es la generación de un procedimiento que sirva como base para la determinación de la vida a fatiga bajo unas condiciones similares, por lo que estos resultados no se deben contemplar como una solución final, sino como la primera iteración de la investigación en este ámbito del proyecto global. Gran parte de las variables contempladas en el estudio dependen de multitud de factores que se ha optado por simplificar con el fin de abarcar un análisis global, que comprende desde el modelado 3D del sistema hasta el cálculo de la vida a fatiga, pasando por la simulación de las cargas aerodinámicas y la obtención de tensiones mediante elementos finitos. Variaciones en la secuencia de cargas o en la severidad de las fases que la componen pueden afectar de manera signitficativa a los resultados, así como la modificación de los parámetros del material empleado. En este caso se ha empleado una secuencia muy simple para facilitar la interpretación de los resultados, pero el estudio de las fases de vuelo es muy importante para la obtención del número de ciclos hasta el fallo. Como ejemplo, si se incluyera una fase de vuelo rápido o de emergencia tras la fase de crucero lento, el número de ciclos hasta el fallo disminuiría hasta los 9422 ciclos, suponiendo en esta fase una velocidad de giro y de avance de 12600 RPM y 20 m/s respectivamente, y un ángulo de paso de 9.38º (manteniendo el resto de parámetros de la Tabla 7-5). Figura 8-1.Diagrama de Rainflow incluyendo una fase adicional tras la fase de crucero lento. Por otra parte, se ha realizado un estudio de vibraciones que, aunque ha sido útil para exponer el tratamiento de la señal necesario, no ha arrojado ninguna influencia significativa en el daño producido a la pieza. Esta afirmación es correcta para el análisis llevado a cabo con los datos disponibles, pero debería ser estudiada con cautela en posteriores investigaciones, ya que las vibraciones originan generalmente una disminución apreciable de la vida a fatiga. Como se ha comentado en el apartado 7.3, la influencia de los fenómenos expuestos no es apreciable para los parámetros de vuelo inicialmente estimados, sin embargo, no se puede llegar a esta conclusión sin antes comprobar sus efectos para los parámetros finales. Se comprueba que para la fase de despegue (más dañina) la deformación máxima es de 0.35 %, por lo que, al ser inferior al 2%, no existe el riesgo de fenómenos relacionados con la no linealidad geométrica. Asimismo, incluyendo las fuerzas de inercia debidas a la rotación de rotor, la variación de tensiones es mínima, siendo el número de ciclos hasta el fallo con esta L Conclusiones 64 situación de 18214 vuelos en lugar de los 18800 obtenidos con anterioridad. Los resultados y procedimientos expuestos en este documento son consecuencia de análisis que no siempre han contribuido de manera directa al avance del estudio. Algunos ejemplos se han mostrado en el documento, como la consideración de una malla hexáedrica o el análisis de vibraciones comentado anteriormente. Sin embargo, otros procedimientos no han sido mencionados durante el documento con el fin de simplificar la compresión de este. Son dignos de mención aspectos como la definición del punto y el plano crítico, que han requerido de de varias iteraciones hasta alcanzar el resultado expuesto, o el análisis con elementos finitos realizado mediante el programa Ansys. Éste último se ha perfeccionado durante la realización del proyecto y ha necesitado de varios procesos de ensayo y error, así como la consulta de diversas fuentes tanto para el desarrollo de un modelo adecuado como para la inclusión de algunos factores analizados como el amortiguamiento, el estudio de vibraciones o los fenómenos considerados en el párrafo anterior. En este aspecto también cabe resaltar las modificaciones realizadas durante el proyecto sobre los métodos y zonas de aplicación de esfuerzos. Otro de los procedimientos cuyo análisis se ha extendido durante el estudio es la cosimulación, que ha requerido de varias iteraciones en el modelo de Adams y diversas modificaciones en el modelo de fuerzas al tiempo que la comprensión de los datos y la precisión de los parámetros iba aumentando como consecuencia de la investigación en este ámbito, que se puede observar comparando los datos de vuelo empleados en los primeros apartados y los parámetros finales expuestos en el capítulo 7. De manera similar, el procedimiento de cálculo de la vida a fatiga ha evolucionado desde un modelo sencillo (Findley) y unos ciclos de vuelo simples, hasta la aplicación de un modelo más complejo (F-S) y la introducción del procedimiento de ‘Rainflow’ en el análisis de los ciclos de carga. Por útlimo cabe resaltar que no todos los procedimientos explicados durante el análisis son aplicables de manera general. Al igual que se ha realizado un estudio de las posibles influencias de diversos factores, puede que en otros casos el análisis deba expandirse y contemplar otros o, por el contrario, descartar algunos de los comentados. Varias consideraciones a este respecto se presentan en el siguiente capítulo. 65 9 FUTURAS LÍNEAS DE INVESTIGACIÓN as suposiciones realizadas a lo largo del proyecto y las limitaciones encontradas en algunos de los aspectos que se han tratado impiden que el estudio llevado a cabo se pueda considerar como un análisis completamente fiable para la determinación de la vida a fatiga. Algunos de estos condicionantes que requieren de un estudio a mayor profundidad, así como ciertas posibilidades de expansión del procedimiento tratado en este trabajo se exponen a continuación: • Sistema y modelado. Para simplificar y agilizar la realización del estudio se ha optado por estudiar únicamente una de las piezas que componen el sistema completo. Aunque se trata de un elemento muy importante y susceptible al fallo, sería necesario realizar un análisis completo del resto de piezas para determinar aquella cuya vida a fatiga sea limitante, es decir, que el número de ciclos hasta el fallo en esta pieza sea el menor de todas las estudiadas. Además, se ha procurado simular la presencia del tornillo, pero sería necesario establecer si las suposiciones realizadas son correctas. • Modelo de elementos finitos. En el apartado 5.1.2 se comenta la limitación numérica debida a la licencia académica de Ansys empleada. La influencia de esta limitación en los resultados debe ser analizada a través de un modelo que permita indicar de manera precisa la convergencia de la malla. Por otra parte, la aplicación de las condiciones de contorno se ha realizado de manera aislada sin simular los contactos y sujeciones entre las distintas partes y por lo tanto sin establecer la influencia de éstos. • Modelo de determinación del fallo. El modelo de Fatemi-Socie empleado es ampliamente utilizando en la industria, pero su implementación requiere de un estudio de la aeronave real que refleje que los resultados obtenidos concuerdan con la vida a fatiga del modelo. La determinación del punto y plano crítico se ha realizado desde una perspectiva analítica que necesita ser confirmada bien por modelos más avanzados o por análisis experimentales. Por otra parte, durante el trabajo se han obviado aspectos tales como el acabado superficial, defectos de fabricación o la temperatura de servicio cuya influencia debería reflejarse y analizarse en futuras investigaciones. También se ha comentado durante el desarrollo del trabajo la existencia de modelos modificados, como por ejemplo el modelo de Kim & Park [5], cuya aplicación podría considerarse más adecuada en función de las características del análisis. • Modelo de fuerzas aerodinámicas. El modelo de fuerzas proporcionado ha permitido simular las condiciones básicas de un vuelo para poder obtener así los esfuerzos generados, sin embargo, se trata de un modelo genérico que puede no concordar de manera exacta con los esfuerzos que se generarían en el caso real. Por lo tanto, se requiere del estudio del modelo para, si fuera necesario, modificar ciertos aspectos de éste. También hay que tener en cuenta que el modelo de pala utilizado durante el análisis puede no ser el más adecuado para este tipo de condiciones de servicio debido, por ejemplo, a los fenómenos de autorrotación que se producen, lo que implica que necesariamente ha de ser modificado el modelo en este aspecto y a su vez el programa de generación de esfuerzos. • Fases de vuelo. Con el fin de simplificar la interpretación de los resultados se ha optado por la creación de un patrón de vuelo reducido y simple. En el apartado anterior se ha comentado la gran influencia que tiene la modificación de las fases de vuelo y por ende se debe hacer hincapié en el estudio del comportamiento de la aeronave a lo largo de un ciclo de vuelo real para reflejar de manera más adecuada los parámetros de vuelo y los esfuerzos que se generan en cada una de ellas. • Vibraciones. En este estudio no ha sido posible realizar un estudio in situ sobre las vibraciones ocurridas durante un vuelo genérico de la aeronave a estudio. Aunque se han empleado los datos de las pruebas de un rótor de mayor envergadura para una primera estimación, es necesario analizar las vibraciones reales que se producen en el sistema para determinar su influencia en el fallo de los distintos elementos. L Futuras Líneas de Investigación 66 • Fuerzas centrífugas y no-linealidades geométricas. Aunque se ha comentado que estos aspectos no influyen de manera significativa en el estudio llevado a cabo, no se debe extrapolar esta conclusión para todos los casos y la influencia de los fenómenos mencionados debería ser analizada en cada caso particular. • Palas empleadas. En el apartado 7.4.1.2 se introduce el fenómeno de autorrotación que se produce a causa del empleo de palas planas por lo que se debería analizar la geometría alar más adecuada que optimice el empuje y minimice los esfuerzos bajo estas condiciones. • Amortiguamiento y modos de vibración. Para la determinación del rango de frecuencias en el que se aplica el amortiguamiento estructural se han empleado los dos primeros modos de vibración del elemento obtenidos a partir del programa Ansys. Sin embargo, estos modos se han calculado para elemento aislado, sin tener en cuenta la posible influencia del resto de elementos como las palas, el rotor o incluso la propia aeronave. Aunque finalmente este parámetro no ha influido de manera sustancial en el análisis, es conveniente estudiar con una visión más global la forma más adecuada para incorporarlo ya que puede suponer una gran diferencia a la hora de analizar los posibles transitorios que se produzcan en el sistema. Una forma de estudiar el valor del amortiguamiento, así como las frecuencias para las que se produce es el denominado como método ‘Moving-block’ [17], utilizado ampliamente en la industria aeronáutica a este efecto. Otro aspecto a tener en cuenta en este ámbito es el efecto denominado como ‘whirl flutter’. Este fenómeno ocurre en rotores que tienen la capacidad de modificar el eje de rotación de las palas y conduce a una perturbación inicial desde el rótor hacia las alas que se retroalimenta y puede conducir a modos de vibración inestables [17]. • Optimización topológica. Además de determinar el fallo de una pieza o sistema, una de las posibles aplicaciones del procedimiento descrito puede ser la optimización topológica. El objetivo del empleo de esta técnica es la eliminación durante el diseño de las zonas que soportan menores tensiones, con el consecuente ahorro de peso y material durante la fabricación. Se trata de una ciencia en auge que incluso algunos programas ya incluyen por defecto pero que no suele tener en cuenta el fallo a fatiga y, por lo tanto, el procedimiento descrito en este trabajo puede ser de gran utilidad a la hora de determinar la viabilidad de un diseño optimizado. • Fabricación. Este ámbito del desarrollo del rotor estudiado no se ha analizado durante el documento y puede influir de manera considerable en ciertos aspectos tratados ya que las características propias de la fabricación (calidad, método de fabricación…) pueden afectar significativamente a la hora de estudiar la viabilidad del diseño, como se ha expuesto en este mismo capítulo. Por lo tanto, se requeriría de un estudio en profundidad a este respecto, incluso la fabricación de un prototipo, que proporcione los datos necesarios para su análisis. • Modelo tripala. El diseño del rotor utilizado está basado en el modelo bipala expuesto en el capítulo 2 pero, aunque se podría implementar en la aeronave manteniendo esta geometría, podría ser interesante analizar los beneficios e inconvenientes que supondría el empleo de un modelo similar pero tripala. 67 A. ANEXO A n este anexo se indican las tensiones principales, normales y tangenciales obtenidas como resultado de la simulación de cada fase de vuelo mediante el modelo de elementos finitos para cada punto critico correspondiente, calculado de la forma expuesta en 7.1. Estas tensiones son empleadas para el cálculo de la vida a fatiga mediante los métodos de Fatemi-Socie y uniaxiales cuyos resultados se encuentran en la Tabla 7-11. Tabla A-1. Tensiones del punto en crítico en cada fase estudiada de forma aislada. σ1 σ2 σ3 σx σy σz τxz τxy τyz Despegue 2.13 108 1.33 107 6.30 106 1.39 107 2.12 108 6.70 106 -4.67 105 -1.09 107 9.14 106 Transición 1ª parte 1.82 108 1.3 107 7.9 106 1.38 107 1.81 108 8.91 106 1.89 106 7.99 106 1.02 107 Transición 2º parte 1.49 108 1.16 107 6.82 106 1.14 107 1.49 108 7.61 106 1.40 106 3.63 106 8.22 106 Crucero 1.50 108 1.16 107 6.83 106 1.14 107 1.49 108 7.63 106 1.41 106 3.71 106 8.25 106 Crucero lento 1.45 108 1.13 107 6.67 106 1.11 107 1.44 108 7.44 106 1.33 106 3.04 106 7.95 106 Crucero ante aterrizaje 8.48 107 6.49 106 3.13 106 6.36 106 8.45 107 3.60 106 8.63 105 1.86 106 4.96 106 E Tensión [Pa] Fase 75 C. ANEXO C ste anexo contiene los programas de Matlab empleados en el modelo de fuerzas aerodinámicas. Además del programa utilizado en la co-simulación, también se expone el modelo utilizado para el cálculo de las variables de vuelo del apartado 7.4.1. Como ya se ha comentado, estos programas se tratan como una caja negra y, por lo tanto, sólo incluyen breves anotaciones para facilitar su uso práctico y las expresiones utilizadas para el cálculo de las variables. Código del modelo de fuerzas de la co-simulación La primera parte de este código, que se encuentra dentro del bloque de fuerzas del programa creado en Simulink (Figura 3-14), define las variables de entrada (velocidad de giro, ángulo de paso y tiempo de simulación) y de salida (fuerzas y pares en cada división de pala). function [T1,T2,T3,T4,T5,T6,T7,T8,T9,T10,T11,T12,T13,T14,T15,T16,T17,Q1,Q2,Q3,Q4,Q5,Q6 ,Q7,Q8,Q9,Q10,Q11,Q12,Q13,Q14,Q15,Q16,Q17,vc] = modelofuerzas(Omega_degs,theta_deg,time) c = 0.03; Ct=0; delta0 = 0.0085; delta1 = 0; delta2 = 0.263; Cl0 = 0; Cl1 = 5.73; A continuación, se introduce el valor de diversas características de la aeronave: masa (m), número de palas (b), mitad de la envergadura total del rótor (R) y densidad del fluido colindante (ρ, ‘rho’ en el código). m = 2.600; b = 2; R = 0.13585; rho=1.225; g = 9.81; % Gravedad A = pi*R^2; % Área x = linspace(0.1,1-1e-6,100); M_tip_max = 0.314; Para definir las condiciones de vuelo, se introducen los parámetros a partir a las entradas de este bloque para el ángulo de paso (ϑ, ‘theta’ en el código) y la velocidad de giro (ω, ‘Omega’ en el código ), y se introduce manualmente la velocidad de avance (VC). % User input conditions theta = theta_deg*pi/180; rad_s2RPM = 60/(2*pi); Omega=Omega_degs*pi/180; Omega_RPM=Omega*rad_s2RPM; E Anexo C 76 vc=5; En esta parte del código se realizan los cálculos necesarios para la determinación de una primera solución de CT (coeficiente de tracción) mediante una teoría simplificada (TCMEP) que servirá como valor inicial para la teoría empleada (TCMEP_P). Las ecuaciones (C-1), (C-2) y (C-3) y la Tabla C-1 reflejan las operaciones y variables utilizadas para conseguir este primer valor. 𝜆= −(𝜎(𝑥)𝐶𝑙𝛼(𝑥) 16 −𝜆𝑐 2)+√(𝜎(𝑥)𝐶𝑙𝛼(𝑥) 16 −𝜆𝑐 2)2+𝜎(𝑥)𝐶𝑙𝛼(𝑥) 8𝜃(𝑥)𝑥 (C-1) 𝑑𝐶𝑇=𝜎𝐶𝑙𝛼 2(𝜃−𝜆 𝑥)𝑥2𝑑𝑥 (C-2) 𝐶𝑇=∫𝑑𝐶𝑇 1 0 (C-3) Tabla C-1. Parámetros empleados en TCMEP. Variable Definición 𝜆 Relación de entrada de corriente 𝜎 Solidez 𝜆𝑐=𝑉𝑐(𝜔𝑅) ⁄ Relación de entrada de corriente debido a velocidad axial 𝐶𝑙 Coeficiente de sustentación 𝜃 Ángulo de paso 𝑥 Distancia adimensionalizada desde el eje de giro lambda_c = vc/(Omega*R); %TEORÍA TCMEP lambda = @(x) -(solidez*Cl1/16 - lambda_c/2)+ sqrt((solidez*Cl1/16 - lambda_c/2)^2 + solidez*Cl1/8*theta.*x); dCT_dx = @(x) 0.5*solidez*Cl1*(theta.*x.^2-lambda(x).*x); Ct = quadgk(@(x) dCT_dx(x),1e-6,1-1e-6) %TEORÍA TCMEP_P lambda = @(x) lambda_Prandtl(x,solidez,Cl1,theta,b,Ct,lambda_c); La teoría ‘TCMEP_P’ incluye las pérdidas de punta de pala y ráiz propuestas en el modelo de Prandtl. Según este módelo, debe corregirse dCT por un factor F, cuya expresión se muestra en la ecuación (C-4). 𝐹(𝑥)=(2 𝜋)cos−1(𝑒−𝑓(𝑥)) (C-4) 77 Obtenida λ mediante la ecuación (C-6), localizada al final del código, se procede a calcular el resto de variables. Las ecuaciones (C-7), (C-8), (C-9), (C-10), (C-11) y (C-12) reflejan las expresiones empleadas en esta sección para el cálculo de los coeficientes de tracción y par (CQ) de cada división de pala. 𝐶𝑇=∫ 1 2𝜎𝐶𝑙(𝜃−𝜆 𝑥)𝑥2 1 0𝑑𝑥 (C-7) 𝑐𝑜𝑛𝑣𝑇=𝜌𝐴(𝜔𝑅)2 (C-8) 𝑑𝐶𝑃𝑖= 𝜆𝐶𝑇 (C-9) 𝑑𝐶𝑃0= 𝜎 2(𝛿0+𝛿2(𝜃−𝜆 𝑥)2)𝑥3 (C-10) 𝑑𝐶𝑃=𝑑𝐶𝑄=𝑑𝐶𝑃𝑖+𝑑𝐶𝑃0 (C-11) 𝑐𝑜𝑛𝑣𝑃 𝜔=𝑐𝑜𝑛𝑣𝑄= 𝐴𝜔2𝑅3 (C-12) C_T_TEP = quadgk(@(x) 1/2*solidez*Cl1*(theta.*x.^2-lambda(x).*x),1e-6,1- 1e-6) alpha = @(x) theta - lambda(x)./x; dCT_dx = @(x) 0.5*solidez*Cl1*(theta.*x.^2-lambda(x).*x); conv_T = rho*A*(Omega^2)*(R^2); 𝑓(𝑥)=𝑏 2(1−𝑥 𝜆(𝑥)) (C-5) 𝜆2+(𝜎𝐶𝑙𝛼 8𝐹 −𝜆𝑐)𝜆−𝜎𝐶𝑙𝛼 8𝐹 𝜗𝑥=0 (C-6) Anexo C 78 dT_dx = @(x) conv_T*0.5*solidez*Cl1*(theta.*x.^2-lambda(x).*x); dCPi_dx = @(x) lambda(x).*dCT_dx(x); dCPo_dx = @(x) 0.5*solidez*(delta0 + delta2*(theta - lambda(x)./x).^2).*x.^3; dCP_dx = @(x) dCPi_dx(x)+ dCPo_dx(x); dCQ_dx = @(x) dCPi_dx(x) + dCPo_dx(x); conv_P = rho*A*(Omega^3)*(R^3); dP_dx = @(x) conv_P*(dCPi_dx(x)+ dCPo_dx(x)); conv_Q = rho*A*(Omega^2)*(R^3); dQ_dx = @(x) conv_Q*(dCPi_dx(x)+ dCPo_dx(x)); % Divisiones de pala x1=1/17; x2=2*x1; x3=3*x1; x4=4*x1; x5=5*x1; x6=6*x1; x7=7*x1; x8=8*x1; x9=9*x1; x10=10*x1; x11=11*x1; x12=12*x1; x13=13*x1; x14=14*x1; x15=15*x1; x16=16*x1; Ct1 = quadgk(@(x) dCT_dx(x),1e-6,x1-1e-6); Ct2 = quadgk(@(x) dCT_dx(x),x1,x2-1e-6); Ct3 = quadgk(@(x) dCT_dx(x),x2,x3-1e-6); Ct4 = quadgk(@(x) dCT_dx(x),x3,x4-1e-6); Ct5 = quadgk(@(x) dCT_dx(x),x4,x5-1e-6); Ct6 = quadgk(@(x) dCT_dx(x),x5,x6-1e-6); Ct7 = quadgk(@(x) dCT_dx(x),x6,x7-1e-6); Ct8 = quadgk(@(x) dCT_dx(x),x7,x8-1e-6); Ct9 = quadgk(@(x) dCT_dx(x),x8,x9-1e-6); Ct10 = quadgk(@(x) dCT_dx(x),x9,x10-1e-6); Ct11 = quadgk(@(x) dCT_dx(x),x10,x11-1e-6); Ct12 = quadgk(@(x) dCT_dx(x),x11,x12-1e-6); Ct13 = quadgk(@(x) dCT_dx(x),x12,x13-1e-6); Ct14 = quadgk(@(x) dCT_dx(x),x13,x14-1e-6); Ct15 = quadgk(@(x) dCT_dx(x),x14,x15-1e-6); Ct16 = quadgk(@(x) dCT_dx(x),x15,x16-1e-6); Ct17 = quadgk(@(x) dCT_dx(x),x16,1-1e-6); Cq1 = quadgk(@(x) dCQ_dx(x),1e-6,x1-1e-6); Cq2 = quadgk(@(x) dCQ_dx(x),x1,x2-1e-6); Cq3 = quadgk(@(x) dCQ_dx(x),x2,x3-1e-6); Cq4 = quadgk(@(x) dCQ_dx(x),x3,x4-1e-6); Cq5 = quadgk(@(x) dCQ_dx(x),x4,x5-1e-6); Cq6 = quadgk(@(x) dCQ_dx(x),x5,x6-1e-6); Cq7 = quadgk(@(x) dCQ_dx(x),x6,x7-1e-6); 79 Cq8 = quadgk(@(x) dCQ_dx(x),x7,x8-1e-6); Cq9 = quadgk(@(x) dCQ_dx(x),x8,x9-1e-6); Cq10 = quadgk(@(x) dCQ_dx(x),x9,x10-1e-6); Cq11 = quadgk(@(x) dCQ_dx(x),x10,x11-1e-6); Cq12 = quadgk(@(x) dCQ_dx(x),x11,x12-1e-6); Cq13 = quadgk(@(x) dCQ_dx(x),x12,x13-1e-6); Cq14 = quadgk(@(x) dCQ_dx(x),x13,x14-1e-6); Cq15 = quadgk(@(x) dCQ_dx(x),x14,x15-1e-6); Cq16 = quadgk(@(x) dCQ_dx(x),x15,x16-1e-6); Cq17 = quadgk(@(x) dCQ_dx(x),x16,1-1e-6); Cp = quadgk(@(x) dCQ_dx(x),1e-6,1-1e-6); Calculados los coeficientes adimensionales obtenidos como resultado de la integración en cada división de pala, se procede a la obtención de las fuerzas y pares mediante las ecuaciones (C-13) y (C-14). Nótese que los resultados son divididos por el número de palas (b) ya que los valores resultantes reflejan la fuerza y par total de todo el sistema. 𝑇=𝑐𝑜𝑛𝑣𝑇𝐶𝑇 (C-13) 𝑄=𝑐𝑜𝑛𝑣𝑄𝐶𝑄 (C-14) % Total values P = conv_P*Cp % total power (Watts) T1 = conv_T*Ct1/b; % total thrust (N) T2 = conv_T*Ct2/b; T3 = conv_T*Ct3/b; T4 = conv_T*Ct4/b; T5 = conv_T*Ct5/b; T6 = conv_T*Ct6/b; T7 = conv_T*Ct7/b; T8 = conv_T*Ct8/b; T9 = conv_T*Ct9/b; T10 = conv_T*Ct10/b; T11 = conv_T*Ct11/b; T12 = conv_T*Ct12/b; T13 = conv_T*Ct13/b; T14 = conv_T*Ct14/b; T15 = conv_T*Ct15/b; T16 = conv_T*Ct16/b; T17 = conv_T*Ct17/b; Q1 = conv_Q*Cq1*1000/b; % total torque (N-mm) Q2 = conv_Q*Cq2*1000/b; Q3 = conv_Q*Cq3*1000/b; Q4 = conv_Q*Cq4*1000/b; Q5 = conv_Q*Cq5*1000/b; Q6 = conv_Q*Cq6*1000/b; Q7 = conv_Q*Cq7*1000/b; Q8 = conv_Q*Cq8*1000/b; Q9 = conv_Q*Cq9*1000/b; Anexo C 80 Q10 = conv_Q*Cq10*1000/b; Q11 = conv_Q*Cq11*1000/b; Q12 = conv_Q*Cq12*1000/b; Q13 = conv_Q*Cq13*1000/b; Q14 = conv_Q*Cq14*1000/b; Q15 = conv_Q*Cq15*1000/b; Q16 = conv_Q*Cq16*1000/b; Q17 = conv_Q*Cq17*1000/b; function lambda = lambda_Prandtl(x,solidez,Cl1,theta,b,C_T,lambda_c) for i = 1:length(x) f = @(lambda) b/2*((1-x(i))/lambda); F = @(lambda) 2/pi*acos(exp(-f(lambda))); lambda(i) = fzero(@(lambda) -(solidez*Cl1/16/F(lambda)-lambda_c/2) + sqrt((solidez*Cl1/16/F(lambda) - lambda_c/2)^2+... solidez*Cl1/8/F(lambda)*theta*x(i))-lambda,[0,10*sqrt(C_T/2)]); end end end Código del modelo de fuerzas para el cálculo de parámetros Este código expone el método empleado para el cálculo de los parámetros durante la definición de las fases de vuelo. A diferencia del método anterior, en este caso no se calcula el empuje, sino que se impone el necesario para contrarrestar el peso y la resistencia aerodinámica, obteniéndose como resultado el ángulo de paso requerido para la velocidad angular utilizada. Las teorías de cálculo son idénticas a las expuestas en el apartado anterior al igual que la nomenclatura empleada, pero incluyendo las variables necesarias en este caso, mostradas en el apartado 7.4.1. c = 0.02753; c = 0.03; delta0 = 0.0085; delta1 = 0; delta2 = 0.263; Cl0 = 0; Cl1 = 5.73; m = 2.600; b = 2; R = 0.13585; g = 9.81; A = pi*R^2; ISA = ISAmodel; theta_1 = 0; rho = ISA.density(0); x = linspace(0.1,1-1e-6,50); solidez = b*c/(pi*R); delta_t = 0.90; A continuación, se indican las variables de vuelo y se calcula la resistencia aerodinámica. Nótese que en este código aparecen las variables necesarias tanto para el vuelo vertical como para el vuelo horizontal. 81 D_in =(R*2)*100/2.54; RPM_max_APC = 150000/D_in; n_max_APC = (RPM_max_APC/60); Omega_max_APC = 2*pi*(RPM_max_APC/60); V_tip_max_APC = Omega_max_APC*R; M_tip_max_APC = V_tip_max_APC/ISA.a0 M_tip_max = M_tip_max_APC; V_tip = ISA.a0*M_tip_max; Omega = delta_t*V_tip/R; rad_s2RPM = 60/(2*pi); Omega_RPM = Omega*rad_s2RPM; %V_c = 5; % Vuelo vertical V_c = 20; % Vuelo horizontal S_ref = 0.45/4; %Superficie alar CD = 0.5; %Coeficiente de resistencia aerodinámica CD0 = 0.02; % Coeficiente de resistencia parasitaria CL = 0.5; % Coeficiente de sustentación %CL = 2; % Coeficiente de sustentación en tansición K = 0.04; % Coeficiente de resistencia inducida CD_cruise = CD0 + K*CL^2; q = 0.5*rho*(V_c^2); D = q*CD*S_ref; D_cruise = q*CD_cruise*S_ref; % Vuelo vertical %T = m*g + D; % Vuelo horizontal T = D_cruise; lambda_c = V_c/(Omega*R); theta = @(theta_0,x) theta_0 + theta_1*x; Es en las siguientes dos líneas de código es donde se encuentra la principal diferencia con el caso anterior, ya que se define el coeficiente de empuje necesario en lugar de ser calculado por el programa. CT_ideal_per_engine = T/(2*rho*A*(Omega^2)*R^2) C_T = CT_ideal_per_engine; for i = 1:1 if i == 1 metodo = 'TCMEP'; else metodo = 'TCMEP_P'; end switch metodo case 'TCMEP' lambda = @(theta_0,x) -(solidez*Cl1/16 - lambda_c/2)+ sqrt((solidez*Cl1/16 - lambda_c/2)^2 + solidez*Cl1/8*theta(theta_0,x).*x); case 'TCMEP_P' lambda = @(theta_0,x) lambda_Prandtl(theta_0,x,solidez,Cl1,theta,b,C_T,lambda_c); otherwise error('El método seleccionado no existe') end C_T_TEP = @(theta_0) quadl(@(x) 1/2*solidez*Cl1*(theta(theta_0,x).*x.^2- lambda(theta_0,x).*x),1e-6,1-1e-6); Anexo C 82 F = @(theta_0) C_T-C_T_TEP(theta_0); guess = 6*C_T/solidez/Cl1+3/2*sqrt(C_T/2)-3/4*theta_1; theta_0 = fzero(F,0.7); alpha = @(theta_0,x) theta(theta_0,x)-lambda(theta_0,x)./x; dCT_dx = @(theta_0,x) 0.5*solidez*Cl1*(theta(theta_0,x).*x.^2- lambda(theta_0,x).*x); conv_T = rho*A*(Omega^2)*(R^2); dT_dx = @(theta_0,x) conv_T*0.5*solidez*Cl1*(theta(theta_0,x).*x.^2- lambda(theta_0,x).*x); dCPi_dx = @(theta_0,x) lambda(theta_0,x).*dCT_dx(theta_0,x); dCPo_dx = @(theta_0,x) 0.5*solidez*(delta0+delta2*(theta(theta_0,x)- lambda(theta_0,x)./x).^2).*x.^3; dCP_dx = @(theta_0,x) dCPi_dx(theta_0,x)+ dCPo_dx(theta_0,x); dCQ_dx = @(theta_0,x) dCPi_dx(theta_0,x)+ dCPo_dx(theta_0,x); conv_P = rho*A*(Omega^3)*(R^3); dP_dx = @(theta_0,x) conv_P*(dCPi_dx(theta_0,x)+ dCPo_dx(theta_0,x)); conv_Q = rho*A*(Omega^2)*(R^3); dQ_dx = @(theta_0,x) conv_Q*(dCPi_dx(theta_0,x)+ dCPo_dx(theta_0,x)); Ct = quadl(@(x) dCT_dx(theta_0,x),1e-6,1-1e-6) Cp = quadl(@(x) dCP_dx(theta_0,x),1e-6,1-1e-6) Cq = quadl(@(x) dCQ_dx(theta_0,x),1e-6,1-1e-6) P_1 = rho*A*(Omega^3)*(R^3)*Cp P_2 = quadl(@(x) dP_dx(theta_0,x),1e-6,1-1e-6) T_1 = conv_T*Ct Q_1 = conv_Q*Cq lambdaio = sqrt(C_T/2); solidezId = b*c/(pi*R); theta_pId = 4*C_T/(solidezId*Cl1)+lambdaio; dCT_dxId = @(x) 0.5*solidezId*Cl1*(theta_pId-lambdaio)*x; theta_xId = @(x) theta_pId./x; dCPi_dxId = @(x) lambdaio*dCT_dxId(x); dCPo_dxId = @(x) 0.5*solidezId*(delta0*x.^3+delta2*(theta_pId- lambdaio)^2*x); C_pId = quadl (@(x) dCPi_dxId(x)+dCPo_dxId(x),0,1) alphaId = @(x) theta_xId(x) -lambdaio./x; alphaOp = sqrt(delta0/delta2); solidez_pOp = 4*C_T/Cl1/alphaOp; solidezOp = @(x) solidez_pOp./x; theta_xOp = @(x) alphaOp+lambdaio./x; dCT_dxOp = @(x) 0.5*solidez_pOp*Cl1*alphaOp*x; dCPi_dxOp = @(x)lambdaio*dCT_dxOp(x); dCPo_dxOp = @(x) 0.5*solidez_pOp*(delta0+delta2*alphaOp^2)*x.^2; C_pOp = quadl (@(x) dCPi_dxOp(x)+dCPo_dxOp(x),0,1) figure (1) hold on 83 plot(x,lambda(theta_0,x)) xlabel('x') ylabel('\lambda') grid on if i == 2 hold on plot(x,lambdaio*ones(size(x)),x,lambdaio*ones(size(x))) legend('TCMEP','TCMEP_P','Ideal','Óptimo') end figure(2) hold on plot(x,theta(theta_0,x)*180/pi) xlabel('x') ylabel('\theta') grid on if i == 2 hold on plot(x,theta_xId(x)*180/pi,x,theta_xOp(x)*180/pi) legend('TCMEP','TCMEP_P','Ideal','Óptimo') end figure(3) hold on plot(x,dCT_dx(theta_0,x)) xlabel('x') ylabel('dC_T/dx') grid on if i == 2 hold on plot(x,dCT_dxId(x),x,dCT_dxOp(x)) legend('TCMEP','TCMEP_P','Ideal','Óptimo') end figure(4) hold on plot(x,dCPi_dx(theta_0,x)) xlabel('x') ylabel('dC_P_i/dx') grid on if i == 2 hold on plot(x,dCPi_dxId(x),x,dCPi_dxOp(x)) legend('TCMEP','TCMEP_P','Ideal','Óptimo') end figure (5) hold on plot(x,dCPo_dx(theta_0,x)) xlabel('x') ylabel('dC_P_o/dx') grid on if i == 2 hold on plot(x,dCPo_dxId(x),x,dCPo_dxOp(x)) legend('TCMEP','TCMEP_P','Ideal','Óptimo') end figure(6) Referencias 90 Tunnerl Models, Pennsylvania State University, 2015. [18] C. S. Gómez, Diseño y Optimización Topológica de Rockers para un modelo de Formula Student, Sevilla, 2017. [19] C. S. Fernández, Co-simulación para el sistema de extensión/retracción hidráulico de un tren de aterrizaje, Sevilla, 2013. [20] «eFatigue,» Altair, [En línea]. Available: https://www.efatigue.com/. [21] ANSYS, «ANSYS e-learning,» [En línea]. Available: https://caeai.com/resources/tips-tricks-hex-brick- meshing-ansys-e-learning.