Desarrollo de herramienta de diagnóstico para el turbohélice PT6-A

Proyecto Fin de Carrera Ingeniería de Telecomunicación Formato de Publicación de la Escuela Técnica Superior de Ingeniería Autor: F. Javier Payán Somet Tutor: Juan José Murillo Fuentes Dep. Teoría de la Señal y Comunicaciones Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2013 Proyecto Fin de Carrera Ingeniería Aeroespacial Desarrollo de Herramienta de Diagnóstico para el Turbohélice PT6-A Autor: Inés Zamora Fernández Tutor: Francisco J.Jiménez-Espadafor Aguilar Dpto. Ingeniería Energética Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2022 Proyecto Fin de Carrera Ingeniería Aeroespacial Desarrollo de Herramienta de Diagnóstico para el Turbohélice PT6-A Autor: Inés Zamora Fernández Tutor: Francisco J.Jiménez-Espadafor Aguilar Catedrático de Universidad Dpto. Ingeniería Energética Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2022 Proyecto Fin de Carrera: Desarrollo de Herramienta de Diagnóstico para el Turbohélice PT6-A Autor: Inés Zamora Fernández Tutor: Francisco J.Jiménez-Espadafor Aguilar El tribunal nombrado para juzgar el trabajo arriba indicado, compuesto por los siguientes profesores: Presidente: Vocal/es: Secretario: acuerdan otorgarle la calificación de: El Secretario del Tribunal Fecha: Agradecimientos Querría comenzar dando las gracias a mi tutor del proyecto, Francisco J.Jiménez Espadafor Aguilar. Su paciencia y guía han sido de gran valor y ayuda para mí en esta trayectoria que ha sido el final de la carrera y el trabajo fin de grado. Su implicación y disponibilidad han sido siempre de extrema ayuda en situaciones complicadas y problemas que parecían no tener salida. Ha sabido enseñarme infinidad de conceptos nuevos y su tutoría me ha facilitado ver con otra perspectiva algunos obstáculos en el camino. Seguidamente, me gustaría agradecer el apoyo de mi compañero de trabajo José Carlos. Con él he compartido muchas horas, ideas y reflexiones. Entre los dos, hemos podido superar las dificultades que se nos han presentado, en forma de códigos numéricos, y complementar las ideas de uno con las del otro aunque estuviéramos en lugares distintos. A mi familia y amigos, querría agradecerles su paciencia y apoyo. Sin sus ánimos durante los altibajos del curso, probablemente no habría sido capaz de llegar a donde estoy hoy. En especial doy gracias a mis padres que, desde la distancia, han sabido siempre valorar mi esfuerzo y celebrar mis éxitos. A Almudena, mi amiga y compañera en mi experiencia en el extranjero. Con ella he compartido todo y nada, lágrimas y risas, preocupaciones y evoluciones. A Alejandro, que siempre me recuerda el valor de mi trabajo. Con él he podido evolucionar, reír y compartir este curso tan especial. A las personas que me han hecho tan feliz en Italia y que me han recordado que, aun cuando uno pone lo académico en primer plano, los pequeños ratos en familia son medicina. Inés Zamora Fernández Sevilla, 2022 I Resumen Es te proyecto fin de grado surge del interés del diseño, modelado y posterior diagnóstico de un motor aeronáutico para su posible posterior utilización en la práctica, como ayuda al proceso de revisión y verificación tan imprescindible en los componentes que conforman un avión de transporte civil. El objetivo principal del proyecto es, a parte del modelaje y el estudio del motor, el desarrollo de una herramienta digital que permita diagnosticar un fallo interno y su magnitud correspondiente, desarrollando así la capacidad de poder detectar la degradación de una parte del motor sin ni siquiera comenzar a inspeccionarlo de forma física. En primera instancia se llevará a cabo un estudio preliminar de los sensores del motor, sus características principales y el estudio del avión al cual propulsa. También se explicará, en detalle, el modelo utilizado para la determinación de las magnitudes termodinámicas partícipes en el cálculo. El proyecto se irá desarrollando de manera ordenada y progresiva ya que constituye, en su totalidad, un proceso que se ha ido construyendo en la misma manera en la que se presenta en esta memoria. Es entonces, objeto de estudio la resolución del ciclo en punto de diseño, los mapas de las turbomáquinas y el posterior ciclo fuera de diseño, con el objetivo de modelar el funcionamiento del motor bajo consideración digitalmente. En la última etapa, en donde reside el principal objetivo del proyecto, se hará un estudio a cerca de cómo se introduce la degradación en el motor, la utilización de la herramienta de diagnóstico y, a continuación, un estudio estadístico explorativo sobre su precisión y bondad. Debido a la magnitud de estudio que está bajo el alcance de este proyecto, durante el mismo se han hecho suposiciones y simplificaciones, las cuales son siempre justificadas y explicadas. Además, el estudio estadístico está acotado únicamente a la degradación de un cierto componente pero, aun así, se reflexionará acerca de las implicaciones que esto tiene y las posibles mejoras que podrían aplicarse. III XNotación IFkFactor de interferencia del componente k cCuerda del perfil tEspesor del perfil τiTrabajo específico de ”i” hEntalpía CpCalor específico a presión constante ¯cCuerda media Swet Superficie bañada RekNúmero de Reynolds del componente k AÁrea TnTracción necesaria LSustentación Ar Argón NNitrógeno H2OAgua pura CO2Dióxido de carbono OOxígeno VVolumen KKelvin PrRelación de presiones T00 Temperatura de remanso en vuelo T01 Temperatura de remanso a la entrada del difusor T02 Temperatura de remanso a la entrada del compresor T21 Temperatura de remanso a la salida del primer escalonamiento de axial T22 Temperatura de remanso a la salida del segundo escalonamien- to de axial T23 Temperatura de remanso a la salida del tercero escalonamiento de axial T24 Temperatura de remanso a la salida del cuarto escalonamiento de axial T25 Temperatura de remanso a la salida del escalonamiento centrífugo T03,Tcc Temperatura de remanso a la salida de la cámara de combustión T04 Temperatura de remanso a la salida del primer escalonamiento de turbina de alta T045 Temperatura de remanso a la salida del primer escalonamiento de turbina de baja/potencia T05 Temperatura de remanso a la salida del segundo escalonamien- to de turbina de baja/potencia T06 Temperatura de remanso a la salida del difusor de salida T0is Temperatura de remanso a la salida del componente ”i” especificado, en condiciones isentrópicas. TiTemperatura a la salida del componente ”i” especificado NVelocidad de giro del eje Nc Velocidad corregida ˙mGasto másico del motor ˙mfGasto másico del carburante mcGasto corregido ηRendimiento de un componente ΦEntropía a presión constante sEntropía πcRelación de compresión P00 Presión de remanso en vuelo P0iPresión de remanso a la salida del componente ”i” Notación XI PiPresión a la salida del componente ”i” πcc Relación de presiones de remanso en la cámara de combustión Wpar Potencia parásita EeEmpuje específico TSFC Thrust Specific Fuel Consumption cEConsumo específico fDosado εPorcentaje de gasto de aire destinado a sangrado Xr,iFracción molar del componente ”i” de los reactivos X p,i Fracción molar del componente ”i” de los productos después de la combustión Tre f Temperatura de referencia Pre f Presión de referencia Nre f Revoluciones de giro de referencia γre f Coeficiente adiabático de referencia γiCoeficiente adiabático del componente ”i” β Parámetro de identificación de posición en mapas de turbomáquinas JParámetro de avance de la hélice CpCoeficiente de potencia de la hélice ai,j Coeficiente dependiente de la temperatura para el cálculo de las magnitudes termodinámicas DhDiámetro de la hélice DDiámetro ΨFactor de carga WeiTrabajo específico de compresor ”i” P.DEn punto de diseño O.DFuera del punto de diseño (Off Design) rmedio Radio medio de un álabe RxError introducido en la degradación del ciclo PDF Probability Density Function (Función de Densidad ee Probabilidad) ηrel Error relativo 1 Introducción Un proyecto de diagnóstico y análisis tecnológico Con este proyecto se intenta conseguir una herramienta efectiva de diagnóstico para un motor de la gama PT6-A o similar. Esta herramienta permitirá analizar y localizar, según unos datos de entrada al motor, qué zonas del mismo están fallando debido a diferentes grados de degradación que se impondrán previamente. El esquema principal del proyecto se puede claramente visualizar en la Figura (1.1), en donde se deja como objeto de mejora futura la estrategia de diagnóstico GPA y de redes neuronales, utilizando el método de optimización. El proyecto consiste de dos grandes bloques, necesarios para el desarrollo de los códigos en Matlab que permitan resolver ciclos y la posterior implementación de estos para analizar casos de fallo del motor. Como primer paso se necesitarán conocer los sensores en servicio del motor: dónde se encuentran, qué miden, cómo funcionan; para poder tener una idea de en qué puntos del ciclo del motor tendremos datos y en cuáles no. Debido a la escasez de datos experimentales de motores de este tipo de libre acceso al público, esta información se utilizará de forma simplificativa y con datos generales. A continuación se construirá en Matlab un modelo del motor elegido aproximado del ciclo de diseño (Design Point) que, con unas condiciones de entrada que se especificarán, proporcione todas las salidas del motor necesarias para su estudio. Complementario a este, se realizará otro código que resuelva el ciclo de este mismo motor pero fuera de diseño (Off Design) en donde, con unas entradas requeridas que serán distintas, veremos cómo reacciona el motor a estar bajo unas condiciones fuera de su punto de operación. A continuación, comenzará el segundo bloque del trabajo en donde, a partir del modelo del motor de referencia que reproduce las prestaciones de nuestro motor en cualquier condición operativa, se creará un modelo perturbado. Las curvas características del motor se parametrizarán para poder distorsionarlas y generar fallos en el motor modelado, que representen cómo se degrada un motor después de muchas horas en servicio. Con los datos de operación que nos devuelva el motor, a través de una estrategia de diagnóstico, podremos evaluar e identificar la degradación que sufre cada uno de los elementos. 1 2Capítulo 1. Introducción Figura 1.1 Diagrama conceptual del trabajo. La revisión y diagnóstico de un motor aeronáutico hoy en día se compone de distintas fases, después de su fabricación. El motor suele ser insertado en una célula de prueba cerrada en donde se le hacen pruebas de banco y se registran los datos que devuelve. Es también necesario someter a los rotores a pruebas de equilibrio y vibración, además de controles de aceite, combustible y fugas de gas. Se le añade a esto pruebas en situaciones desfavorables como la ingesta de granizo, polvo, nieve y desechos, además de pruebas en donde el motor ya es instalado en el ala y es testado en el hangar exterior. Por último, una vez en servicio, el motor debe ser revisado cada cierto tiempo, 5 años usualmente, para su correcta manutención y puesta a punto. Aun siendo este tipo de pruebas completamente necesarias y correctas para hacer del motor un elemento seguro, no es tan común aplicar técnicas digitales para probar la eficiencia del motor. Es cierto que, las pruebas experimentales tienen el beneficio de contar con datos reales de la máquina, pero careciendo del componente inteligente que tiene una herramienta de diagnóstico digital. Es por ello que, a estas pruebas experimentales, sería de gran ayuda agregarle una herramienta de tipo computacional como la que introducimos en este proyecto fin de grado. No solo sería útil para la primera inspección del motor sino, para las frecuentes revisiones a las que se someterá este a lo largo de su vida operativa, así haciendo la labor de inspección y búsqueda de fallos en el interior mucho menos laboriosa y más focalizada. Figura 1.2 Prueba en banco de motor aeronáutico. 2 Modelo en el punto de diseño 2.1 Caracterización del motor y sensores en servicio 2.1.1 Familia de motores PT6A La familia de motores elegida, de la que forma parte nuestro motor modelo, es la de Pratt & Whitney PT6A. Los motores que la componen son turbohélices aeronáuticos, que son conocidos particularmente por su elevada fiabilidad y tiempo de duración medio. Su variante principal es la A, en donde podemos encontrar motores de aviones en una amplia gama de potencias y tamaños. Adicionalmente, pueden encontrarse variantes que gozan de menos auge como la B o la C, los cuales son turboejes para helicópteros. Contando con más de 60 modelos, cada motor de la gama PT6A se utiliza para propósitos distintos: transporte de personas comercial y privado, entrenamiento militar, transporte de mercancías en condiciones climatológias adversas y misiones especiales. Esta en actual servicio para ambos la Fuerza Aérea y el Ejército de los Estados Unidos, además de estar en uso en más de 170 países. Debido a la gran variedad de motores que ofrecen, las buenas prestaciones y fiabilidad, se ha decidido escoger el motor bajo estudio dentro de esta familia. 2.1.2 Motor bajo estudio El motor que se ha elegido para el análisis y estudio posterior es el Pratt & Whitney PT6A-65B. Este es un tipo de motor que comenzó su aplicación en el abril del año 1982 y que recibió la certificación de la EASA en septiembre de 1984. Debido a sus prestaciones, nivel tecnológico y alta potencia, lo hemos considerado adecuado para nuestro estudio. Este motor es clasificado por la EASA como un modelo III. En términos generales, es un turbo-prop de turbina libre que cuenta con 4 escalonamientos de compresor axial y uno centrífugo, una única cámara de combustión anular, un solo escalonamiento de turbina del generador de gas y dos de la de potencia. A continuación se explicitarán algunas de sus características y prestaciones. •Longitud total:1490 mm •Diametro total: 464 mm •Peso seco específico (specific dry weight): 225 kg •Potencia contínua máxima: 875 KW •Sistema de control del combustible: hidromecánico •Potencia en despegue (5 minutos): 875 KW 3 4Capítulo 2. Modelo en el punto de diseño •Sangrado máximo 1. Externo: -5.25 2. Durante el arranque: -0.68 kg/min •Temperaturas ITT (Inter Turbine Temperature) 1. Temperatura máxima contínua: 810 ºC 2. Temperatura en despegue (5 minutos): 820 ºC 3. Starting temperature: 1000ºC 4. Temperatura maxima contínua de operación del aceite: 110 ºC 5. Temperatura minima de operación del aceite (suelo): 110 ºC •Velocidad del rotor del generador de gas: 39000 rpm •Revoluciones de la turbina de potencia: 1700 rpm (100%). Correspondencia con 29894 rpm. •Revoluciones de la turbina de potencia en transient: 1870 (110%) •Datos importantes al nivel del mar 1. Potencia de eje equivalente: 1249 hp 2. Potencia de eje: 1173 hp 3. Empuje del chorro: 189 lb 4. Potencia final de turbina (output): 1700 rpm 5. Revoluciones del generador de gas: 39000 rpm •Datos importantes en despegue 1. Potencia de eje equivalente: 1249 hp 2. Potencia de eje: 1173 hp 3. Empuje del chorro: 189 lb 4. Potencia final de turbina (output): 1700 rpm 5. Revoluciones del generador de gas: 39000 rpm •Reversa maxima 1. Eje: 900 hp 2. Output: 1650 rpm •Aeronaves que utilizan este motor 1. Beechcraft 1900 2. Beechcraft 1900C •Caja de reducción 1. Ratio= 0.0568:1 2.1.3 Sensores en servicio Con el objetivo de hacer un estudio paramétrico del ciclo del motor tanto en el punto de diseño como fuera, necesitamos conocer la posición de los sensores en servicio del motor. La medida de las temperaturas, presiones, velocidades y demás magnitudes en el interior de la máquina es crucial para hacer un estudio por etapas lo más cercano a la realidad posible. Sin embargo, no es tarea sencilla encontrar toda esta información de manera precisa, por tanto, se estimará, de manera lógica y coherente, la posición de aquellos sensores de los que no pueda encontrarse información suficiente. Para comenzar insertaremos un esquema del interior del motor completo para poder visualizarlo por partes (Figura 2.1). Principalmente, se pueden apreciar en el esquema dos zonas: la azul y la roja. 2.1 Caracterización del motor y sensores en servicio 5 La zona azul es la sección de compresión en donde primero encontramos la caja de reducción accesoria la cual es dirigida o rotada por el eje y movida por el compresor y que tiene como función dar potencia a sistemas auxiliares o "accesorios" como a las bombas de aceite y combustible, el FCU y el generador de esfuerzo de torsión (torque). A su izquierda está el compresor que, en nuestro caso, cuenta con 4 escalonamientos axiales y uno centrífugo final para proporcionar mayor relación de compresión total. En color blanco, entre el compresor y la turbina de baja (o del compresor) se encuentra la cámara de combustión y todo el sistema de ignición. A continuación se encuentra la turbina de baja, la cual solo tiene un escalonamiento. Esta está en color azul ya que alimenta a la zona de compresión del generador de gas. La zona roja, de derecha a izquierda, cuenta con dos escalonamientos de turbina de potencia los cuales proporcionan potencia a la hélice de nuestro turbo-prop. Es fácil visualizar en el esquema el eje (shaft) por dónde se trasmite dicha potencia y las dos etapas de caja de reducción. Estas reducen la velocidad rpm de la turbina a una más conveniente para la hélice en cuestión. También puede observarse que la tobera de salida es casi perpendicular al eje del motor (en color blanco) debido a que todo el empuje se generará por la tracción de la hélice y no por el empuje a chorro de la tobera de salida, teniendo, como comentaremos más adelante el parámetro αya definido: α=1. Figura 2.1 Esquema interno del motor PT6A-65B. En la Figura (2.2) también se pueden apreciar bien cada una de las partes, que están señalizadas y nombradas. 6Capítulo 2. Modelo en el punto de diseño Figura 2.2 Esquema interno del motor PT6A-65B con nombres. De manera que, una vez visto el esquema principal, presentamos una imagen Figura (2.3) en donde se indica la ubicación de los sensores principales del motor. Podemos observar que en la zona de compresión, en la primera sección, se encuentra un sensor de la velocidad del rotor de compresión (compressor speed) y, con él, un tacómetro generador (”tacho generator”). El tacómetro generador es un pequeño generador eléctrico que contiene dos polos dentro del campo magnético generado por un rotor magnético que gira a velocidad proporcional a la velocidad del compresor (o turbina de potencia). El generador produce una onda sinusoidal de un ciclo por revolución y da un voltaje de aproximadamente 21 voltios cuando gira a 4200 RPM, que es un estándar en la industria. Contando los ciclos junto con la proporcionalidad de los engranajes sabremos las RPM del compresor y turbina de potencia. En la misma sección inicial encontramos un indicador de la temperatura y presión del aceite utilizado para la lubricación del sistema. El sensor ITT (”Inter Turbine Temperature”) es un indicador que se encuentra justo entre la sección azul y roja de la Figura (2.1) y, por tanto, entre la turbina del generador de gas y la turbina de potencia. Este dispositivo no mide la máxima temperatura del motor (que correspondería con la que tendríamos a la salida de la cámara de combustión y, por tanto, entrada a la turbina del generador) si no a una temperatura interna dentro de la zona de turbina. Está compuesto de 10 pares termoeléctricos de cromel/alumel unidos en paralelo a través del "bus bar" y con salida al exterior a través de un cable blindado ("wiring harness") con conexión para la instrumentación en la cabina del piloto. Cuando la temperatura aumenta, se genera un voltaje mayor en cada uno de los pares termoeléctricos. Debido a que la distribución de calor no es uniforme, los pares ven temperaturas diferentes y generan voltajes distintos. Para obtener una medida media, estos se conectan en paralelo y por tanto la medida generada es la media de temperatura en 10 posiciones distintas. Por esto mismo la temperatura que se lee en el ITT puede no ser justo la de la sección 5, y por ello se corrige con unos ciertos tests, finalmente teniendo la temperatura medida corregida en un indicador en cabina, cómo mencionado antes. Otro indicador es el de esfuerzo de torsión (torque) el cual es un sistema hidromecánico para medir la potencia generada por el motor; está situado dentro de la carcasa trasera de la caja de reducción. El mecanismo consiste de un cilindro, un pistón, sellos, válvulas y muelle. El funcionamiento de este indicador se basa en la medición de la torsión aplicada a la hélice. Esta torsión produce un pequeño movimiento rotacional en la primera caja de reducción anular de la hélice. Esta empuja la válvula de control y el pistón y, por tanto, permite que aumente la presión de aceite y que el pistón vuelva a presionar sobre la caja de reducción. Este movimiento rotacional se para cuando hay un balance entre la fuerza de la caja y la presión de aceite. 2.1 Caracterización del motor y sensores en servicio 7 Por último, nos encontramos el indicador de la velocidad de la hélice (Propeller speed), que nos da la medida en revoluciones por minuto. Figura 2.3 Posición de sensores en el esquema de un motor genérico de la serie PT6-A según el manual oficial de entrenamiento de Pratt and Whitney.. A parte de los indicadores mostrados en la imagen, hay algunos otros sistemas indicadores y adicionales que podemos mencionar y que acompañaremos con algunas imágenes. A diferencia de los anteriores estos no serán objeto de estudio del trabajo ya que no proporcionan información útil para el ciclo del motor o el control de las variables principales. Por ello, solo se mencionarán brevemente: • La regulación y filtración del aceite lubricante se hace a través del termómetro de bypass y la válvula de chequeo del aceite (Figura 2.4). • El sistema del combustible el FCU (Fuel Control Unit) controla la temperatura del combustible (Figura 2.5). • El chip detector proporciona indicación de partículas metálicas contaminantes en el sistema de aceite. Se encuentra en la carcasa o carenado exterior situado cercano a las cajas de reducción. •En la Figura (2.6) se encuentra el sistema de ignición automático. • Se utiliza una válvula de sangrado para evitar problemas en el compresor y liberar la presión 2.5 al exterior como se observa en Figura (2.7) en el número 2. • El sistema "pre-swirl" se utiliza para cambiar el ángulo de ataque del aire que entra al compresor. Mejora el margen de "surge" a velocidades del rotor bajas (Figura 2.8 ) 14 Capítulo 2. Modelo en el punto de diseño Figura 2.13 Esquema de ala con estrechamiento y sus partes. ¯ctotal =¯c1A1+¯c2A2 A1+A2 (2.4) Figura 2.14 Cuerda media aerodinámica total para un estrechamiento doble. ¯c=2 3·(c3 r−c3 t) (c2 r−c2 t)(2.5) Figura 2.15 Cálculo de cuerda media aerodinámica siendo cr y ct valores de cuerda en la raiz y la punta del tramo considerado. Ala de doble estrechamiento Primer estrechamiento •Longitud=2.4813 m •b=2.316 m •¯c1=2.316 m •A1=5.746 m2 Segundo estrechamiento •Longitud=4.9626 m •b2=1.0752 m •¯c2=1.7713 m •A2=8.4146 m2 Ala completa •¯cw=1.9923 m •t cw=0.7444 2.316 =0.3214m Timón vertical •Longitud=2.399 m 2.2 Modelo de avión en donde se integra el motor 15 •cr=3.06 m •ct=2.315 m •¯cv=2.357 m •t cv=0.331 3.06 =0.10816m Timón horizontal •Longitud=2.4813 m •cr=1.489 m •ct=0.827 m •¯ch=1.1895 m •t ch=0.3 1.489 =0.20147m Fuselaje •Longitud =lB=17.7m •NB=lN=15.715=longitud de referencia de fuselaje •NBtop =2.150m=longitud cabeza del fuselaje •Dmax=1.654 m •l d=9.501 Motor •Nengine =4.466m=longitud referencia del motor •Dmax=0.496 m •l d=9.004 2.2.2 Polar parabólica del avión: Cálculo de CD0 El cálculo de la polar parabólica de una aeronave concreta es complejo y, para poder hacerlo de manera precisa, necesitaríamos datos reales de la propia aeronave. Sin embargo, haré una aproximación de la misma, debido a la escasez de datos disponibles. La expresión de la polar parabólica en la Figura (2.16), pone en evidencia la necesidad de calcular el término cD0, el cual se calculará con el método del área parásita equivalente. CD=CD0+kC2 L(2.6) Figura 2.16 Ecuación de la polar parabólica. El método del área parásita equivalente calcula el coeficiente de resistencia a sustentación nula con la expresión (2.17). CD0=1.1 Swing ∑fk(2.7) Figura 2.17 Coeficiente CD0 donde en la expresión se multiplica por 1.1 para contar con contribuciones misceláneas. 16 Capítulo 2. Modelo en el punto de diseño fk=Cf k ·CFk·IFk·FFk·Swetk(2.8) Figura 2.18 Factores debidos a contribuciones varias a la resistencia aerodinámica. Cada uno de estos términos a calcular se definen a continuación: Cfk=0.455 log10 ·Rek 2 .58 (2.9) Figura 2.19 Coeficiente de fricción para flujo turbulento en placa plana. Rek=ρVlk µ(2.10) Figura 2.20 Número de Reynolds para el flujo para cada componente del avión k, siendo lk la longitud de referencia de cada componente k del avión. Longitud de referencia de cada componente •B (Body) ⇒lB •W (Wing) ⇒¯cw •H (Horizontal Tail) ⇒¯ch •V (Vertical Tail) ⇒¯cv •Nwor NB(Wing or Body Nacelles) ⇒lN • T (Tip tanks) ⇒lT (En nuestro caso, el avión no tiene "tip tanks" o tanques de carburante en la punta del ala) Figura 2.21 Factores de interferencia de cada componente. 2.2 Modelo de avión en donde se integra el motor 17 Figura 2.22 Calculo de los factores de forma FFk . Estos factores incluyen los efectos de grosor donde el término t c es el "thickness ratio" del perfil (grosor máximo dividido entre la cuerda) y l d es el "fineness ratio" del perfil (longitud dividida entre diámetro máximo). Figura 2.23 Esquema del avión simplificado vista en alzado. Figura 2.24 Esquema del avión simplificado vista en planta. Generatriz del cono= 2.315m y radio del cono= 0.827m. Finalmente, el área Swetk es el área externa de la superficie del avión que es tocada por al aire. Puede ser calculada utilizando figuras geométricas simples: las superficies sustentadoras se pueden considerar placas 18 Capítulo 2. Modelo en el punto de diseño planas y los cuerpos son conos o cilindros. Una aproximación de la simplificación de estas áreas se puede visualizar en la Figura (2.23, 2.24). De esta forma, las áreas consideradas serán consideradas como escrito en las ecuaciones (2.11, 2.12, 2.13, 2.14, 2.15, 2.16, 2.17) Swettotal =Swetv+SwetH+SwetT+Swetnacellesnody +Swetnacelleswing +Swetw(2.11) Swetv=2.399 ·3.06 =7.3409m2(2.12) SwetH=2.4813 ·1.489 ·2=7.3893m2(2.13) SwetT=0(2.14) Swetnacellesbody =2·π·0.827 ·2.315 =12.029m2(2.15) Swetnacelleswing =2·π·0.248 ·4.466 +π·0.248 ·0.661 =7.474m2(2.16) Swetw=28.8m2(2.17) 2.2.3 Cálculo del coeficiente de Oswald (e) y el coeficiente de resistencia inducida (k) Para completar la expresión de la polar parabólica, necesitamos calcular el coeficiente k, llamado coeficiente de resistencia inducida. En las expresiones (2.18, 2.19, 2.20) se calculan tanto este k, como el coeficiente e (coeficiente de Oswald), necesario en el cálculo del precedente. Debemos tomar bajo consideración que en la expresión (2.19), "b" se refiere a la envergadura completa del avión y Salar a la superficie alar completa en planta, considerando el área que estaría ocupando el fuselaje también. e=1.78(1−0.045(AR)0.68)−0.64 =0.79987 (2.18) AR =b2 Salar =8.3858 (2.19) K=1 πARe =0.04745 (2.20) 2.2.4 Modelo de cálculo de empuje demandado por aeronave Una vez calculados los coeficientes requeridos para la polar parabólica, podemos plantear el modelo que utilizaremos en Matlab para calcular el empuje que nos demanda nuestro avión para cada cota y altura. Asignado un vehículo (W,S, CD0 , k), una cota de vuelo (h o ρ ) y una velocidad de vuelo, podemos calcular la tracción que necesita el avión para un vuelo en configuración de crucero equilibrado y estacionario. Nos basamos en un esquema del avión como en el de la Figura (2.25). Las ecuaciones que gobiernan la dinámica de esta condición de vuelo son las especificadas en las ecuaciones que se muestran abajo. De esta forma, calculando previamente CD0 con el método que se describió anteriormente y calculando CL con los datos que 2.3 Modelo ’Gas Path’con fluidos reales: programación del ciclo en el punto de diseño 19 tenemos asignados, se puede calcular la tracción necesaria ( Tn ) lo que, en nuestro caso, es lo mismo que el empuje demandado. Esta función es añadida en los cógidos de Matlab. Figura 2.25 Fuerzas aerodinámicas en vuelo en crucero. Tn=D=1 2ρV2S(CD0+kC2 L) = Å1 2ρSCD0ãV2+Ç2kW2 ρSå1 V2=AV2+B V2=E(2.21) L=W=1 2ρV2SCL(2.22) CL=2(W/S) ρV2(2.23) 2.3 Modelo ’Gas Path’con fluidos reales: programación del ciclo en el punto de diseño Para poder obtener un mapa parametrizado de prestaciones del turbohélice en el punto de diseño, debemos primero modelar numéricamente en Matlab el ciclo de un motor para unos datos de entrada. Esto lo llevaremos a cabo con un modelo ’gas path’con fluidos reales. 2.3.1 Bases del modelo El modelo que se va a utilizar está basado en unas ciertas suposiciones y simplificaciones que se indican a continuación: • Se suponen los gases ideales y, por tanto, válida la ecuación de los gases ideales (pV=nRT) con propiedades dependientes de la temperatura. R será una constante ya que el fluido no alcanzará condiciones que requieran la consideración de la presión reducida. • Se considerará el peso molecular real de los fluidos de la reacción de combustión, diferente para reactivos y productos. • La composición del aire seco (reactivos) será (% en peso): N2 (75,47), O2 (23,19), Ar(1,29), CO2 (0,05) •El combustible utilizado en nuestro programa es C10H17 •El poder calorífico es 43MJ/kg referido a 298 K •Los gases de combustión productos son CO2,H2O,O2,Ar,N2 •Consideramos flujo estacionario y unidimensional •Configuración de motor en eje único 20 Capítulo 2. Modelo en el punto de diseño 2.3.2 Ecuaciones y cálculos para el modelo A continuación se escribirán una serie de ecuaciones que definen las variables de trabajo necesarias para montar el ciclo: h(T) = ZT Tre f Cp(T)dT (2.24) Figura 2.26 Entalpía. Φ(T) = ZT Tre f Cp(T) TdT (2.25) Figura 2.27 Φ, Entropía a presión constante. Pr=exp(Φ−Φ0 R)(2.26) Figura 2.28 Relación de compresión. TdS =dh −vd p (2.27) Figura 2.29 Relación de Gibbs. Φ2−Φ1=Rln P2 P1 ... > P2 P1s=cte =Pr2 Pr1 (2.28) Figura 2.30 Diferencia de Φsolo si el proceso es isoentrópico. Prc=exp(ηpc Φ−Φ0 R)(2.29) Figura 2.31 Relación de compresión. PrT=exp(Φ−Φ0 RηpT )(2.30) Figura 2.32 Relación de expansión. Podemos observar que todas estas variables bajo estudio son dependientes de la temperatura. Para modelar 2.3 Modelo ’Gas Path’con fluidos reales: programación del ciclo en el punto de diseño 21 dicho comportamiento se van a definir una serie de coeficientes que dependen de la temperatura y que modelan las propiedades termodinámicas que tenemos que calcular. Estos coeficientes se presentan en la Figura (2.33) para condiciones estándar de 298,15 K y 100 KPa. Figura 2.33 Coeficientes para propiedades termodinámicas en condiciones estándar . De esta manera, las integrales (2.24, 2.25, 2.26), pueden escribirse con las ecuaciones (2.31, 2.32, 2.33) ˜cp,i ˜ R=ai1+ai2T+ai3T2+ai4T3+ai5T4(2.31) ˜ hi ˜ RT =ai1+ai2 2T+ai3 3T2+ai4 4T3+ai5 5T4+ai6 T(2.32) ˜si ˜ R=ai1lnT+ai2T+ai3 2T2+ai4 3T3+ai5 4T4+ai7(2.33) El programa que modele nuestro ciclo del motor se basará en la definición de estas propiedades del fluido de trabajo, haciendo uso de los coeficientes correspondientes según la temperatura en cada estación, separando los casos T≤1000 yT>1000. Se programarán, por tanto, dos códigos en Matlab para la definición de estas variables. El primero coeftermo.m, que tendrá en su interior para cada elemento, los coeficientes de la tabla (2.33), y el segundo Cp_h_aire.m, en donde llamando al fichero anterior, para una temperatura y fracción molar de las especies, nos proporcione a la salida la h (KJ/Kg) y el Cp (tanto en KJ/KmolK como en KJ/KgK). Con esta base clara, procedemos a programar, en orden, todos los ficheros de cada estación del ciclo: difusor, compresor, cámara de combustión, turbina del compresor, turbina de potencia (de baja o de la hélice) y tobera. Todos los códigos que se mencionan en la totalidad del trabajo se encuentran en los códigos de Matlab que se proporcionarán con el trabajo. 2.3.3 Estaciones del ciclo y su resolución Resolución de forma secuencial Una primera resolución que hicimos para este ciclo en el punto de operación fue secuencial. Esta se basa en la típica resolución utilizada en lecciones teóricas y que, en teoría se podría utilizar para resolver a mano. 22 Capítulo 2. Modelo en el punto de diseño Esto es posible en el ciclo On Design y no en el Off Design, como se verá posteriormente. En esta resolución se pueden ver más claramente las incógnitas y ecuaciones del problema y cómo se resuelve paso por paso. Comenzamos con el DIFUSOR con fichero difusorciclo.m que constituye la estación (0-1). De entrada tenemos las condiciones de vuelo (P0,M0,T0) y la relación de pérdida de presión de remanso ( πdmax ). Se definen los elementos a utilizar y la composición del aire seco, así como las fracciones molares de cada especie. Se calculan los valores de γ0 , v0 y a0 con las propiedades del fluido a la temperatura exterior y se define la entalpía de remanso como en (2.34). ZT00 T0 Cp(T)dT =U2 2⇒h00 =h(T0)+ v2 0 2(2.34) Para definir las condiciones 00, aplicamos la definición de entalpía del modelo donde h00 =h(T00) , definiendo así un comando fsolve que nos daría T00 . Para resolver P00 utilizamos la ecuación del proceso isentrópico (2.35) ZT00 T0 Cp(T) TdT =ZP00 P0 RlogÅP00 P0ã(2.35) Por último, para las condiciones 01, escribiríamos las ecuaciones (2.36, 2.37). T00 =T01 (2.36) P01 =P00 ·πdmax (2.37) El COMPRESOR , que constituye la estación (1-2), se resuelve con el código compresor.m. Tenemos de entrada (T01,P01) del fichero anterior y también (πc,ηc pol). De nuevo, para el aire seco exterior, se calcula h01 y, a continuación se plantea la ecuación del proceso de compresión que se muestra en (2.38), la cual se resuelve con un fsolve para T02 . Esto se hace presuponiendo que T02 ≤1000 , en caso contrario, el programa plantea el mismo cálculo pero para los coeficientes superiores a 1000 K. Esto se hace con un comando if y , a veces, ifelse, que se utilizará de manera recurrente en los sucesivos códigos. ZT02 T01 Cp(T) Tdt =ZP00 P0 RlogÅP01 P02 ã(2.38) Con T02 , calculamos h02 y el trabajo específico del compresor(2.39). Por último P02 se calcula a través de la relación de compresión (2.39,2.40). Así estarían ya calculadas las condiciones 02. τc=h02 −h01 (2.39) P02 =P01πc(2.40) La CÁMARA DE COMBUSTIÓN constituye la estación (2-3), a la cual le corresponden dos ficheros en el programa: ccθ.mycalculo_f.m A la entrada del primer fichero, que sería el principal tendremos ( T0 , T02 , P02 , θ , ηcc , πcc ) siendo: θ el grado de calentamiento (2.31), ηcc el rendimiento de la cámara de combustión y πcc la pérdida de presión de remanso en la cámara de combustión. 2.3 Modelo ’Gas Path’con fluidos reales: programación del ciclo en el punto de diseño 23 Immediatamente ya podemos resolver las condiciones 03 (T03,P03) con las ecuaciones (2.41, 2.42) T03 =T0θ(2.41) P03 =P02πcc (2.42) A parte de estas condiciones 03, necesitamos como salida de la función, las fracciones molares de los productos de la combustión y el dosado (f). Primero, definimos la fórmula molecular del combustible a utilizar que, en nuestro caso, será: C10H17 . También, se definen los coeficientes (n=10 y m=17) como los índices del carbono y el hidrógeno del carburante. La fórmula de la reacción de combustión es necesaria para resolver las fracciones moleculares de las especies que han salido de la cámara de combustión como producto. La reacción de combustión y los flujos másicos de reactivos y productos según la estequiometría se pueden ver en las expresiones siguientes. CnHm+[n+m 4](02+3.717N2+0.04456Ar +0.001568CO2)λrel = [n+λrel[n+m 4]0.001568]CO2+ m 2H2O+λrel3.717[n+m 4]N2+λrel[n+m 4]0.04456Ar +[λrel −1][n+m 4]O2 (2.43) ˙mN2=PMN2 PMf uel ˙mf uelλrel3.717[n+m 4] ˙mAr =PMAr PMf uel ˙mf uelλrel0.04456[n+m 4] ˙mCO2=PMCO2 PMf uel ˙mf uelλrel0.001568[n+m 4] ˙mO2=PMO2 PMf uel ˙mf uelλrel[n+m 4] ˙mN2=PMN2 PMf uel ˙mf uelλrel3.717[n+m 4] ˙mAr =PMAr PMf uel ˙mf uelλrel0.04456[n+m 4] ˙mCO2=PMCO2 PMf uel ˙mf uel[n+λrel[n+m 4]0.001568] ˙mO2=PMO2 PMf uel ˙mf uel[λrel −1][n+m 4] ˙mH2O=PMH2O PMf uel ˙mf uel m 2 Siendo λ=Aire Fuel yλrel =(Aire/Fuel) (Aire/Fuel)esteq Con estas expresiones definidas, sabiendo los pesos moleculares de cada especie y adimensionalizando cada flujo másico con el flujo de aire total ( ˙ma ) se pueden calcular las fracciones molares de los productos. Las expresiones para todo esto se hallan en el fichero ccθ.m . Se puede deducir también que el dosado 30 Capítulo 2. Modelo en el punto de diseño Difusor de salida(5-6) P6=P04 expÇ−(φ04(T04)−φ6s(T6s)) Rgå(2.79) P05 P6 =expÇφ05(T05)−φ6(T6) Rgå(2.80) 0.95 =P05 P04 ·P04 P03 ·πcc P02 P01 ·πtd ·P00 P6 (2.81) P05 =P06 (2.82) Las 22 incógnitas que resolveremos en el ciclo serán: T00,P00,T01,P01,T02,P02,Wcomp,T03,ε,P03,f,T04,P04,πth,T05,P05,πtl,T05sT6s,P6,P06,T6. De esta forma, el sistema tiene 22 ecuaciones y 22 incógnitas el problema queda cerrado. Los iterantes iniciales que se le proporcionan a la función, como comentamos anteriormente, son valores lógicos de cada una de estas variables en punto de operación de un motor de nuestro calibre. Las funciones de Matlab usadas para esta resolución son: ciclo.m, en donde están recogidas todas las ecuaciones del ciclo y ciclosol.m, en donde se llama a la herramienta fsolve para resolver la función anterior, donde se dan los iterantes iniciales y donde, a posteriori, se calculan datos de salida necesarios para poder evaluar las prestaciones del motor. Al ser un ciclo en punto de diseño, no trabajamos con las dimensiones del motor, por lo que todas las salidas serán variables específicas. Las salidas de interés calculadas por la función son, como mencionados anteriormente: Ee, TSFC, CE. Sus expresiones se ilustran en las Figuras (2.83, 2.84, 2.85) Ee=−v0+ηprophηboxηmec ·(1+f)·1−ε 1000·(1−Wparasita ·τtba ja·1 v0 (2.83) TSFC =Åf·(1−ε 1000 ·1 Eeã(2.84) cE=Çf (1+f)·(1−Wparasita)·τhelice å(2.85) De esta forma, el ciclo en punto de diseño nos da como salida estas variables, las presiones y temperaturas en todos los puntos del ciclo,el dosado y las relaciones de expansión de las turbina de baja y alta. 2.3.4 Resultados en diseño y elección del punto de diseño Una vez resuelto el ciclo "On Design", desarrollaremos una serie de gráficas para la correcta visualización de los resultados. Los mapas que se han obtenido no son mapas de operación, los cuales se presentarán más adelante. Los resultados de este ciclo, son mapas en punto de diseño que, escogidas unas variables de entrada ( α,θ,πc ), nos presentan unos datos de salida concretos. Considerando alguna de estas variables constantes y asignando un rango de valores a las restantes, podemos obtener gráficos muy útiles para la interpretación de los resultados. Observaremos que, en el caso de muchos de los gráficos mostrados, α siempre se toma omo α=1 debido a que hemos quitado la tobera y hemos impuesto una relación de presiones constante, que nos hace prescindir de este parámetro de entrada, suponiendo que el chorro de la tobera en realidad no proporciona empuje. 2.3 Modelo ’Gas Path’con fluidos reales: programación del ciclo en el punto de diseño 31 Por otra parte, estos mapas nos servirán para escoger un punto de diseño que optimice las variables calculadas en el ciclo: Ee , TSFC y cE . Este punto elegido debe dar un valor lógico de las mismas para la cota de vuelo escogida en punto de operación y para el avión elegido de referencia. De aquí que los mapas que consideramos no son aún de operación ya que no nos darán información de las prestaciones y comportamiento del motor en cualquier operación, sino que representarán un único punto dentro de estos mapas anteriormente mencionados. La elección del punto de diseño se basa en la optimización de alguna variable del ciclo, o bien en la operación de un avión escogido, que lleve a bordo el motor bajo estudio. Habrá casos en donde, para vuelos cortos o de menor rango, interesará un empuje específico alto y, para aplicaciones de aviación civil y de rango mayor, un TSFC minimizado. Nuestro motor en concreto es comúnmente utilizado en aeronaves regionales, aunque también ha sido utilizado en aeronaves de carga, ejecutivas y de uso militar. Adoptando una actitud generalista, escogemos como aeronave de referencia el Beechcraft 1900C, aeronave de vuelo regional. Por tanto, la elección del punto de diseño se basa en las condiciones de vuelo en crucero de este avión, siendo esta condición operativa la que constituye el mayor tiempo de uso del avión y del motor. De este modo, utilizamos como datos de entrada P0=37.7634KPa , M0=0.4125(v0=127,Û 7m/s,a0= 309.79m/s) , T0=238.83K , correspondientes a los valores de presión, velocidad y temperatura según la Atmósfera Estándar Internacional a una altura de 7600 metros que es la indicada como "cruising altitude" del 1900C. La primera gráfica a analizar es la Figura (2.38) en donde, para un alpha unidad, se representa el empuje específico con respecto a la relación de compresión con θ constante. Se aprecia que las curvas son contínuas y que, para πc bajos presentan un empuje específico relativamente alto, llegando a su máximo en una relación de compresión de 7 u 8 aproximadamente. A partir de aquí, todas las curvas comienzan a decrecer suavemente, debido a que para relaciones de compresión mayores comienza el sistema a no ser tan eficiente y disminuir empuje específico y aumentar consumo. También, a mayor valor de el parámetro θ , mayor es el empuje específico debido a que, si aumentamos la temperatura a la salida de la cámara de combustión, obtendremos más potencia en la turbina y por tanto mayor tracción en la hélice. A continuación, en la Figura (2.39) podemos observar el grafico de TSFC con respecto a πc con los mismos parámetros constantes. En este caso se puede ver un comportamiento inverso al empuje específico; el TSFC decrece de un valor inicial hasta un mínimo que se alcalnza alrededor de πc=11 para después aumentar. De manera contraria al caso del empuje, cuando aumentamos θ disminuye el TSFC ya que en Ee ocurre el efecto inverso y este se encuentra en el denominador de la expresión en cuestión. También se observa que, a mayor θ , es decir, a mayor temperatura en cámara de combustión, más lento es el crecimiento del TSFC una vez alcanzado el máximo. La Figura (2.40) representa el empuje específico frente a θ teniendo como parámetros constantes α=1 y πc . En esta gráfica se puede observar cortes entre las curvas de diversos valores de la relación de compresión, lo que se debe a comportamientos distintos según la elección del parámetro. Para una relación de πc=4 , existe una relación lineal creciente entre Ee y θ , de forma que cuando la temperatura en cámara de combustión crece, el empuje específico crece proporcionalmente. Sin embargo, para relaciones de compresión más altas (12, 16, 20...), se observa una tendencia lineal pero con mayor pendiente, de forma que concluimos que, teniendo un gas más comprimido a la entrada de cámara, cuanto mayor es la temperatura, mayor aún será Ee . Este razonamiento es lógico ya que, un flujo de aire más denso o con mayor presión en cámara de combustión, tiene una zona de inflamabilidad más grande, como se observa en la Figura (2.42). Así, mayor parte del gas entrará en combustión y antes en el tiempo, provocando un incremento del Ee más pronunciado. Pero, cuanto mayor es πcmenor es el corte en ordenadas de la curva, por tanto el mínimo Eese reduce. La Figura (2.41) enfrenta el TSFC con θ , con los demás parámetros constantes. Aquí, de nuevo, se observa dos tendencias dependiendo de la relación de compresión escogida. Para relaciones de compresión bajas ( πc=4,8 ) observamos una débil dependencia del TSFC con θ , en donde, para temperaturas de la cámara de combustión mayores el TSFC y por tanto, el consumo de carburante específico, se reduce de forma muy atenuada. Para relaciones de compresión más elevadas (12,16,20,etc), podemos ver que la tendencia a decrecer 32 Capítulo 2. Modelo en el punto de diseño es mucho más pronunciada al comienzo, cuando los valores de θ son aún pequeños. Esto se debe a que, con unas temperaturas de cámara estándares, si nuestro compresor comprime mucho el gas, no necesitaremos tanto combustible para producir la ignición y la reacción de combustión, como se puede de nuevo observar en el gráfico de limitaciones de inflamabilidad. Figura 2.38 Gráfico empuje específico con relación de compresión. Figura 2.39 Gráfico TSFC con relación e compresión. 2.3 Modelo ’Gas Path’con fluidos reales: programación del ciclo en el punto de diseño 33 Figura 2.40 Gráfico empuje específico con parámetro theta. Figura 2.41 Gráfico TSFC con parámetro theta. En la Figura (2.43) se representa una superficie 3D en donde se dibuja el empuje específico con respecto a los parámetros θ y πc . Se observa que, con relación de compresión constante, al aumentar θ aumenta Ee debido a que la cámara de combustión tendrá mayor temperatura. A su vez, a θ constante, si aumentamos la relación de compresión, el empuje específico aumenta pero en menor medida. Otra gráfica de interés es la Figura (2.44), en donde se vuelve a representar una malla 3D pero esta vez con la variable "z" siendo el TSFC. Observamos que, para valores altos de πc y θ simultáneamente, el TSFC 34 Capítulo 2. Modelo en el punto de diseño Figura 2.42 Límites de inflamabilidad de las mezclas fuel-aire. es bajo. Mientras que si, dejando el valor de θ alto reducimos la compresión, aumentaremos un poco este consumo TSFC debido a que, aún teniendo una alta temperatura en la cámara, el aire que entra en ella para ser parte de la combustión está poco comprimido por lo que la reacción no es muy efectiva. De la otra parte, si dejamos el parámetro πc constante en valores altos, y disminuimos θ nos encontraremos con un aumento mucho más notable del TSFC, y por tanto, del dosado. Esto nos indica que, la temperatura que se alcance en la cámara es un factor muy importante en el consumo de carburante que se necesita. Esto mismo puede observarse en la gráfica (2.45), construida a partir de los datos del ciclo también. En esta gráfica se puede ver la dependencia de la temperatura con el dosado y, por tanto, con el consumo de comsbutible. A mayor temperatura en la cámara, mayor es el corte en la ordenada, significando que a dosado nulo, tendríamos una temperatura mayor. Sin embargo, lo más significativo es que a mayor dosado, mayor temperatura se consigue en la cámara y viceversa, validando entonces el razonamiento anterior de que, a altos dosados se obtienen temperaturas altas (θ). 2.3 Modelo ’Gas Path’con fluidos reales: programación del ciclo en el punto de diseño 35 Figura 2.43 Superficie 3D empuje específico. Figura 2.44 Superficie 3D TSFC. 36 Capítulo 2. Modelo en el punto de diseño Figura 2.45 Gráfico temperatura de la cámara de combustión con el dosado. Para realizar un estudio más profundo de la correlación de las variables de interés, se dibuja la gráfica representada en la Figura (2.46), en donde se enfrentan los valores de Ee en abcisas y TSFC en ordenadas, para un α=1 . Esta gráfica es clásica en el estudio de ciclo de motores y tiene una forma cualitativa característica en donde, las líneas pintadas en azul son curvas de iso-grado de calentamiento (iso θ ) y las dibujadas en rojo curvas de iso-relación de compresión (iso πc ). La relación de compresión aumenta a medida que vamos descendiendo por las lineas rojas y el grado de calentamiento a medida que avanzamos hacia la derecha, siendo los valores representados: πc= [4,8,12,16,20,24]yθ= [4.5,5,5.5,6], en ese orden. Adentrándonos en el análisis de esta gráfica, en una misma línea de iso-grado de calentamiento (por ejemplo,iso θ =5.5), vemos que hay un máximo de Ee a TSFC medio y que, en ambos extremos, cuando el TSFC aumenta el empuje específico disminuye, como nos indica la expresión del mismo: TSFC= f∗(1−ε) Ee . Recorriendo las lineas de iso-relación de compresión (iso πc ), vemos que al aumentar Ee , disminuye suavemente el TSFC, pendiente que es cada vez más acusada para relaciones de compresión mayores. Por último, podemos resaltar que este tipo de gráficas son útiles para la elección de un punto de diseño según la aeronave bajo estudio. La correlación con los distintos tipos de aviones se puede ver claramente en la gráfica cualitativa representada en la Figura (2.47). Para aviones "business jet" tomamos TSFC altos y Ee bajos ya que, nos importa tener una disponibilidad para volar para el usuario y no nos preocupa tanto el consumo que este pueda tener o el empuje específico ya que no necesita volar a altas velocidades. Si nos movemos hacia la derecha, a TSFC y Ee altos simultaneamente, estaremos modelando un avión de uso militar en donde, lejos de ser el consumo lo que nos importe, necesitamos un gran empuje para las maniobras y el combate. Bajando por la curva como indica la flecha, nos encontraremos con aviones de aplicación civil o de "long-range" en donde se necesita, a parte de un empuje específico alto, una minimización del TSFC debido a la importancia de los costes en el funcionamiento de las compañías aéreas. Dicho esto, procedemos a la selección de un punto de diseño en esta gráfica que sirva como referencia para el escalado de los mapa de turbomáquina. Si nuestra decisión fuera tomar el Ee maximo, estaríamos forzados a escoger un punto de tangencia vertical en el mapa, mientras que si necesitáramos minimizar el TSFC se estudiarían los puntos de tangencia horizontal con estas curvas. En nuestro caso, escogemos un punto de operación según el vuelo más común del avión Beechcraft 1900C; el crucero. Los valores buscados se encuentran en un punto medio entre los extremos de aviones antes mencionados. El avión necesitará un TSFC suficientemente bajo, pero lógico según sus características y prestaciones, complementando con un empuje específico medianamente alto. Este punto es escogido, como se muestra en las Figuras (2.48, 2.49), en α=1 , 2.3 Modelo ’Gas Path’con fluidos reales: programación del ciclo en el punto de diseño 37 θ=5.5 y πc=12 , los cuales son valores congruentes con el Manual de Entrenamiento del PT6A de nuestro tipo de motor. Para ello tendremos en el punto de operación un Ee=1266.96 N kg/s , TSFC =60.8091 kg kN·h y un cE=323.438.695 g KW·h . Este punto de diseño se utilizará en los pasos posteriores y se comprobará su validez debido a la lógica de los resultados posteriores. Figura 2.46 Curvas isoπe isoθ. Figura 2.47 Curvas isoπe isoθcualitativas genéricas. 38 Capítulo 2. Modelo en el punto de diseño Figura 2.48 Curvas isoπe isoθ. Figura 2.49 Curvas isoπe isoθrespecto al consumo específico. 2.4 Mapa de la hélice y cálculo de los parámetros N y D Para completar el estudio en punto de diseño, construiremos el mapa de la hélice, al igual que hemos procedido con cada una de las turbomáquinas del motor. Esto se hará para calcular tanto el número de revoluciones (N) como el diámetro (D) y comprobar, por tanto, su validez. De esta manera se verificará que el escalado y la selección de datos que hemos hecho previamente en nuestro motor, son lógicos para el escalado de la hélice 2.4 Mapa de la hélice y cálculo de los parámetros N y D 39 que utiliza. Para proceder, tendremos que digitalizar el mapa de referencia genérico que se puede consultar en la Figura (2.50). El mapa digitalizado se ve en la Figura (2.51) Figura 2.50 Mapa de la hélice de referencia. A continuación se explicitan las fórmulas para el cálculo de dos parámetros de la hélice: el coeficiente de potencia (Cp) y el parámetro de avance (J) 1 . El primero es un coeficiente adimensional relacionado con la potencia necesaria para mover la hélice. El segundo es el paso efectivo adimensionalizado con el diámetro de la hélice, siendo el paso efectivo la velocidad entre el número de revoluciones por unidad de tiempo. Cp =HP 2000σ(N/1000)3(D/10)5(2.86) J=V0 ND (2.87) El rendimiento de la hélice se puede expresar en términos de estos parámetros o con la expresión a continuación: ηprop =E·V0 HP (2.88) El procedimiento que seguiremos para calcular N y D será el que se especifica a continuación. En primer lugar, se calculan los resultados del ciclo en punto de diseño, los cuales ya tenemos, además de conocer la velocidad de vuelo ( V0 ) y la potencia a la entrada de la hélice (HP). Fijamos un punto en la gráfica, ya digitalizada, el cual coincida con la zona del rendimiento ηproph escogido en el ciclo y se calculan los valores de CpPD,JPD en el punto de diseño. Teniendo los valores de V0,JPD se puede obtener N·D en el punto de 1HP(Horsepower), N(rpm), D(ft), V0(velocidad de vuelo en ft/min), σ=δ θ,δ=P Pstd ,θ=T Tstd ,Pstd =100KPa,Tstd =288K 46 Capítulo 3. Modelo fuera de diseño magnitud. Una vez calculados estos factores en una función de Matlab, que se encuentra también en el apéndice de códigos, tendremos que multiplicar cada tabla de valores que conforman los gráficos generalizados por su correspondiente factor de escalado. Como excepción, para el caso de PR tendremos que restar la unidad, multiplicar el factor y volver a sumar la unidad. Una vez hecho este paso, tendremos los mapas escalados de cada turbomáquina y, por tanto, su ley de comportamiento. En algunos casos, las tablas han sufrido retoques debido a que presentaban picos cuando se graficaban y, por tanto, estos se suavizaron para asemejarse más fielmente al comportamiento real del escalonamiento. Los mapas de operación escalados de cada turbomáquina se exponen en las Figuras: 3.4, 3.5, 3.6, 3.7, 3.8, 3.9, 3.10, 3.11,3.12. Estas leyes de comportamiento serán luego "llamadas" por el código del ciclo fuera de diseño, para poder resolver nuevas ecuaciones que se especificarán en la sección posterior. Figura 3.4 Mapa escalado del primer escalonamiento de axial. 3.2 Escalado de los mapas de las turbomáquinas 47 Figura 3.5 Mapa escalado del segundo escalonamiento de axial. Figura 3.6 Mapa escalado del tercer escalonamiento de axial. 48 Capítulo 3. Modelo fuera de diseño Figura 3.7 Mapa escalado del cuarto escalonamiento de axial. 3.2 Escalado de los mapas de las turbomáquinas 49 Figura 3.8 Mapa escalado del compresor centrífugo: gasto corregido. 50 Capítulo 3. Modelo fuera de diseño Figura 3.9 Mapa escalado del compresor centrífugo: rendimiento. 3.2 Escalado de los mapas de las turbomáquinas 51 Figura 3.10 Mapa escalado de la turbina de alta. Figura 3.11 Mapa escalado del primer escalonamiento de turbina de baja. 52 Capítulo 3. Modelo fuera de diseño Figura 3.12 Mapa escalado del segundo escalonamiento de turbina de baja. Como se observa de las figuras, los gráficos tienen señalados distintos puntos de colores que se refieren a distintas condiciones operativas. Estas se identificarán y explicarán en la siguiente sección. 3.3 Resolución del ciclo fuera de diseño La resolución del ciclo fuera de diseño se consigue, en nuestro caso, a través de un sistema de 60 ecuaciones con 60 incógnitas que será explicitado a continuación 12 3 4 Las entradas que tenemos, en este caso, son distintas que en el ciclo anterior. Le daremos al programa como entrada de datos la condicón de vuelo ( P0,T0 , M0 ) y la temperatura T03 de salida de cámara de combustión. Como comentado antes, los parámetros ( θ,πc ) ya no son entrada si no solución del ciclo y, por tanto, la variable que podemos variar a la entrada es esta temperatura de combustión la cual es intercambiable con otras magnitudes como: Ee,˙ m f o,Ngeneradorgas. El esquema del motor que se tratará en este apartado será parecido al anterior pero considerando todos los escalonamientos quue componen los compresores y las turbinas. De esta forma, nos regiremos por el dibujo esquemático presentado en la Figura(3.13). 1 Las funciones que se explicitan como fson códigos en Matlab de matrices de datos correspondientes a los gráficos escalados que, para una entrada de Nc,βdan el valor pedido de mc, pr o η 2 La letra X indica el valor de las fracciones de todos los componentes de productos o reactivos correspondientes en cada ecuación, siendo rpara reactivo y ppara producto. Ra y Rg indican los valores de la constante de gases para el aire antes y después de la combustión. La letra hindica una función que calcula la entalpía correspondiente a una temperatura y composición indicadas. 3 Cuando se escriben las presiones de remanso P0i para i=1,21,22,23,etc, estas se calculan como producto de fracciones de manera que, por ejemplo,P023 =P00 ·πtb ··pr2 4La letra γrepresenta una función que calcula el índice gamma para una temperatura y composición dadas 3.3 Resolución del ciclo fuera de diseño 53 Figura 3.13 Diseño del motor para el ciclo off design. Difusor de entrada (etapa 00-01) T00 =T01 (3.23) Escalonamiento de axial 1 (etapa 01-021) mc1=f(βcomp1,Nccomp1)(3.24) pr1=f(βcomp1,Nccomp1)(3.25) η1=f(βcomp1,Nccomp1)(3.26) φ(T021s,Xr) = φ(T01,Xr)+Ra·log(pr1)(3.27) η1=h(T021s,Xr)−h(T01,Xr) h(T021,Xr)−h(T01,Xr)(3.28) mc1=˙mqT01 Tre f qP01 Pre f pRaγre f pRre f γ01,Xr (3.29) Escalonamiento de axial 2 (etapa 021-022) mc2=f(βcomp2,Nccomp2)(3.30) pr2=f(βcomp2,Nccomp2)(3.31) 54 Capítulo 3. Modelo fuera de diseño η2=f(βcomp2,Nccomp2)(3.32) φ(T022s,Xr) = φ(T021,Xr)+Ra·log(pr2)(3.33) η2=h(T022s,Xr)−h(T021,Xr) h(T022,Xr)−h(T021,Xr)(3.34) mc2=˙mqT021 Tre f qP021 Pre f pRaγre f pRre f γ021,Xr (3.35) Escalonamiento de axial 3 (etapa 022-023) mc3=f(βcomp3,Nccomp3)(3.36) pr3=f(βcomp3,Nccomp3)(3.37) η3=f(βcomp3,Nccomp3)(3.38) φ(T023s,Xr) = φ(T022,Xr)+Ra·log(pr3)(3.39) η3=h(T023s,Xr)−h(T022,Xr) h(T023,Xr)−h(T022,Xr)(3.40) mc1=˙mqT022 Tre f qP022 Pre f pRaγre f pRre f γ022,Xr (3.41) Escalonamiento de axial 4 (etapa 023-024) mc4=f(βcomp4,Nccomp4)(3.42) pr4=f(βcomp4,Nccomp4)(3.43) η4=f(βcomp4,Nccomp4)(3.44) φ(T024s,Xr) = φ(T023,Xr)+Ra·log(pr4)(3.45) η4=h(T024s,Xr)−h(T023,Xr) h(T024,Xr)−h(T023,Xr)(3.46) mc4=˙mqT023 Tre f qP023 Pre f pRaγre f pRre f γ023,Xr (3.47) Escalonamiento de centrífugo (etapa 024-025) mc5=f(βcomp5,Nccomp5)(3.48) pr5=f(βcomp5,Nccomp5)(3.49) η5=f(βcomp5,Nccomp5)(3.50) 3.3 Resolución del ciclo fuera de diseño 55 φ(T025s,Xr) = φ(T024,Xr)+Ra·log(pr5)(3.51) η5=h(T025s,Xr)−h(T024,Xr) h(T025,Xr)−h(T024,Xr)(3.52) mc5=˙mqT024 Tre f qP024 Pre f pRaγre f pRre f γ024,Xr (3.53) Cámara de combustión (etapa 025-03) ˙m·(1−ε 100)·h(T025,Xr)+ ηcc ·H p ·˙mf= ( ˙m(1−ε 100)+ ˙mf)·h(T03,Xp)(3.54) Turbina de alta (etapa 03-04) 5 mc6=f(βcomp6,Nccomp6)(3.55) pr6=f(βcomp6,Nccomp6)(3.56) η6=f(βcomp6,Nccomp6)(3.57) φ(T04s,Xp) = φ(T03,Xp) +Rg·log(pr6)(3.58) η6=h(T03,Xp)−h(T04,Xp) h(T03,Xp)−h(T04s,Xp)(3.59) mc6=˙mtqT03 Tre f qP03 Pre f pRgγre f »Rre f γ03,Xp (3.60) Compatibilidades compresor-turbina de alta ˙m(1−ε 100)+ ˙mf(h(T03,X p)−h(T04,Xp)·ηmec =˙m(h(T025,Xr)−h(T01,Xr)) (3.61) Nre f 6·Nc6T03 Tre f sγ(T04,Xp)·Rg γre f Rre f =Nre f 1·Nc1T01 Tre f sγ(T01,Xr)·Ra γre f Rre f (3.62) Nre f 6·Nc6T03 Tre f sγ(T04,Xp)·Rg γre f Rre f =Nre f 2·Nc2T021 Tre f sγ(T021,Xr)·Ra γre f Rre f (3.63) Nre f 6·Nc6T03 Tre f sγ(T04,Xp)·Rg γre f Rre f =Nre f 3·Nc3T022 Tre f sγ(T022,Xr)·Ra γre f Rre f (3.64) Nre f 6·Nc6T03 Tre f sγ(T04,Xp)·Rg γre f Rre f =Nre f 4·Nc4T023 Tre f sγ(T023,Xr)·Ra γre f Rre f (3.65) Nre f 6·Nc6T03 Tre f sγ(T04,Xp)·Rg γre f Rre f =Nre f 5·Nc5T024 Tre f sγ(T024,Xr)·Ra γre f Rre f (3.66) 5˙mtes el gasto másico que llega a la turbina que se puede ver en las ecuaciones de compatibilidad compresor-turbina 62 Capítulo 3. Modelo fuera de diseño Figura 3.20 Mapa del primer escalonamiento de turbina de baja con 4 condiciones de vuelo. Figura 3.21 Mapa del segundo escalonamiento de turbina de baja con 4 condiciones de vuelo. En general se observa, que todas las condiciones de vuelo planteadas están dentro del rango esperable, sin superar la línea de seguridad o margen por efecto surge. Se observa que en la turbina de alta los puntos se 3.3 Resolución del ciclo fuera de diseño 63 solapan, entendiendo que las condiciones de vuelo no afectarán de manera notable a la relación de expansión o al gasto en la turbina de alta. En la turbina de baja, sin embargo, se observan cambios mayores, por ejemplo, para el punto azul que sería el vuelo al nivel del mar, PR disminuye, es decir la presión de remanso, después de pasar por un escalonamiento de turbina, no se reduce tanto. Finalmente, he de comentar que, para evitar problemas en Matlab con respecto a ordenes de magnitud muy distintos, los cuales provocan a veces imposibilidad de resolver un sistema de muchas ecuaciones, hemos insertado en el código algunos factores de amplificación y reducción, para que algunas variables se calcularan en unas magnitudes que hicieran que su resultado estuviera en un orden de magnitud próximo al de las demás. 4 Condiciones de degradación Con el objetivo de estudiar cómo se degrada el motor, intentar localizar de dónde viene el fallo interno y analizar la eficacia y eficiencia de la herramienta de diagnóstico utilizada, tendremos que hacer un estudio previo de las condiciones de degradación. Debido a la falta de un motor real con una degradación o fallo real, tendremos que, a nuestro modelo digital del motor, introducirle unos errores que son conocidos a priori por nosotros. En condiciones reales la degradación de un motor no es normalmente conocida y se intenta estimar por medio de pruebas experimentales o numéricas. Cómo uno de los objetivos del proyecto es el análisis de la eficacia de nuestro método para identificar la degradación de una parte del motor de forma numérica, necesitaremos conocer el error introducido en nuestro modelo en Matlab para así poder compararlo con el resultado del optimizador y evaluar la bondad del método bajo estudio. Sin embargo, como es lógico, esto solo será un dato para comparar los resultados de la herramienta de diagnóstico, ya que a esta le vendrán dados unos datos que simulen unos reales, sin información previa de cuál es la degradación que, a priori, hemos introducido en el motor. Se podrá ver, a lo largo de los próximos capítulos, que la casuística de fallos y degradaciones que pueden acaecer en un motor son irrisorias, casi incuantificables. Para hacer un estudio completo habría que introducir degradaciones en todos los escalonamientos de turbomáquinas, en diferentes órdenes y cuantías, con porcentajes de degradación distintas y para las tres magnitudes consideradas en las leyes de comportamiento de las turbomáquinas: mc,η,pr . Esto supone una cuantía de información enorme y, debido al limitado alcance del proyecto no es posible abarcar el estudio completo. Por ello, como explicaremos a continuación, se tratará de hacer un estudio completo de la casuística de fallo para el primer escalonamiento de axial, el cual es uno de los lugares con más probabilidad de fallo en un motor, para distintos grados de degradación en las tres variables a considerar. También se estudiarán casos aislados para poder evaluar la variedad de datos que puede tratar la herramienta de diagnóstico. 4.1 Procedimiento sistemático para definir las condiciones de degradación 4.1.1 Magnitudes con sensor en el motor PT6A Para poderle introducir a la herramienta de diagnóstico datos cercanos a la realidad, tendremos que definir las magnitudes que son medidas por los sensores del motor. Cómo no es completamente pública toda la información interna de los motores de nuestra gama se ha intentado estimar qué magnitudes son las que se tiene medida conocida, a partir del estudio de los sensores que hicimos previamente. Se estima que, para un motor de nuestras características, las magnitudes medidas en el motor que serían de interés para nuestro estudio y que corresponderían al estudio de sensores previo serían: las revoluciones por minuto del generador de gas ( Ncent ), temperatura ITT ( T04 ), temperatura de salida de la turbina de potencia ( T05 ), par de la hélice y, por tanto, potencia de la hélice, gasto másico de combustible ( ˙mf ), presión y temperatura a la salida del compresor centrífugo y entrada en cámara de combustión ( T025 , P025 ), la velocidad corregida del compresor centrífugo que es una medida indirecta de la presión P024 y la temperatura T024 ( Nc5 ) 65 66 Capítulo 4. Condiciones de degradación 4.1.2 Condiciones de vuelo A continuación tendremos que definir las condiciones de vuelo en donde se hará el estudio. Nosotros, en general para el estudio del primer escalonamiento, vamos a imponer condiciones de cota nominal y velocidad nominal aunque se pueden definir otras tres condiciones de vuelo representativas en el estudio que luego se incluirán en el apartado de casos aislados. Estas cuatro condiciones se especifican a continuación: •Cota nominal y velocidad nominal •Cota nominal y 70% velocidad nominal •70% cota nominal y velocidad nominal •70% cota nominal y 70% velocidad nominal En el primer estudio del escalonamiento del primer axial, el ciclo off-design será utilizado para, no estudiar la eficacia en otras condiciones de vuelo que no sean el punto de diseño, sino estudiar puramente la degradación de los componentes en un vuelo normal. 4.1.3 Matriz de degradación Para poder plantear de manera ordenada y efectiva un estudio estadístico completo de la degradación del escalonamiento„ debemos definir lo que se denomina una matriz de degradación, que pone en evidencia todas las posibles degradaciones que se consideran en el estudio. Al tener 3 variables que se pueden degradar ( mc , pr , η ), vamos a considerar, para tener un número razonable de datos, tres grados de degradación ( 1% , 3%,5%). De esta forma, la matriz o tabla de degradación, que ilustra los casos bajo estudio sería: Degradación (%) Gasto corregido (mc) 1 3 5 Rendimiento (pr) 1 3 5 Relación de presión (η) 1 3 5 De esta forma, definimos una terna de degradación a un trío de valores de degradación para las tres variables consideradas. Para un estudio estadístico completo del primer escalonamiento se necesitarán entonces, por combinatoria, 33=27 ternas de degradación. Para cada terna de degradación, se genera un nuevo mapa de compresor del escalonamiento. Esto se obtiene degradando el gasto corregido, el rendimiento y la relación de presión, multiplicando cada tabla por un factor de reducción. Por ejemplo, para la terna (1, 3, 5), se multiplicará la matriz de gasto corregido por 0.99, la de relación de presión por 0.97 y la de rendimiento por 0.95. Así, el programa resuelve el ciclo fuera de diseño para la condición de vuelo considerada y con el error considerado. 4.1.4 Adición de ruido Como hemos visto, los sensores que tenemos en el motor son limitados, y las magnitudes que miden son las que nosotros, en una situación real, tendríamos a nuestra disposición. Estas magnitudes, cuando el motor está degradado, se puede simular con las ternas que hemos definido anteriormente, siendo las variables medidas por los sensores las salidas del ciclo fuera de diseño una vez introducido este error. Sin embargo, en la realidad, los sensores tienen errores de medida, y en específico, es importante considerar el ruido de señal que se acopla al sensor cuando este mide y que, por tanto, altera el resultado que nosotros vemos. Así, a cada magnitud generada por el modelo fuera de diseño hay que añadirle ruido de señal que se contabiliza como una magnitud aleatoria de valor ±Rx% del fondo escala del sensor. Este valor Rx depende de la calidad del sensor y suele ser un valor comprendido entre el 0.5% y el 1% . En nuestro caso, al estudiar y probar ambos errores en el código, encontramos que la solución total del método es muy tosca con un error del 1% y, por tanto, consideramos un error del 0.5% , valor completamente válido para sensores internos de un motor en aeronáutica. 4.1 Procedimiento sistemático para definir las condiciones de degradación 67 Por otro lado, el fondo escala del sensor lo estableceremos a un 20% de la magnitud máxima esperada. Por ejemplo, si esperamos que en el compresor centrífugo la máxima presión alcanzada en despegue sea de 12 bar (situación fuera del punto de diseño), el sensor correspondiente, para dar un margen a la medida de esta magnitud, se seleccionará con un fondo escala de Px=12 ·1.2=14.4bar. Para la generación de la señal de ruido, utilizaremos la función randi en Matlab, de forma que, para que su valor oscile entre 1 y -1, plantearemos el siguiente código: r(i) = 0.1∗randi([−10,10],1,1)(4.1) De esta forma, para definir la magnitud que a nosotros nos llega del sensor, tendremos que sumarle a la medida obtenida del modelo OD (PP), el error r(i) multiplicado por su valor Rx y por el valor del fondo escala del sensor, de manera que el cálculo quedaría: Magnitud sensor =PP +0.1·randi([−10,10],1,1)·Px·Rx 100 (4.2) Para poder comprobar que el ciclo OD funcionaba de manera correcta con estas degradaciones, hemos probado a comparar la diferencia entre las magnitudes de salida del ciclo con y sin degradación. Los resultados son satisfactorios ya que, para una terna aleatoria de degradación (1, 3, 5) obtenemos unos resultados degradados y sin degradar que son diferentes pero no excesivamente debido a que el fallo no es demasiado grande.Los resultados se pueden ver en la tabla (4.1). Como se puede observar, el gasto de combustible, es mayor cuando tenemos el primer escalonamiento degradado, ya que el motor no está funcionando al 100% de sus capacidades. Tabla 4.1 Comparativa de valores degradados y sin degradar. Terna de degradación (1,3,5) Sin degradación Ncent =39160rpm Ncent =39083rpm T04 =1016.8K T04 =1016.8K Par hélice= 3080.4 N·mPar hélice= 3140.0 N·m T05 =779.8649K T05 =778.1512K ˙mf=165.1645kg/h˙mf=167.1693kg/h T025=582.7489 K T025=582.7246 K P025=479.7564 KPa P025=485.4881 KPa Nc5=0.9749 Nc5=0.9722 He de comentar que esta fase de estudio del primer escalonamiento tiene en cuenta el compresor completo, es decir, cuando se da la terna de degradación, se da al programa un vector con los tres primeros valores de la terna y los 12 restantes nulos (ternas nulas de el escalonamiento 2,3 4 y centrífugo), ya que, luego, para el proceso de optimización se tendrá en cuenta el efecto del compresor como conjunto. 4.1.5 Generación de datos de vuelo con motor el motor degradado A la hora de crear una muestra de datos que escenifique un estudio estadístico completo y fiable, tendremos que generar, no solo un paquete de datos de los sensores para un ruido en un instante de tiempo, sino una matriz de datos que genere los datos recogidos en un tiempo considerable, es decir, los datos de vuelo. Nosotros, para definir el estudio, consideraremos 100 instantes de tiempo, con 100 señales de ruido distintas, como ocurriría en la realidad. Bien, para ello, cada condición de vuelo (que en nuestro caso es considerada solo la nominal) y para cada terna de degradación (27 ternas) se genera una matriz S(X,Z) donde X es el número de datos del vector aleatorio de ruido y Z las magnitudes generadas por los sensores (8 magnitudes en nuestro caso).Así, para generar datos de entrada para la herramienta de diagnóstico, habrá que crear 27 matrices de dimensión 100x8. Todo este proceso está implementado en un código en Matla con la ayuda de un bucle for 68 Capítulo 4. Condiciones de degradación 4.1.6 Herramientas para la estimación de la degradación Una vez definida la entrada, o lo que es lo mismo, las magnitudes medidas por los sensores, tendremos que escoger una metodología o herramienta de estimación de degradación y aplicar la metodología oportuna. Sabemos que para este problema y otros de la misma índole, hay muchas metodologías disponibles, las cuales se comentarán en el capítulo posterior. El algoritmo que nosotros utilizaremos está ilustrado en el esquema de la Figura (4.1), que se basa en el uso de funciones de optimización de gradientes conjugados. Figura 4.1 Metodología de la herramienta de optimización. Se puede ver que estos algoritmos o herramientas que proporciona Matlab como funciones, se basan en un proceso iterativo simple. Una vez decidida la terna de degradación, y tener la matriz correspondiente S(X,Z) para ella, esta se mete como entrada en el modelo. A su vez, la herramienta tiene como incógnitas un vector t= ( t1 , t2 , t3 , t4 , t5 , t6 , t7 , t8 , t9 , t10 , t11 , t12 , t13 , t14 , t15 ) que es el vector que identificaría las ternas de error del compresor completo, el cual, para nuestro caso ejemplo, tendría la solución exacta de (0.99,0.97,0.95,1,1,1,1,1,1,1,1,1,1,1,1). Para poder resolver, se introduce este vector incógnita como un error en el modelo OD, teniendo un iterante inicial para comenzar que será (1,1,1,1,1,1,1,1,1,1,1,1,1,1,1). El modelo OD calculará una matriz de datos M(X,Z) con respecto al error inicial, que irá cambiando debido al proceso iterativo. La función optimizadora, irá cambiando los valores del vector incógnita de forma que, en cada iteración, evaluará la diferencia entre la matriz entrada de los sensores y la evaluada con el vector t, fijándose si esta diferencia es menor que una cierta tolerancia predefinida. Si la diferencia es mayor que la tolerancia, se cambiará oportunamente el valor de los valores del vector t, se volverá a evaluar el modelo OD con estos errores, y se volverá a comparar esta matriz con la de los sensores. Esta metodología se explicará en más profundidad en el siguiente capítulo. Para volver a evitar el problema de los ordenes de magnitud que comentamos en la resolución del ciclo fuera de diseño, se multiplicarán los valores de algunas magnitudes medidas por sensores por factores que hagan que estos estén en unidades más adecuadas. Por ejemplo: el par de la hélice vendrá dividido por 10 midiendose así en dN ·m , el valor de Nc5 se multiplicará por 1000, el Ncent se cambiará de rpm a rad/s. El objetivo de todo esto es poder tener todas estas magnitudes en un orden de 102aproximadamente. 5 Herramienta de diagnóstico y resultados Para identificar dónde está el fallo interno en el motor y poder discernir la cuantía de degradación del mismo, hay varias herramientas de diagnóstico que pueden ser usadas. La que nosotros hemos utilizado es una basada en optimización, sin embargo, debido a la limitación del alcance del proyecto, no se han podido evaluar otras como son las estrategias de diagnóstico basadas en GPA (Grade Point Average) y en redes neuronales. El objetivo de este proyecto no es solo crear el código para la utilización de la herramienta de diagnóstico del motor sino también evaluar la bondad y eficacia de la misma. Es posible que un enfoque basado en redes neuronales hubiera dado resultados más precisos pero, estos conllevan un tiempo, coste computacional y cantidad de estudio exhaustivo que estaba fuera de nuestros límites. La dirección tomada ha sido en base a algoritmos de optimización, los cuales, como veremos son de concepto sencillo pero útil. Las dos funciones utilizadas en Matlab han sido fmincon y lsqnonlin , las cuales están basadas en el concepto de algoritmo de optimización CG (Conjugate Gradient), es decir, de gradientes conjugados. Este tipo de optimizadores son de carácter local, lo que significa que, no nos dan la seguridad de encontrar un mínimo absoluto de una función, sino que, el mínimo que encuentran depende de las condiciones iniciales de entrada. Estos funcionan mejor en entornos donde la función objetivo que se quiere minimizar, en nuestro caso, la diferencia entre la matriz de modelo OD y la de sensores, no cambie mucho, es decir no tenga muchos mínimos y máximos locales. Desgraciadamente, esta información no la conocemos ya que la función objetivo que se evalúa aquí, es una resta de matrices construidas a partir de valores aleatorios. Además, el problema bajo estudio tiene una incógnita de magnitud 15, es decir, el optimizador, solo para resolver este primer escalonamiento de axial, tiene que lidiar con 15 parámetros libres. Esto, en principio, solo complica el problema, pero no hace que nos salgamos de los límites de validez de este tipo de algoritmos ya que el mayor número de parámetros libres que estos usualmente suelen permitir son 25. El funcionamiento de estos métodos es sencillo, y se basan en algo similar al algoritmo de Newton para encontrar ceros en una función. Supongamos que tenemos una función objetivo, de la cual queremos encontrar el mínimo (es decir, extrapolado a nuestro caso, queremos encontrar la solución del vector t que me minimice la resta de las dos matrices, problema que tiene mínimo absoluto y solución exacta), y unos cuantos grados de libertad que podemos variar para este cometido. Bien, lo primero que se hace es evaluar las derivadas direccionales de esta función para cada uno de los parámetros libres y ver cuál de ellas es mayor (esto se hace evaluando una pequeña perturbación de la función con respecto a un pequeño cambio del parámetro, de manera numérica). Una vez tengamos cuál es mayor,el algoritmo se tomará la libertad de variar la solución de manera que se varíe ese parámetro en la dirección que haga que la función disminuya (esto podrá ser que el parámetro aumente en cada iteración o disminuya, dependiendo de su relación con la función objetivo). Así, el algoritmo continuará hasta el punto que vea que la función ya no disminuye sino que, en cambio, esta comienza a aumentar. Esto es indicación de un mínimo local, y por tanto, la iteración cesa, y el programa devuelve la solución obtenida. Este concepto se ilustra en la Figura (5.1). 69 70 Capítulo 5. Herramienta de diagnóstico y resultados Figura 5.1 Cálculo de un mínimo con algoritmo de optimización. Un punto en contra es que, dependiendo de dónde comencemos el estudio, es decir, cuál sea nuestra condición inicial, dependiendo de cómo de cerca esta este del mínimo absoluto y de la forma de la función, el resultado que encontraremos será el mínimo absoluto o solo uno local. Esto también puede visualizarse en la Figura (5.2). Es por ello que, como comentaremos posteriormente, nuestra herramienta siempre funcionará mejor si le damos como entrada un iterante inicial cercano a la solución exacta que conocemos a priori. Obviamente, en la realidad esto podría implementarse dependiendo de las condiciones de estudio en las que estemos. Si no sabemos absolutamente nada del fallo, pondremos como iterante inicial el vector t=(1,1,1,1,1,1,1,1,1,1,1,1,1,1,1) pero si, por el contrario, sabemos, ya sea por experiencia, manutención o exploración física del motor, que la degradación debe ser en un escalonamiento en específico, esto se puede plasmar en el iterante inicial dado a nuestra herramienta. Como, en principio, desconocemos esta situación real, nosotros impondremos un iterante inicial representativo de una situación de fallo completamente desconocida. Figura 5.2 Problema del cálculo de mínimos con optimización basada en gradientes conjugados. Es de interés explicar en unas pinceladas, de qué trata cada función en Matlab utilizada, las cuales implementan un algoritmo de optimización. Comenzamos por la función fmincon, cuyo objetivo es encontrar un mínimo de una función multivariable no lineal restringida. Con restringida se refiere a que, a la función y al proceso de optimización se le pueden imponer restricciones expresadas en inecuaciones. Por ejemplo, si nosotros queremos que una cierta magnitud csea menor que 0, se puede imponer esto como vinculo y, a continuación, proceder con la optimización, siempre teniendo en cuenta esta condición. Debido a que nuestra función es multivariable y, mediante prueba y error, descubrimos que la función no resolvía de manera correcta si le decíamos que nos minimizara la resta de cada fila de la matriz diferencia, impusimos que nos minimizara el error mínimo cuadrático de esta resta (ecuación 5.1). EMC =1 n n ∑ i=1Ymod,i−Yi2(5.1) 71 Si este error, por ejemplo, hubiéramos querido que valiera más que un cierto número, por razones externas, se puede imponer en la llamada a la función, la cual se formula así: x=f mincon(f un,x0,A,b,Aeq,beq,lb,ub,nonlcon,options) , de forma que a la entrada de la misma se llama a la función multivariable de la que se requiere el mínimo, al iterante inicial (x0), a las desigualdades que permitan restringir la función (A,B), a las igualdades que permiten restringir (Aeq, Beq) los límites superiores (ub) e inferiores (lb) de la solución, a la limitación de la minimización a las desigualdades/igualdades no lineales con un parámetro (c o ceq) y a la adición de opciones. La segunda función se llama lsqnonlin, la cual tiene el objetivo de resolver problemas de mínimos cuadrados no lineales. Es decir, esta función nos hace el mínimo de error mínimo cuadrático comentado antes, de manera automática. Una diferencia notable entre esta función y la anterior es que esta acepta que se le entregue como función la resta de ambas matrices, como un vector, es decir, no hace falta entregarle el error mínimo cuadrático ya calculado para que este sea un escalar si no que puede entregársele un vector. Esto, como veremos, hace que esta función tenga resultados más satisfactorios en comparativa. 5.0.1 Comparación de resultados de "fmincon" y "lsqnonlin" para casos dispersos Con objetivo de optar por una de las dos funciones de Matlab para realizar un estudio estadístico de los resultados, hemos hecho una comparativa entre los resultados que dan una y otra función para casos de degradación y opciones distintas. El procedimiento ha sido de carácter exploratorio, es decir, definiendo una serie de tolerancias, opciones y ternas de degradación diversas, se puede analizar a simple vista cuál de las dos opciones nos resulta más satisfactoria a la hora de encontrar un mínimo del vector t, que se acerque a la solución, que es conocida. Para comenzar, se discernió desde el principio que el error que se utilizaría era del Rx=0.5% , esto es debido a que, cuando probamos una de las funciones para el mismo iterante inicial, terna y opciones, pero cambiábamos el error a un valor del 1% , se observaba que los resultados eran bastante más erróneos, partiendo de la noción de que ya es complicado encontrar una solución precisa. Esto lo podemos ver en la tabla ( ?? ). Se hizo una comparativa de, suponiendo un vector incógnita de solo tres componentes, es decir, solo considerando el primer escalonamiento de axial y no el compresor completo, los resultados para un error Rx=0.5% y Rx=1% . Se puede razonar de los datos de la tabla que, para un error mayor la función tarda menos en finalizar la optimización, por el encontrar un mínimo local antes, pero que la solución que da se parece mucho menos a la solución Xsol que la que tiene un error menor. Esto es un resultado predecible y esperable ya que, cuanto más ruido acepte el sensor, la señal se ve más distorsionada y, por tanto, el optimizador tiene menos información para poder discernir, de entre el ruido, cuál es la solución mejor. Xsol = [0.98,0.99,0.96]Rx=0.5% Rx=1% Tiempo cálculo t=261.37s t =196.1176s Solución X=[0.9856, 0.9529, 1.0274] X=[0.9874, 0.9308, 1.0731] Por otro lado, está la comparativa entre ambos métodos, que ya se hizo para un vector t de 15 incógnitas,es decir, considerando el compresor completo, aunque no las turbinas. En la tabla (5.1), se presentan los valores comparativos entre los resultados dados por "fmincon" y "lsqnonlin" en diferentes situaciones. Debemos matizar el título dela columna "iteraciones permitidas". Al ser un problema multivariable, con ruido y con un iterante inicial que no es igual a la solución deseada, el programa es, en ocasiones, incapaz de encontrar una solución al problema con la tolerancia que se le pide por defecto, en un tiempo razonable. Por ello, para poder tener una solución en todo caso, debemos de imponerle un máximo número de iteraciones. Esto mismo es lo que se le impone a la función "fmincon". Mientras que, a la función "lsqnonlin" lo que se impone es el número máximo de evaluación de la función objetivo. 78 Capítulo 6. Estudio explorativo y estadístico de la capacidad de resolución de la herramienta de diagnóstico (a) a. (b) b. Figura 6.3 Histograma con diferentes anchos para t3 y t4 de terna 335. Se puede observar que, las figuras que posteriormente iremos señalando, están colocadas en pareja. Esto es, para cada ”t” y terna, se representan dos gráficos que se ponen uno al lado del otro para su comparación visual: PDF (Probability Density Function) de los datos a la izquierda y el histograma a la derecha. El histograma, con el ancho de barra comentado, tiene dos gráficas diversas, solapadas entre sí pero pudiendo verse ambas debido a que se ha añadido un efecto de transparencia. Las barras azules constituyen la función de densidad de probabilidad pero representada de forma discreta, es decir, el área de las barras constituye la probabilidad de que un dato salga como solución. La barras amarillas, las más importantes, son aquellas que representan el porcentaje de un dato en el estudio. Es decir, las veces que ha salido un valor dentro del rango de la barra, pero en forma de porcentaje, dividido entre los datos totales. Comenzamos con los resultados de el gasto corregido ( mc ) del primer escalonamiento de compresor, en la Figura (6.4). Se observa que, en la gráfica de la izquierda, la tendencia de datos es similar a una normal o Gaussiana con la excepción de los datos extremos, tanto a izquierda como a derecha. Solo puede observarse el pico central y su posterior caída. Este fenómeno es debido a que, el algoritmo de optimización usado centra sus solución en un cierto dato que cree que es correcto, por ejemplo el 0.982. Tiene lógica que, en todos los casos con ruido aletorio, la herramienta siempre dé valores cercanos o entorno a este resultado, ya que intenta minimizar la función objetivo. El programa, por tanto, es raro que de un valor de, por ejemplo 0.95, que constituiría un extremo de la gráfica. Por ello, los resultados se concentran en el pico de ”máxima probabilidad” y solo encontramos ciertos casos aislados en donde el resultado se encuentra fuera del rango central. Esto se ve en la Figura comentada donde, el gran grueso de datos está en el intervalo 0.98-0.986 y solo hay un valor, de toda la muestra, que daba un resultado de 0.9745. Además, a la derecha, en el histograma, podemos ver que la PDF en azul representa una PDF que ha construido Matlab para hacer ”fitting” con los datos. Sin embargo, como hemos comentado, la información importante, reside en las barras amarillas. Aquí podemos ver que, el rango de datos que más ha salido como solución en el antepenúltimo, teniendo de valor central 0.984 aproximadamente, seguido de valores alrededor de 0.983. Bien, sabiendo que la solución real debería haber sido 0.97, sabemos que el programa no nos ha dado un dato muy equivocado pero a la vez no muy preciso. Esto se debe a lo explicado anteriormente de la no linealidad del sistema, el acoplamiento entre las distintas variables y que el algoritmo suele parar cuando encuentra un mínimo local, no necesariamente un mínimo absoluto. 6.3 Resultados 79 (a) a. (b) b. Figura 6.4 PDF e histograma. Algo parecido a lo explicado ocurre en las gráficas de la Figura (6.5), que representan la solución para la degradación de la relación de presiones, es decir, el parámetro ”t2”. Aquí, las distribuciones PDF están desplazadas a la izquierda, al contrario que las anteriores. Esto quiere decir que la tendencia de los datos es, acumularse en el pico, aunque también se sobreestimen los resultados algunas veces, sin aparecer casos en donde se estimen valores por debajo de los del rango principal. Aquí, puede verse claramente, en las gráficas amarillas, que hay una moda clara, el valor de 0.9722 aproximadamente. Este, como ya veremos en próximas tablas, coincide con la media, lo cual indica un estudio estadístico correcto. Aquí, el valor se acerca muchísimo al real (0.97), pareciendo ser el método en esta terna más efectivo para la degradación en ”pr” que en el gasto corregido. (a) a. (b) b. Figura 6.5 PDF e histograma. La siguiente pareja de Figuras (6.6), nos muestran los resultados para el rendimiento en el primer escalonamiento de axial. En este caso, la PDF tiene un pico mucho más pronunciado, dando el método un valor muy repetidamente, cosa que se ve también en el histograma. El valor, con diferencia más repetido, es el 0.99 aproximadamente. En este caso, el algoritmo no ha hecho un buen trabajo ya que la solución real estaría en 0.95. Es cierto que, ha podido detectar una degradación pero no en la magnitud que deseamos, lo cual es útil pero no lo eficaz que querríamos. En esta terna, parece entonces ser, que el valor de degradación peor estimado, ya sea por las ecuaciones en las que está envuelto o su dependencia con otras magnitudes, es en el rendimiento del escalonamiento. 80 Capítulo 6. Estudio explorativo y estadístico de la capacidad de resolución de la herramienta de diagnóstico (a) a. (b) b. Figura 6.6 PDF e histograma. Para tener una visión más profunda de lo que ocurre en el estudio, he querido representar también en esta sección, los resultados de esta misma terna pero para el segundo escalonamiento (t4,t5,t6), ya que, este está cerca del primero y, por tanto, como veremos, sufre interferencia de la degradación del mismo, creyendo así el optimizador que este segundo escalonamiento también sufre degradación aunque esta, en la realidad sea nula. En la Figura (6.7), se observa que la moda de los datos que nos devuelve el programa es de una degradación de gasto corregido de 0.9825. Este valor, obviamente, no es real ya que debería devolver un 1, pero,se observa, curiosamente, que la distribución es mucho más plana e uniforma. Es decir, no indica un pico claro en donde se encuentra la mayoría de valores, sino que hay muchos casos en donde salen muchas soluciones distintas. Esto, aun no siendo correcto, puede ser una buena señal debido a que, parece que el optimizador intenta, en cada señal de ruido y prueba, minimizar el valor del parámetro t4, sin llegar nunca a un intervalo concreto en donde presente un pico pronunciado, ya que, no existe la degradación. (a) a. (b) b. Figura 6.7 PDF e histograma. Esto mismo ocurre con la Figura (6.8), en donde la distribución es más plana, y nos da un valor de 0.975 aproximadamente. De nuevo, resultado erróneo, pero tendencia correcta ya que la distribución se apoya en el lado derecho, que cubre también el valor 1. Este no es el que más se repite, pero, constituye solución en algunos casos. 6.3 Resultados 81 (a) a. (b) b. Figura 6.8 PDF e histograma. Sin embargo, en el caso del rendimiento en la Figura (6.9) no se puede decir lo mismo. Aquí, la distribución tiene un pico pronunciado, es decir, que para la mayoría de casos, el optimizador ha encontrado una solución parecida. El valor medio de la misma es de 0.991, lo cual es incorrecto, pero es cierto que se acerca mucho a uno, por lo que podríamos decir que es medianamente satisfactorio. (a) a. (b) b. Figura 6.9 PDF e histograma. Como añadido a estos gráficos, nos faltarían las parejas de Figuras que constituyen los parámetros ( t7 , t8 , t9 , t10 , t11 , t12 , t13 , t14 , t15 ), de los escalonamientos tercero, cuarto y centrífugo. Estos datos se encuentran completos en el apéndice para su consulta. Sin embargo, he considerado oportuno, mostrar diversos datos calculados de los 15 parámetros que son representativos del estudio. En primer lugar, observamos la tabla primera de resultados. Aquí se muestra si el optimizador, utilizando los valores medios que devuelve, da un resultado positivo, negativo, falso positivo o falso negativo. Se observa que, para el primer escalonamiento, independientemente de los valores que consiga, los resultados son todos positivos verdaderos. Es decir, habiendo degradación en la primera corona, este la detecta sin problemas. En el segundo escalonamiento podemos observar en todas las columnas un ”3”, lo que indica un falso positivo, es decir, detecta una degradación cuando no la hay. Esto se debe a la complejidad del problema y a la no linealidad del mismo. Se puede ver una mejora en el tercer escalonamiento donde, a parte de falsos positivos, en la relación de presiones nos da un negativo, lo cual es un resultado correcto, no habiendo degradación en este sector. Una mejora parecida de los resultados aparece en el escalonamiento de centrífugo, en donde el programa 82 Capítulo 6. Estudio explorativo y estadístico de la capacidad de resolución de la herramienta de diagnóstico detecta, para la relación de presiones y el rendimiento, un negativo en la degradación, siendo esta información verídica. La tendencia de estos 4 resultados que se muestran en la gráfica, es muy parecida en casi todas las ternas estudiadas. Sin embargo ponemos solo este ejemplo para que pueda verse en general cómo funciona. Como mencionamos antes, los resultados de esta tabla deben de analizarse para el conjunto de los tres valores de cada escalonamiento. Así, veríamos que el primer escalonamiento en su conjunto da un valor positivo, mientras que el segundo daría un falso positivo en su completo, al igual que el cuarto. Considerando el resultado mayoritario en cada pareja de tres, se observa que para el escalonamiento 3 se considera que nuestra herramienta da un falso positivo peo que el escalonamiento 5 es analizado acertadamente, ya que podemos considerar como si nos devolviera un negativo en degradación. Resultados mcprη Escalonamiento 1 1 1 1 Escalonamiento 2 3 3 3 Escalonamiento 3 3 2 3 Escalonamiento 4 3 3 3 Escalonamiento 5 3 2 2 Es de interés también la tabla de resultados de ”media” ya que se nos muestran el valor medio de la solución que da el optimizador. De nuevo, se ve cómo la degradación está relativamente bien detectada en los parámetros ”t1” y ”t2” en la primera corona, mientras que, para el rendimiento el valor de 0.9882 es bastante equivocado. No obstante, pueden compararse también estos valores con los histogramas anteriores, para verificar si la moda coincide con la media, lo cual da fiabilidad al estudio. En todos los casos, la moda y la media residen en el mismo rango de ancho de barra, lo cual es buena señal, significando que el dato más repetido es el que constituye también el valor medio de toda la distribución, es decir, puede ser considerado como la solución, a grandes rasgos, de ese escalonamiento. Además, estudio de los 4 resultados explicado anteriormente está basado en los resultados de media de esta tabla. Resultados t1 t2 t3 t4 t5 t6 t7 t8 t9 t10 t11 t12 t13 t14 t15 Media 0.9827 0.9722 0.9882 0.9823 0.9734 0.9911 0.9800 1.0033 0.9929 0.9532 1.0059 0.9927 0.9827 1.0037 0.9960 En la tabla siguiente se muestran resultados del denominado ”error relativo”. Con esto pretendemos estudiar cuánto es el error que se comete en magnitud de la degradación calculada versus la real. Es decir, simplemente, se ha calculado el error relativo entre el valor medio que devuelve la herramienta, con el valor verdadero. Por ejemplo, en el caso de el parámetro t1, tendríamos que εrel =0.97 −0.9827 0.97 ·100 en porcentaje. Así, el resultado sería del −1.3090% lo cual significa que ha fallado un −1.3090% del valor verdadero y, al tener signo negativo, ha estimado una degradación menor de la que realmente se tiene. Observando todos los resultados de error relativo en conjunto, podemos ver que el mayor error que se comete es de 4.6847% , valor que ya rozaría lo inaceptable pero que no es descabellado para ser el máximo. Resultado importante es también que hay tanto errores positivos como negativos. Este resultado es bueno ya que significa que el optimizador, en el intento de encontrar la solución más óptima, se mueve tanto por encima como por debajo del valor real. Si solo viéramos errores positivos o solo negativos, esto diría que el método tiende a sobreestimar o subestimar datos,lo cual no es positivo para un estudio en donde se requiere ser precisos. También es interesante ver, para un mismo parámetro t, si el error relativo en cada solución (de las 50 que hay) se mueve entre valores positivos y negativos. Esto, aun habiendo sido calculado en Matlab, puede verse con más claridad en el histograma correspondiente al primer parámetro. Se ve que, todos los datos recogidos para la solución ”t1” son mayores que el mínimo valor que se representa, que es 0.97. Es decir, todas las soluciones dan una degradación menor de la que realmente hay, teniendo una tendencia en este caso a subestimar el error, lo cual no es bueno. El vector que se calcula en Matlab del error relativo sería, por tanto, todo resultados de porcentajes negativos de longitud 50 (este se puede consultar en los códigos del proyecto). Bien, aun así, esto no ocurre siempre, por ejemplo, en el parámetro 2, vemos que se calculan datos por debajo de 0.97 y por encima. En este caso el optimizador presentaría valores a izquierda y derecha, no sólo subestimando la degradación sino también sobreestimándola, siendo esto un resultado deseado desde el punto de vista estadístico. 6.3 Resultados 83 Resultados t1 t2 t3 t4 t5 t6 t7 t8 t9 t10 t11 t12 t13 t14 t15 Error relativo -1.3090 -0.2236 -4.0187 1.7669 2.6634 0.8883 2.0043 -0.3266 0.7057 4.6847 -0.5898 0.7317 1.7289 -0.3686 0.3596 Por último se recogen los datos del nivel de confianza en la siguiente tabla. En este caso, no hemos representado una PDF en la cual se representen los anchos de banda ara un cierto nivel de confianza pero si que presentamos cálculos que pueden ser de interés. El nivel de confianza es un dato en estadística que indica la probabilidad máxima con la que podemos asegurar que el parámetro que se estima (el valor de t) se encuentra dentro del intervalo estimado. Es decir, el porcentaje de probabilidad de que nuestro dato esté dentro de un cierto rango de valores en nuestra distribución PDF. Nosotros definimos el nivel de confianza estándar ( 95% ), que se calcula como, el doble del valor de la desviación estándar, calculada en Matlab por ”fitting” de una distribución normal. Esto constituiría un ancho de banda, desde el centro de la distribución PDF hacia los lados extremos de longitud 2σ , a cada lado. La tabla que se presenta, calcula un 0 si el valor verdadero de degradación no se encuentra en el intervalo determinado por el nivel de confianza y un 1 si sí lo está. Se observa que hay muchos casos en donde, como por ejemplo los parámetros t1 y t3, el valor real no se encuentra en nuestro intervalo. Esto es debido al fenómeno del pico que hemos comentado antes. Al dar el optimizador solo valores cercanos dentro de un rango estrecho (formando un pico de distribución y no una distribución normal), el ancho determinado por el nivel de confianza, si el pico no está centrado cerca del valor real, no contendrá la solución verdadera , aun estando esta muy cerca de la estimada. Es así que también, en algunos casos, la banda de nivel de confianza contiene a la solución verdadera, como en el caso del parámetro ”t2”. Esto es debido a que el pico de la distribución está más o menos centrado en el valor correcto. Esto nos da otra indicación más de que la relación de compresión estimada en los escalonamientos es mejor , normalmente, que las otras dos variables. Esto podría deberse a que el cálculo de los valores del ciclo está estrechamente y directamente relacionado con los ratios de presiones y que, las otras variables ( mc,η ) se calculan a través de más relaciones y aparecen de manera más indirecta en las ecuaciones. Esto último dificultaría el proceso del programa para optimizar los valores de las variables. Resultados t1 t2 t3 t4 t5 t6 t7 t8 t9 t10 t11 t12 t13 t14 t15 Nivel de confianza 0 1 0 0 1 0 0 1 0 0 1 1 0 0 0 A continuación,vamos a hacer mención a algunas gráficas que se adjuntan en el apéndice, las cuales son fáciles de consultar ya que están clasificadas por ternas y títulos, en orden. Comenzamos analizando algunos gráficos de la terna (1,1,1), la cual es la de menos degradación, para ilustrar sus buenos resultados. Para las tres primeras parejas de gráficas, las correspondientes a las soluciones para mc1,pr1 e η1 , observamos muy buenos resultados. En todas, como se observa en la tabla que se presenta a continuación con las medias, los valores encontrados mayoritariamente se encuentran en el entorno de 0.99, el cual es el valor solución. Se observa que, para las demás coronas de compresor, también se detecta una degradación, lo cual es incorrecto, pero esta en casi todos los casos es muy pequeña. Podríamos decir que cuando la solución de valores por encima de 0.995 y por debajo de 1.005 consideraremos que el programa ha supuesto degradación nula, por tanto estos resultados son ejemplo de una terna satisfactoriamente resuelta por el programa. Resultados t1 t2 t3 t4 t5 t6 t7 t8 t9 t10 t11 t12 t13 t14 t15 Media 0.9944 0.9906 0.9967 0.9945 0.9960 0.9975 0.9928 1.0008 0.9981 0.9822 1.0014 0.9981 0.9936 1.0017 0.9982 También presentamos los valores de los errores relativos en una tabla. En esta se evidencia que los errores para esta terna son mucho menores que los presentados para la anterior (3, 3, 5), indicación clara del buen funcionamiento del optimizador para cuando el nivel de degradación supuesto es igual en las tres variables consideradas. De hecho el mayor error relativo es de -0.6752 por ciento, valor completamente válido si suponemos un error notable a partir del 3% . Esta tabla y la anterior también muestran, de nuevo, la capacidad del programa de calcular con mayor precisión la degradación de la presión, ya que,puede observarse un error de cálculo del −0.0598% que es extremadamente bajo.También, en las gráficas puede verse que la clara moda de datos está en este valor de 0.99. 84 Capítulo 6. Estudio explorativo y estadístico de la capacidad de resolución de la herramienta de diagnóstico Resultados t1 t2 t3 t4 t5 t6 t7 t8 t9 t10 t11 t12 t13 t14 t15 Error relativo -0.4466 -0.0598 -0.6752 0.5548 0.4031 0.2503 0.7227 -0.0844 0.1862 1.7810 -0.1403 0.1912 0.6372 -0.1717 0.1805 Es importante mencionar también los buenos resultados que se extraen de los escalonamientos aguas abajo como el 3 y el 4. Si observamos la gráfica correspondiente a la degradación del gasto corregido en el escalonamiento 3, se puede ver un pico pronunciado en las barras amarillas para un valor cercano a 0.995, lo cual evidencia que el programa detecta la no degradación de este tercer escalonamiento en muchas ocasiones. También para el gasto corregido en el escalonamiento 4, podemos ver una media y una moda que coinciden en el valor 1, solución exacta del problema. Así, mostramos un ejemplo del optimizador siendo eficaz en mostrar un negativo verdadero, es decir, predecir con precisión que, en una corona en específico, no hay fallo. A continuación, se analizan algunos resultados obtenidos de la terna de fallo (5,5,5), la cual evidencia la mayor degradación estudiada en el trabajo estadístico. En este caso, podemos observar que,los resultados en el primer escalonamiento no son tan precisos como en la terna anterior. Se presentan en las gráficas correspondientes al gasto corregido y a la presión distribuciones binomiales, es decir, con dos picos en diferentes rangos de valores, lo cual nos da resultados poco significativos en la práctica. Es cierto que, aun no presentando una tendencia clara, de nuevo, en la predicción de la presión, la media que se calcula con los resultados es de 0.9586, lo cual predice bien la realidad. Sin embargo, los resultados medios para t1 y t3 son 0.9728 y 0.9812 respectivamente, los cuales son un poco errados. De hecho, para el rendimiento, el programa presenta un fuerte pico en su valor medio, lo cual quiere decir que el programa no ha sido capaz, casi en ninguno de los casos, de encontrar una solución cercana a 0.95. Estos resultados muestran un ejemplo de terna que, aun presentando mucha degradación (lo cual debería ser beneficioso para el resultado), nos da una solución no muy buena. Esto es parte del análisis de cualquier estudio estadístico y de la posibilidad de encontrar ternas que sean resueltas más fácilmente que otras. Interesantes son también los resultados presentados en ternas como la (5,3,3), que muestran comportamientos binomiales como los anteriormente mencionados. En la gráfica correspondiente al gasto corregido del primer escalonamiento de la terna (5,3,3), se observa que hay una barra más alta que las demás pero que, a su derecha también encontramos otro pico que no es acorde con una distribución normal Esto significa que la herramienta de diagnóstico, en muchos casos, ha encontrado la solución en 0.975, pero en otros muchos, aunque en menor medida, ha sido encontrada en 0.979. Ninguno de los resultados encontrados coincide con el verdadero 0.95 que deseamos pero es interesante ver este tipo de comportamientos para así comprender que el método es sensible al ruido aleatorio que le es introducido, dando así resultados de naturaleza distinta. El resultado de la relación de presión en la primera corona vuelve a ser mejor que los demás, presentando un valor medio entorno a 0.966, siendo el correcto 0.97. Mientras que el rendimiento se estima en 0.985 cuando es 0.97. De acuerdo con lo esperado de un estudio estadístico, también podemos encontrar entre la gran cantidad de resultados, algunos que presentan comportamientos extraños, fuera de tendencia y completamente errados. Este es el caso de la detección del parámetro Pr2 en la terna (5,3,3). Si observamos el histograma correspondiente nos daremos cuenta de que, además de no presentar una tendencia clara, la moda se presenta en un valor medio de 0.905, siendo la solución deseada 1. Este es un caso en donde, la relación de presiones, que suele estar tener resultados buenos, es completamente errada. Este resultado, como comentado al principio, es solo la evidencia de que en un estudio estadístico siempre tendremos casos fuera de la norma general y que pondrán en cuestión la validez de una regla o método. Sin embargo, este tipo de resultados no invalida los resultados satisfactorios observados con anterioridad. 6.3 Resultados 85 (a) a. (b) b. Figura 6.10 PDF e histograma. 6.3.1 Análisis de casos sueltos En este apartado analizaremos casos sueltos de degradación, esta vez teniendo en cuenta el efecto de la turbina y sus tres escalonamientos. Introduciremos degradación de las mismas magnitudes ( 1%,2%,5% ) trabajadas anteriormente, pero haciendo fallar la turbina también. De esta manera, el programa trabaja con un vector de incógnitas ”t” de longitud 24, ya que contamos con 8 turbomáquinas y 3 variables que se consideran objeto de fallo en cada una. Este análisis ya no constituye un estudio estadístico completo, pero nos permitirá obtener información sobre la validez de aplicación del método en el caso en el que se considera el motor completo y nos proporcionará con una idea de los errores que se cometen al añadir más incógnitas. En las tablas que se presentan posteriormente se han analizado 4 casos sueltos. El primero corresponde a la degradación de el primer escalonamiento de axial con la terna (1,1,1) y la degradación de la primera turbina de alta con la terna (1,1,1), como puede verse en la columna llamada ”solución”, donde aparece la solución real de la degradación. En el segundo caso se ha supuesto una degradación del primer escalonamiento con la terna de fallo (5,3,3) y otra de (1,1,1) en la primera turbina de alta. A continuación, en el tercer caso hemos considerado solo degradación en los dos escalonamientos de la turbina de baja, para poder ver los resultados que da el optimizador solo considerando errores al final del ciclo. Esta degradación constituye dos ternas iguales de valor (5,5,3). Por último, hemos analizado un caso en donde toda la turbina, tanto escalonamiento de alta como ambos de baja, se degradan, con ternas iguales de (5,5,3). En las tablas se puede observar que, los errores relativos son similares a los vistos en el caso del primer escalonamiento, aunque pueden llegar a valores de alrededor del −5% , como en el caso 3. Esto puede ser debido al aumento de incógnitas en el problema, que suma dificultad al trabajo del optimizador. Puede verse que, en general, el caso que mejor ha sido resuelto es el caso 1, para el cual tenemos menos elementos degradados, los cuales están en la zona anterior del motor. Este buen resultado puede deberse a que al optimizador le es más sencillo, desde la condición inicial que se le da (todo unos), dar resultados cercanos a 1, ya que busca un mínimo local, lo que hace que una degradación vecina a la inexistente, de resultados mejores. Sí es cierto también que con mayor degradación el sistema detecta mejor dónde se focaliza el error, ya que el fallo es más distintivo cuando el programa intenta comparar el caso real con ruido de su problema incógnita, que comienza, en la primera iteración, con degradación 0. También puede verse que en la columna de ”resultados” aparecen valores variopintos, dependiendo del caso que observemos. En el cuarto caso, por ejemplo, observamos que, en la zona primera del motor, el optimizador da muchos negativos (2) lo cual es correcto, además de dar, en la mayoría de casos donde hay degradación, un positivo (1). Esto significa que este tipo de terna con gran valor de degradación es bien detectada, incluso si se produce en una turbina y no en el compresor. 86 Capítulo 6. Estudio explorativo y estadístico de la capacidad de resolución de la herramienta de diagnóstico Sin embargo, en casos como el 3, la no degradación se detecta bien mientras que la degradación da bastantes falsos negativos, los cuales no supondrían un beneficio en la práctica. Además, cerca de donde se produce la degradación aparecen falsos positivos (3), los cuales pueden indicar que la herramienta de diagnóstico predice una degradación en un punto anterior a donde se va a producir, teniendo, a veces, un problema de ”desfase”. Es también interesante ver que, cuando las degradaciones están repartidas entre el compresor y la turbina, habiendo dos zonas degradadas bien diferenciadas, el optimizador hace un peor trabajo de detección del fallo. Esto puede verse en los casos 1 y 2, donde en la zona media no aparecen degradaciones. Aquí, el optimizador nos devuelve no solo negativos (2), sino falsos positivos (3), los cuales no deberían aparecer. Esto se debe a que, la no linealidad y acoplamiento del sistema, hacen difícil para la herramienta detectar dos zonas de degradación claras, dejando la zona entre ellas continua y con valor igual a 1. El optimizador parece querer seguir detectando fallos del mismo orden y por tanto, traslada la tendencia de degradación a la zona central. En general, se puede decir que, el optimizador hace una labor parecida al caso de la degradación de una corona de compresor pero presentando errores un poco mayores y problemas de detección cuando se degradan zonas diversas en un mismo instante, de manera que tiene problemas en individuarlas e aislarlas. Era de esperar que un ligero empeoramiento del comportamiento de la herramienta se diera al añadir más incógnitas. Resultados Caso 1 Error relativo Valor solución (X) Solución real Resultados t1 1.5290 0.9749 0.9900 1 t2 -0.5942 0.9959 0.9900 4 t3 -0.4385 0.9943 0.9900 1 t4 -1.4610 0.9854 1.0000 3 t5 -0.4182 1.0042 1.0000 2 t6 0.5683 0.9943 1.0000 3 t7 0.4359 0.9956 1.0000 2 t8 -0.8017 1.0080 1.0000 3 t9 0.5089 0.9949 1.0000 3 t10 0.0494 0.9995 1.0000 2 t11 -0.9384 1.0094 1.0000 3 t12 0.4892 0.9951 1.0000 2 t13 1.4400 0.9856 1.0000 3 t14 0.2432 0.9976 1.0000 2 t15 -0.5951 1.0060 1.0000 3 t16 1.1163 0.9789 0.9900 1 t17 -1.6314 1.0062 0.9900 3 t18 -0.1918 0.9919 0.9900 1 t19 2.0596 0.9794 1.0000 3 t20 -0.4345 1.0043 1.0000 2 t21 -1.3654 1.0137 1.0000 3 t22 0.1599 0.9984 1.0000 2 t23 -0.3027 1.0030 1.0000 2 t24 -1.2168 1.0122 1.0000 3 6.3 Resultados 87 Resultados Caso 2 Error relativo Valor solución (X) Solución real Resultados t1 -1.6445 0.9656 0.9500 1 t2 -0.8727 0.9785 0.9700 1 t3 -3.1603 1.0007 0.9700 4 t4 0.8828 0.9813 0.9900 1 t5 -0.1402 0.9914 0.9900 1 t6 -1.3120 1.0030 0.9900 4 t7 0.6547 0.9935 1.0000 3 t8 0.5703 0.9943 1.0000 3 t9 -0.3850 1.0038 1.0000 2 t10 0.3086 0.9969 1.0000 2 t11 0.4096 0.9959 1.0000 2 t12 -0.4889 1.0049 1.0000 2 t13 -1.4200 1.0142 1.0000 3 t14 -0.6765 1.0068 1.0000 3 t15 0.3250 0.9967 1.0000 2 t16 -1.3161 1.0030 0.9900 4 t17 -1.2282 1.0022 0.9900 4 t18 -0.4141 0.9941 0.9900 1 t19 -1.2214 1.0122 1.0000 3 t20 0.6106 0.9939 1.0000 3 t21 1.4954 0.9850 1.0000 3 t22 -0.1609 1.0016 1.0000 2 t23 0.5097 0.9949 1.0000 3 t24 1.2762 0.9872 1.0000 3 Resultados Caso 3 Error relativo Valor solución (X) Solución real Resultados t1 0.9203 0.9908 1.0000 3 t2 0.1819 0.9982 1.0000 2 t3 0.2019 0.9980 1.0000 2 t4 0.4209 0.9958 1.0000 2 t5 -0.3486 1.0035 1.0000 2 t6 0.2178 0.9978 1.0000 2 t7 0.1451 0.9985 1.0000 2 t8 -0.5669 1.0057 1.0000 3 t9 0.1955 0.9980 1.0000 2 t10 0.0319 0.9997 1.0000 2 t11 -0.6720 1.0067 1.0000 3 t12 0.1827 0.9982 1.0000 2 t13 0.7042 0.9930 1.0000 3 t14 -0.2145 1.0021 1.0000 2 t15 -0.5711 1.0057 1.0000 3 t16 0.9782 0.9902 1.0000 3 t17 -2.5757 1.0258 1.0000 3 t18 2.5719 0.9743 1.0000 3 t19 -3.2105 0.9805 0.9500 1 t20 -5.1232 0.9987 0.9500 4 t21 -1.8820 0.9883 0.9700 1 t22 -5.2042 0.9994 0.9500 4 t23 -5.0833 0.9983 0.9500 4 t24 -2.4987 0.9942 0.9700 1 94 Capítulo A. Anexo Tabla A.3 Compresor axial: rendimiento isentrópico. HHHH H β Nc0,499 0,559 0,635 0,738 0,804 0,862 0,928 1 1,045 1,077 1,109 1,14 00,690 0,695 0,701 0,695 0,680 0,665 0,645 0,622 0,610 0,605 0,600 0,598 0,07692 0,692 0,700 0,709 0,710 0,700 0,680 0,660 0,638 0,624 0,617 0,611 0,610 0,15385 0,692 0,702 0,714 0,724 0,720 0,700 0,676 0,653 0,639 0,631 0,623 0,620 0,23077 0,692 0,703 0,718 0,732 0,734 0,723 0,694 0,669 0,653 0,645 0,637 0,634 0,30769 0,692 0,704 0,720 0,737 0,742 0,740 0,720 0,687 0,667 0,659 0,650 0,648 0,38462 0,691 0,704 0,721 0,739 0,748 0,750 0,741 0,706 0,682 0,672 0,662 0,659 0,46154 0,691 0,703 0,721 0,741 0,752 0,760 0,754 0,730 0,695 0,683 0,671 0,668 0,53846 0,690 0,702 0,721 0,742 0,755 0,764 0,767 0,755 0,718 0,698 0,682 0,679 0,61538 0,689 0,701 0,719 0,742 0,757 0,769 0,779 0,770 0,739 0,711 0,690 0,685 0,69231 0,687 0,699 0,717 0,740 0,758 0,773 0,790 0,787 0,749 0,724 0,700 0,691 0,76923 0,686 0,697 0,715 0,739 0,758 0,775 0,792 0,792 0,760 0,731 0,704 0,695 0,84615 0,684 0,695 0,712 0,736 0,757 0,775 0,794 0,794 0,763 0,734 0,709 0,699 0,92308 0,682 0,693 0,710 0,734 0,755 0,775 0,794 0,795 0,765 0,736 0,710 0,700 10,680 0,690 0,706 0,730 0,753 0,773 0,793 0,795 0,768 0,740 0,711 0,696 Tabla A.4 Compresor axial: surge. mc1,02 1,33 1,81 2,54 3,225 3,85 4,525 5,14 5,605 5,75 5,89 5,94 PR 1,592 1,989 2,605 3,542 4,421 5,223 6,089 6,878 7,474 7,660 7,840 7,904 A.1.1 Turbina de gas Tabla A.5 Turbina de alta: gasto corregido. HHHH H β Nc0,4 0,5 0,6 0,7 0,8 0,9 1 1,1 1,2 011,7900 11,7700 11,7500 11,7300 11,7100 11,6900 11,6900 11,7100 11,7300 0,125 17,5219 17,5036 17,3000 17,0942 16,6992 16,3017 15,9841 15,6505 15,7708 0,25 19,5469 19,4584 19,3600 19,2119 19,0741 18,8091 18,5819 18,3541 18,0475 0,375 19,9941 19,9470 19,8900 19,7916 19,7117 19,5792 19,4266 19,2567 19,0413 0,5 20,1113 20,0919 20,0694 20,0275 19,9919 19,8888 19,7969 19,6850 19,5400 0,625 20,1248 20,1200 20,1152 20,1055 20,0855 20,0267 19,9670 19,8770 19,7673 0,75 20,1131 20,1131 20,1200 20,1200 20,1000 20,0938 20,0506 19,9675 19,8744 0,875 20,0966 20,0983 20,1200 20,1200 20,1000 20,1000 20,0700 20,0017 19,9234 120,0800 20,0900 20,1200 20,1200 20,1000 20,1000 20,0700 20,0100 19,9400 Tabla A.6 Turbina de alta: rendimiento isentrópico. HHHH H β Nc0,4 0,5 0,6 0,7 0,8 0,9 1 1,1 1,2 00,55 0,56 0,56 0,56 0,56 0,55 0,54 0,47 0,39 0,125 0,78391 0,84375 0,88977 0,90297 0,88414 0,82519 0,77945 0,71656 0,60672 0,25 0,875 0,77875 0,82844 0,89375 0,91906 0,92397 0,89622 0,86203 0,79344 0,375 0,75 0,74563 0,78703 0,86063 0,89422 0,92314 0,92852 0,91586 0,88453 0,5 0,70625 0,72625 0,75938 0,8275 0,87075 0,91063 0,93194 0,93463 0,922 0,625 0,55 0,70711 0,74227 0,80211 0,85169 0,89984 0,92584 0,93848 0,93664 0,75 0,78391 0,69156 0,73156 0,78313 0,83581 0,88469 0,91688 0,93594 0,94306 0,875 0,75875 0,68414 0,71914 0,76914 0,81911 0,86656 0,90242 0,92828 0,93825 10,72563 0,68 0,715 0,765 0,81 0,85 0,89 0,88453 0,925 A.2 Anexo B: Gráficas del estudio estadístico 95 Tabla A.7 Turbina de alta: presiones mínimas y máximas. Nc0,4 0,5 0,6 0,7 0,8 0,9 1 1,1 1,2 Presión mínima 1,15 1,15 1,15 1,15 1,15 1,15 1,15 1,15 1,15 Presión máxima 3,8 3,8 3,8 3,8 3,8 3,8 3,8 3,8 3,8 A.1.2 Turbina de potencia Tabla A.8 Turbina de baja: gasto corregido. HHHH H β Nc0,4 0,5 0,6 0,7 0,8 0,9 1 1,1 1,2 03 2,985 2,97 2,96 2,95 2,94 2,92 2,91 2,9 0,08333 3,855 3,84 3,8 3,7625 3,7125 3,64 3,5825 3,54 3,505 0,16667 4,455 4,455 4,395 4,32 4,23 4,14 4,065 4,015 3,97 0,25 4,71 4,71 4,66 4,6 4,545 4,4525 4,39 4,3175 4,2775 0,33333 4,9 4,89 4,85 4,8 4,76 4,71 4,65 4,56 4,54 0,41667 4,995 4,97 4,95 4,905 4,885 4,8575 4,8175 4,7525 4,705 0,5 5,015 5,005 4,99 4,965 4,95 4,925 4,89 4,845 4,8 0,58333 5,0275 5,0175 5,0075 4,995 4,985 4,97 4,94 4,9075 4,8575 0,66667 5,035 5,025 5,015 5,01 5 4,99 4,975 4,95 4,91 0,75 5,04 5,03125 5,02125 5,02125 5,01125 5,0025 5,00125 4,9825 4,95375 0,83333 5,04 5,0375 5,0275 5,0275 5,0175 5,015 5,0075 4,995 4,9725 0,91667 5,04 5,04 5,03 5,03 5,02 5,02 5,01 5 4,98 15,04 5,04 5,03 5,03 5,02 5,02 5,01 5 4,98 Tabla A.9 Turbina de baja: rendimiento isentrópico. HHHH H β Nc0,4 0,5 0,6 0,7 0,8 0,9 1 1,1 1,2 00,745 0,707 0,635 0,56 0,52 0,47 0,4 0,3 0,2 0,08333 0,769 0,77875 0,7535 0,70625 0,67875 0,64125 0,59475 0,5 0,36 0,16667 0,7755 0,805 0,809 0,7945 0,782 0,7575 0,726 0,6545 0,5475 0,25 0,77225 0,8155 0,833 0,83925 0,8425 0,83175 0,81075 0,758 0,68625 0,33333 0,762 0,812 0,842 0,859 0,87 0,87 0,86 0,83 0,802 0,41667 0,746 0,80075 0,84 0,8595 0,877 0,88175 0,87825 0,86225 0,82625 0,5 0,7255 0,783 0,8325 0,856 0,876 0,8805 0,882 0,8755 0,8535 0,58333 0,708 0,76375 0,81675 0,8465 0,86975 0,877 0,882 0,88125 0,871 0,66667 0,695 0,746 0,795 0,828 0,8655 0,8695 0,879 0,8845 0,8805 0,75 0,68362 0,72938 0,77362 0,8085 0,86037 0,86037 0,8745 0,88675 0,88737 0,83333 0,67675 0,71625 0,75675 0,791 0,84725 0,84725 0,867 0,8855 0,88925 0,91667 0,67175 0,70688 0,74213 0,775 0,82438 0,83563 0,85838 0,88237 0,89 10,668 0,7 0,729 0,76 0,795 0,825 0,849 0,878 0,89 Tabla A.10 Turbina de baja: presiones mínimas y máximas. HHHH H β Nc0,4 0,5 0,6 0,7 0,8 0,9 1 1,1 1,2 Presión mínima 1,15 1,15 1,15 1,15 1,15 1,15 1,15 1,15 1,15 Presión máxima 1,65 1,65 1,65 1,65 1,65 1,65 1,65 1,65 1,65 A.2 Anexo B: Gráficas del estudio estadístico 96 Capítulo A. Anexo A.3 Terna 111 (a) a. (b) b. Figura A.1 PDF e histograma. (a) a. (b) b. Figura A.2 PDF e histograma. A.3 Terna 111 97 (a) a. (b) b. Figura A.3 PDF e histograma. (a) a. (b) b. Figura A.4 PDF e histograma. (a) a. (b) b. Figura A.5 PDF e histograma. 98 Capítulo A. Anexo (a) a. (b) b. Figura A.6 PDF e histograma. (a) a. (b) b. Figura A.7 PDF e histograma. (a) a. (b) b. Figura A.8 PDF e histograma. A.3 Terna 111 99 (a) a. (b) b. Figura A.9 PDF e histograma. (a) a. (b) b. Figura A.10 PDF e histograma. (a) a. (b) b. Figura A.11 PDF e histograma. 100 Capítulo A. Anexo (a) a. (b) b. Figura A.12 PDF e histograma. (a) a. (b) b. Figura A.13 PDF e histograma. (a) a. (b) b. Figura A.14 PDF e histograma. A.3 Terna 111 101 (a) a. (b) b. Figura A.15 PDF e histograma. 102 Capítulo A. Anexo A.4 Terna 113 (a) a. (b) b. Figura A.16 PDF e histograma. (a) a. (b) b. Figura A.17 PDF e histograma. A.4 Terna 113 103 (a) a. (b) b. Figura A.18 PDF e histograma. (a) a. (b) b. Figura A.19 PDF e histograma. (a) a. (b) b. Figura A.20 PDF e histograma. 110 Capítulo A. Anexo (a) a. (b) b. Figura A.36 PDF e histograma. (a) a. (b) b. Figura A.37 PDF e histograma. (a) a. (b) b. Figura A.38 PDF e histograma. A.5 Terna 115 111 (a) a. (b) b. Figura A.39 PDF e histograma. (a) a. (b) b. Figura A.40 PDF e histograma. (a) a. (b) b. Figura A.41 PDF e histograma. 112 Capítulo A. Anexo (a) a. (b) b. Figura A.42 PDF e histograma. (a) a. (b) b. Figura A.43 PDF e histograma. (a) a. (b) b. Figura A.44 PDF e histograma. A.5 Terna 115 113 (a) a. (b) b. Figura A.45 PDF e histograma. 114 Capítulo A. Anexo A.6 Terna 131 (a) a. (b) b. Figura A.46 PDF e histograma. (a) a. (b) b. Figura A.47 PDF e histograma. A.6 Terna 131 115 (a) a. (b) b. Figura A.48 PDF e histograma. (a) a. (b) b. Figura A.49 PDF e histograma. (a) a. (b) b. Figura A.50 PDF e histograma. 116 Capítulo A. Anexo (a) a. (b) b. Figura A.51 PDF e histograma. (a) a. (b) b. Figura A.52 PDF e histograma. (a) a. (b) b. Figura A.53 PDF e histograma. A.6 Terna 131 117 (a) a. (b) b. Figura A.54 PDF e histograma. (a) a. (b) b. Figura A.55 PDF e histograma. (a) a. (b) b. Figura A.56 PDF e histograma. 118 Capítulo A. Anexo (a) a. (b) b. Figura A.57 PDF e histograma. (a) a. (b) b. Figura A.58 PDF e histograma. (a) a. (b) b. Figura A.59 PDF e histograma. A.6 Terna 131 119 (a) a. (b) b. Figura A.60 PDF e histograma. 126 Capítulo A. Anexo A.8 Terna 135 (a) a. (b) b. Figura A.76 PDF e histograma. (a) a. (b) b. Figura A.77 PDF e histograma. A.8 Terna 135 127 (a) a. (b) b. Figura A.78 PDF e histograma. (a) a. (b) b. Figura A.79 PDF e histograma. (a) a. (b) b. Figura A.80 PDF e histograma. 128 Capítulo A. Anexo (a) a. (b) b. Figura A.81 PDF e histograma. (a) a. (b) b. Figura A.82 PDF e histograma. (a) a. (b) b. Figura A.83 PDF e histograma. A.8 Terna 135 129 (a) a. (b) b. Figura A.84 PDF e histograma. (a) a. (b) b. Figura A.85 PDF e histograma. (a) a. (b) b. Figura A.86 PDF e histograma. 130 Capítulo A. Anexo (a) a. (b) b. Figura A.87 PDF e histograma. (a) a. (b) b. Figura A.88 PDF e histograma. (a) a. (b) b. Figura A.89 PDF e histograma. A.8 Terna 135 131 (a) a. (b) b. Figura A.90 PDF e histograma. 132 Capítulo A. Anexo A.9 Terna 151 (a) a. (b) b. Figura A.91 PDF e histograma. (a) a. (b) b. Figura A.92 PDF e histograma. A.9 Terna 151 133 (a) a. (b) b. Figura A.93 PDF e histograma. (a) a. (b) b. Figura A.94 PDF e histograma. (a) a. (b) b. Figura A.95 PDF e histograma. 134 Capítulo A. Anexo (a) a. (b) b. Figura A.96 PDF e histograma. (a) a. (b) b. Figura A.97 PDF e histograma. (a) a. (b) b. Figura A.98 PDF e histograma. A.9 Terna 151 135 (a) a. (b) b. Figura A.99 PDF e histograma. (a) a. (b) b. Figura A.100 PDF e histograma. (a) a. (b) b. Figura A.101 PDF e histograma. 238 Capítulo A. Anexo (a) a. (b) b. Figura A.357 PDF e histograma. (a) a. (b) b. Figura A.358 PDF e histograma. (a) a. (b) b. Figura A.359 PDF e histograma. A.26 Terna 535 239 (a) a. (b) b. Figura A.360 PDF e histograma. 240 Capítulo A. Anexo A.27 Terna 551 (a) a. (b) b. Figura A.361 PDF e histograma. (a) a. (b) b. Figura A.362 PDF e histograma. A.27 Terna 551 241 (a) a. (b) b. Figura A.363 PDF e histograma. (a) a. (b) b. Figura A.364 PDF e histograma. (a) a. (b) b. Figura A.365 PDF e histograma. 242 Capítulo A. Anexo (a) a. (b) b. Figura A.366 PDF e histograma. (a) a. (b) b. Figura A.367 PDF e histograma. (a) a. (b) b. Figura A.368 PDF e histograma. A.27 Terna 551 243 (a) a. (b) b. Figura A.369 PDF e histograma. (a) a. (b) b. Figura A.370 PDF e histograma. (a) a. (b) b. Figura A.371 PDF e histograma. 244 Capítulo A. Anexo (a) a. (b) b. Figura A.372 PDF e histograma. (a) a. (b) b. Figura A.373 PDF e histograma. (a) a. (b) b. Figura A.374 PDF e histograma. A.27 Terna 551 245 (a) a. (b) b. Figura A.375 PDF e histograma. 246 Capítulo A. Anexo A.28 Terna 553 (a) a. (b) b. Figura A.376 PDF e histograma. (a) a. (b) b. Figura A.377 PDF e histograma. A.28 Terna 553 247 (a) a. (b) b. Figura A.378 PDF e histograma. (a) a. (b) b. Figura A.379 PDF e histograma. (a) a. (b) b. Figura A.380 PDF e histograma. 254 Capítulo A. Anexo (a) a. (b) b. Figura A.396 PDF e histograma. (a) a. (b) b. Figura A.397 PDF e histograma. (a) a. (b) b. Figura A.398 PDF e histograma. A.29 Terna 555 255 (a) a. (b) b. Figura A.399 PDF e histograma. (a) a. (b) b. Figura A.400 PDF e histograma. (a) a. (b) b. Figura A.401 PDF e histograma. 256 Capítulo A. Anexo (a) a. (b) b. Figura A.402 PDF e histograma. (a) a. (b) b. Figura A.403 PDF e histograma. (a) a. (b) b. Figura A.404 PDF e histograma. A.29 Terna 555 257 (a) a. (b) b. Figura A.405 PDF e histograma. Índice de Figuras 1.1 Diagrama conceptual del trabajo 2 1.2 Prueba en banco de motor aeronáutico 2 2.1 Esquema interno del motor PT6A-65B 5 2.2 Esquema interno del motor PT6A-65B con nombres 6 2.3 Posición de sensores en el esquema de un motor genérico de la serie PT6-A según el manual oficial de entrenamiento de Pratt and Whitney. 7 2.4 Regulación y filtración del aceite 8 2.5 Sistema de combustible 8 2.6 Sistema de ignición 9 2.7 Esquema parte del motor enumeradas 9 2.8 Sistema pre-swirl 10 2.9 Relaciones de reducción 10 2.10 Parámetros usados en el ciclo programado 11 2.11 Temperaturas y presiones de un estudio de simulación computacional para el PT6A-65 12 2.12 Plano del Beechcraft 1900C-1 13 2.13 Esquema de ala con estrechamiento y sus partes 14 2.14 Cuerda media aerodinámica total para un estrechamiento doble 14 2.15 Cálculo de cuerda media aerodinámica siendo cr y ct valores de cuerda en la raiz y la punta del tramo considerado 14 2.16 Ecuación de la polar parabólica 15 2.17 Coeficiente CD0donde en la expresión se multiplica por 1.1 para contar con contribuciones misceláneas 15 2.18 Factores debidos a contribuciones varias a la resistencia aerodinámica 16 2.19 Coeficiente de fricción para flujo turbulento en placa plana 16 2.20 Número de Reynolds para el flujo para cada componente del avión k, siendo lk la longitud de referencia de cada componente k del avión 16 2.21 Factores de interferencia de cada componente 16 2.22 Calculo de los factores de forma FFk . Estos factores incluyen los efectos de grosor donde el término t c es el "thickness ratio" del perfil (grosor máximo dividido entre la cuerda) y l d es el "fineness ratio" del perfil (longitud dividida entre diámetro máximo) 17 2.23 Esquema del avión simplificado vista en alzado 17 2.24 Esquema del avión simplificado vista en planta. Generatriz del cono=2.315my radio del cono= 0.827m17 2.25 Fuerzas aerodinámicas en vuelo en crucero 19 2.26 Entalpía 20 2.27 Φ, Entropía a presión constante 20 2.28 Relación de compresión 20 2.29 Relación de Gibbs 20 2.30 Diferencia de Φsolo si el proceso es isoentrópico 20 2.31 Relación de compresión 20 2.32 Relación de expansión 20 2.33 Coeficientes para propiedades termodinámicas en condiciones estándar 21 259 260 Índice de Figuras 2.34 Ecuación de la cámara de combustión adimensionalizada con el gasto de aire a la entrada del motor 24 2.35 Diagrama h-s zona de la turbina 25 2.36 Diagrama h-s completo 25 2.37 Figura de motor de referencia esquematizada 28 2.38 Gráfico empuje específico con relación de compresión 32 2.39 Gráfico TSFC con relación e compresión 32 2.40 Gráfico empuje específico con parámetro theta 33 2.41 Gráfico TSFC con parámetro theta 33 2.42 Límites de inflamabilidad de las mezclas fuel-aire 34 2.43 Superficie 3D empuje específico 35 2.44 Superficie 3D TSFC 35 2.45 Gráfico temperatura de la cámara de combustión con el dosado 36 2.46 Curvas isoπe isoθ37 2.47 Curvas isoπe isoθcualitativas genéricas 37 2.48 Curvas isoπe isoθ38 2.49 Curvas isoπe isoθrespecto al consumo específico 38 2.50 Mapa de la hélice de referencia 39 2.51 Mapa de la hélice digitalizado 40 3.1 Mapa generalizado de un compresor axial 45 3.2 Mapa generalizado de un compresor centrífugo 46 3.3 Mapa generalizado de una turbina 46 3.4 Mapa escalado del primer escalonamiento de axial 48 3.5 Mapa escalado del segundo escalonamiento de axial 49 3.6 Mapa escalado del tercer escalonamiento de axial 49 3.7 Mapa escalado del cuarto escalonamiento de axial 50 3.8 Mapa escalado del compresor centrífugo: gasto corregido 51 3.9 Mapa escalado del compresor centrífugo: rendimiento 52 3.10 Mapa escalado de la turbina de alta 53 3.11 Mapa escalado del primer escalonamiento de turbina de baja 53 3.12 Mapa escalado del segundo escalonamiento de turbina de baja 54 3.13 Diseño del motor para el ciclo off design 55 3.14 Mapa del primer escalonamiento de axial con 4 condiciones de vuelo 61 3.15 Mapa del segundo escalonamiento de axial con 4 condiciones de vuelo 61 3.16 Mapa del tercer escalonamiento de axial con 4 condiciones de vuelo 62 3.17 Mapa del cuarto escalonamiento de axial con 4 condiciones de vuelo 62 3.18 Mapa del escalonamiento de compresor centrífugo con 4 condiciones de vuelo 63 3.19 Mapa del escalonamiento de turbina de alta con 4 condiciones de vuelo 63 3.20 Mapa del primer escalonamiento de turbina de baja con 4 condiciones de vuelo 64 3.21 Mapa del segundo escalonamiento de turbina de baja con 4 condiciones de vuelo 64 4.1 Metodología de la herramienta de optimización 70 5.1 Cálculo de un mínimo con algoritmo de optimización 72 5.2 Problema del cálculo de mínimos con optimización basada en gradientes conjugados 72 6.1 Prueba de corte de la PDF multivariable 76 6.2 Histograma con diferentes anchos para t1 y t2 de terna 335 79 6.3 Histograma con diferentes anchos para t3 y t4 de terna 335 80 6.4 PDF e histograma 81 6.5 PDF e histograma 81 6.6 PDF e histograma 82 6.7 PDF e histograma 82 6.8 PDF e histograma 83 6.9 PDF e histograma 83 6.10 PDF e histograma 87 Índice de Figuras 261 A.1 PDF e histograma 98 A.2 PDF e histograma 98 A.3 PDF e histograma 99 A.4 PDF e histograma 99 A.5 PDF e histograma 99 A.6 PDF e histograma 100 A.7 PDF e histograma 100 A.8 PDF e histograma 100 A.9 PDF e histograma 101 A.10 PDF e histograma 101 A.11 PDF e histograma 101 A.12 PDF e histograma 102 A.13 PDF e histograma 102 A.14 PDF e histograma 102 A.15 PDF e histograma 103 A.16 PDF e histograma 104 A.17 PDF e histograma 104 A.18 PDF e histograma 105 A.19 PDF e histograma 105 A.20 PDF e histograma 105 A.21 PDF e histograma 106 A.22 PDF e histograma 106 A.23 PDF e histograma 106 A.24 PDF e histograma 107 A.25 PDF e histograma 107 A.26 PDF e histograma 107 A.27 PDF e histograma 108 A.28 PDF e histograma 108 A.29 PDF e histograma 108 A.30 PDF e histograma 109 A.31 PDF e histograma 110 A.32 PDF e histograma 110 A.33 PDF e histograma 111 A.34 PDF e histograma 111 A.35 PDF e histograma 111 A.36 PDF e histograma 112 A.37 PDF e histograma 112 A.38 PDF e histograma 112 A.39 PDF e histograma 113 A.40 PDF e histograma 113 A.41 PDF e histograma 113 A.42 PDF e histograma 114 A.43 PDF e histograma 114 A.44 PDF e histograma 114 A.45 PDF e histograma 115 A.46 PDF e histograma 116 A.47 PDF e histograma 116 A.48 PDF e histograma 117 A.49 PDF e histograma 117 A.50 PDF e histograma 117 A.51 PDF e histograma 118 A.52 PDF e histograma 118 A.53 PDF e histograma 118 A.54 PDF e histograma 119 A.55 PDF e histograma 119 A.56 PDF e histograma 119 A.57 PDF e histograma 120 262 Índice de Figuras A.58 PDF e histograma 120 A.59 PDF e histograma 120 A.60 PDF e histograma 121 A.61 PDF e histograma 122 A.62 PDF e histograma 122 A.63 PDF e histograma 123 A.64 PDF e histograma 123 A.65 PDF e histograma 123 A.66 PDF e histograma 124 A.67 PDF e histograma 124 A.68 PDF e histograma 124 A.69 PDF e histograma 125 A.70 PDF e histograma 125 A.71 PDF e histograma 125 A.72 PDF e histograma 126 A.73 PDF e histograma 126 A.74 PDF e histograma 126 A.75 PDF e histograma 127 A.76 PDF e histograma 128 A.77 PDF e histograma 128 A.78 PDF e histograma 129 A.79 PDF e histograma 129 A.80 PDF e histograma 129 A.81 PDF e histograma 130 A.82 PDF e histograma 130 A.83 PDF e histograma 130 A.84 PDF e histograma 131 A.85 PDF e histograma 131 A.86 PDF e histograma 131 A.87 PDF e histograma 132 A.88 PDF e histograma 132 A.89 PDF e histograma 132 A.90 PDF e histograma 133 A.91 PDF e histograma 134 A.92 PDF e histograma 134 A.93 PDF e histograma 135 A.94 PDF e histograma 135 A.95 PDF e histograma 135 A.96 PDF e histograma 136 A.97 PDF e histograma 136 A.98 PDF e histograma 136 A.99 PDF e histograma 137 A.100 PDF e histograma 137 A.101 PDF e histograma 137 A.102 PDF e histograma 138 A.103 PDF e histograma 138 A.104 PDF e histograma 138 A.105 PDF e histograma 139 A.106 PDF e histograma 140 A.107 PDF e histograma 140 A.108 PDF e histograma 141 A.109 PDF e histograma 141 A.110 PDF e histograma 141 A.111 PDF e histograma 142 A.112 PDF e histograma 142 A.113 PDF e histograma 142 A.114 PDF e histograma 143 Índice de Figuras 263 A.115 PDF e histograma 143 A.116 PDF e histograma 143 A.117 PDF e histograma 144 A.118 PDF e histograma 144 A.119 PDF e histograma 144 A.120 PDF e histograma 145 A.121 PDF e histograma 146 A.122 PDF e histograma 146 A.123 PDF e histograma 147 A.124 PDF e histograma 147 A.125 PDF e histograma 147 A.126 PDF e histograma 148 A.127 PDF e histograma 148 A.128 PDF e histograma 148 A.129 PDF e histograma 149 A.130 PDF e histograma 149 A.131 PDF e histograma 149 A.132 PDF e histograma 150 A.133 PDF e histograma 150 A.134 PDF e histograma 150 A.135 PDF e histograma 151 A.136 PDF e histograma 152 A.137 PDF e histograma 152 A.138 PDF e histograma 153 A.139 PDF e histograma 153 A.140 PDF e histograma 153 A.141 PDF e histograma 154 A.142 PDF e histograma 154 A.143 PDF e histograma 154 A.144 PDF e histograma 155 A.145 PDF e histograma 155 A.146 PDF e histograma 155 A.147 PDF e histograma 156 A.148 PDF e histograma 156 A.149 PDF e histograma 156 A.150 PDF e histograma 157 A.151 PDF e histograma 158 A.152 PDF e histograma 158 A.153 PDF e histograma 159 A.154 PDF e histograma 159 A.155 PDF e histograma 159 A.156 PDF e histograma 160 A.157 PDF e histograma 160 A.158 PDF e histograma 160 A.159 PDF e histograma 161 A.160 PDF e histograma 161 A.161 PDF e histograma 161 A.162 PDF e histograma 162 A.163 PDF e histograma 162 A.164 PDF e histograma 162 A.165 PDF e histograma 163 A.166 PDF e histograma 164 A.167 PDF e histograma 164 A.168 PDF e histograma 165 A.169 PDF e histograma 165 A.170 PDF e histograma 165 A.171 PDF e histograma 166 Bibliografía [1] KLM (2016):¿Cómo se prueba el motor de un avión?, URL:https:jetnews.com.mx/se-prueba-motor-avion/ [2] Avutarda (2016):El papel de la relación Lift y Drag en el diseño de aviones , URL:https://greatbustardsflight.blogspot.com/2016/02/el-papel-de-la-relacion-lift-drag-en-el.html [3] Anónimo (2021):Pratt and Whitney Canada PT6,Wikipedia. [4] Euravia (2021):PT6A-65B Capabilities, URL:https://www.euravia.aero/engine-detail/pt6a-65b?locale=en [5] Anónimo (2022):Artificial Neural Network,Wikipedia. [6] Zhang Yingao (2016-2017):WJ-25 Design Proposal [7] David G.Hull:Fundamentals of Airplane Flight Mechanics, Editorial Springers. [8] Pratt and Whitney Canada:PT6A-60 series training manual, Pratt and Whitney Canada, Inc. [9] United Turbine:PT6: Descriptive course and guide to troubleshooting. [10] European Aviation Safety Agency (2013):Easa Type-Certificate Data Sheet, TDS. [11] Pablo Hinojosa López (2018):Análisis de Turborreactores: Operación Fuera de Diseño, Universidad de Sevilla. [12] Gonzalo Zamora Caballero(2021):Modelo de Aeronave con propulsión híbrida, Universidad de Sevilla. [13] Department of Transportation Federal Aviation Administration(2008):Type-Certificate Data Sheet NO a24CE,John Wiley. [14] U.S Department of Transportation (2012):Type-Certificate Data Sheet NO E4EA 271