Comparación modelos turbojet

Universidad Politécnica de Valencia Aerorreactores y Motores Cohete COMPARACIÓN DE MODELOS DE CICLO TURBORREACTOR Jorge García Tíscar 16 de noviembre de 2011 Resumen En el presente estudio, se ha propuesto un motor real de referencia (Junkers Jumo 004B) a partir de cuyos datos se han aplicado una serie de modelos ter- modinámicos del ciclo turborreactor, realizándose una comparación entre los distintos resultados de empuje y consumo específico predichos por cada uno de ellos con las prestaciones del motor recuperadas de la literatura existente. Índice Nomenclatura 4 1. Introducción 5 1.1. Objetivo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 1.2. Metodología . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 2. Motor de referencia 6 3. Modelo de gas perfecto 8 3.1. Análisis por estaciones . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 3.2. Prestaciones y figuras de mérito . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 3.3. Resultados del cálculo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 4. Modelo de gas semiperfecto 12 4.1. Análisis por estaciones . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 4.2. Resultados del cálculo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 5. Modelo de equilibrio 16 5.1. Análisis por estaciones . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16 5.2. Resultados del cálculo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19 6. Software de cálculo 20 7. Variaciones de los modelos 23 7.1. Modelo de gas perfecto con 2 constantes . . . . . . . . . . . . . . . . . . . . . . . 23 7.2. Modelo de gas semiperfecto (correlación 2) . . . . . . . . . . . . . . . . . . . . . 23 8. Discusión 24 9. Conclusiones 26 Referencias 26 A. Ejemplo de salida CEA 27 B. Ejemplo de ciclo con Mathematica 30 2 Índice de figuras 2.1. Messerschmitt Me 262, tres vistas . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 2.2. Diagrama del Jumo 004B . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 2.3. Diagrama de flujos de aire . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 6.1. Ventana de datos principales de GasTurb . . . . . . . . . . . . . . . . . . . . . . . 20 6.2. Ventana de optimización de GasTurb . . . . . . . . . . . . . . . . . . . . . . . . . . 21 6.3. Nomenclatura y circuito de aire en GasTurb . . . . . . . . . . . . . . . . . . . . . 22 6.4. Diagrama del ciclo TS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22 8.1. Resultados relativos a los datos del motor real . . . . . . . . . . . . . . . . . . . 24 8.2. Comparación de ηt h respecto a GasTurb . . . . . . . . . . . . . . . . . . . . . . . 25 Índice de tablas 2.1. Prestaciones del motor original . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 3.1. Prestaciones con el modelo de propiedades constantes . . . . . . . . . . . . 11 4.1. Prestaciones con el modelo de propiedades variables . . . . . . . . . . . . . . 15 5.1. Resultado de la iteración en CC . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17 5.2. Prestaciones con el modelo de equilibrio . . . . . . . . . . . . . . . . . . . . . . . 19 6.1. Resultados del ciclo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 6.2. Prestaciones con el modelo de GasTurb . . . . . . . . . . . . . . . . . . . . . . . . 22 7.1. Prestaciones con el modelo de propiedades constantes (2 val.) . . . . . . . 23 7.2. Prestaciones con la 2ª correlación. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 8.1. Resultados absolutos de los distintos modelos . . . . . . . . . . . . . . . . . . . 24 8.2. Resultados relativos a los datos del motor real . . . . . . . . . . . . . . . . . . . 24 3 Nomenclatura La nomenclatura que se utilizará se describe a continuación. La numeración de las estaciones del motor corresponde a la habitual (SAE), el subíndice 0 indica magnitu- des totales y el subíndice s indica magnitudes ideales (isentrópicas). Símbolo Significado Unidades c1 Porcentaje m˙a de refrigeración estátor – c2 Porcentaje m˙a de refrigeración rotor – cp c Calor específico a presión constante flujo caliente J/(K kg) cp f Calor específico a presión constante flujo frío J/(K kg) E Empuje neto N Es p Empuje específico m/s f Dosado relativo – h Entalpía J H f Poder calorífico específico de combustible KJ/kg m˙a Flujo másico de aire kg/s m˙ f Flujo másico de combustible kg/s R Constante gas ideal específica aire J/(K kg) rd Coeficiente de recuperación de presión de parada – TSF C Consumo específico respecto al empuje g/(kN s) V0 Velocidad de vuelo m/s V9 Velocidad salida tobera m/s γc Relación calores específicos cp/cv flujo caliente – γ f Relación calores específicos cp/cv flujo frío – ∆pc c Pérdida presión en cámara de combustión Pa ηc Rendimiento compresor – ηc c Rendimiento cámara de combustión – ηd Rendimiento difusor – ηn Rendimiento tobera – ηm Eficiencia mecánica – ηt Rendimiento turbina – ηt h Eficiencia térmica – pic Relación de compresión – 4 1 Introducción La manera más común de evaluar teóricamente las prestaciones de los motores a reacción consiste en analizar su ciclo termodinámico. Sin embargo, existen muchas maneras de realizar dicho análisis, en función de la profundidad y de las hipótesis de trabajo que se realicen con el fin de simplificar el modelo. 1.1 Objetivo El objetivo del presente trabajo es comparar los distintos modelos que se emplean para analizar el ciclo termodinámico de los motores a reacción, desde los más sim- ples a los rigurosos; en este caso en un turborreactor de flujo único. Se pretende de esta manera discutir si el aumento de precisión en el análisis justifica el incremento de complejidad (y por tanto de coste computacional) correspondiente y de esta manera intentar encontrar si algún modelo representa un balance idóneo entre exactitud y simplicidad. 1.2 Metodología Para ello se va a seguir el siguiente método: en primer lugar se va a proponer un mo- tor de ejemplo, cuyos parámetros se conozcan o se puedan estimar. A continuación se va a analizar dicho motor mediante los siguientes modelos de cálculo: Modelo de gas perfecto: se supondrán constantes las características del fluido de trabajo, siendo éstas las correspondientes antes de la combustión. Modelo de gas semiperfecto: el siguiente paso en complejidad es suponer pro- piedades variables con la temperatura, empleando para ello regresiones poli- nómicas. Modelo de equilibro: por último, se emplea un programa de cálculo de compo- sición de equilibro químico para calcular las propiedades termodinámicas. Software de cálculo: a modo de comprobación de los resultados anteriores, se emplea un software de cálculo de ciclos empleando los mismos parámetros. Mediante cada uno de ellos se calcularán las prestaciones de dicho motor de ejem- plo suponiendo que se encuentra en un banco de pruebas (T0 = 288,15 K, p0 = 101,325 kPa, V0 = 0 m/s, combustible JP-4 con H f = 43,1 MJ/kg) y se compararán los distintos resultados. También se considerarás modificaciones o refinamientos de los métodos base con- siderados, de tal manera que se pueda comprobar si alguno de ellos se puede mejorar añadiendo alguna consideración para mejorar el cálculo. 5 2 Motor de referencia El motor que se ha elegido es el clásico Junkers Jumo 004B, el primer motor a reac- ción en ser producido en serie. Dicho motor se comenzó a desarrollar en Alemania el otoño de 1939, siendo su primera prueba en octubre de 1940 y su primer vuelo con el Messerschmitt Me 262 Schwalbe1 en julio de 1942. Figura 2.1: Messerschmitt Me 262, tres vistas A principios de 1944 se comenzó su producción en serie . En mayo de 1945 se ha- bían entregado 6000 motores a la Luftwaffe. El Jumo 004, a pesar de ser uno de los primeros turbojets en entrar en servicio, presenta características muy innovadoras, como un área de tobera variable, turbomaquinaria axial (frente a la maquinaria ra- dial inicial) e incluso, como se ha dicho, una turbina refrigerada por aire. Figura 2.2: Diagrama del Jumo 004B Dicha refrigeración es consecuencia de las limitaciones de la economía alemana durante la guerra, concretamente de la escasez de níquel para fabricar aleaciones re- sistentes al calor [1]. Las primeras versiones de desarrollo sí las empleaban, pero a la hora de producir en serie el motor se hizo evidente que las pocas reservas con las que contaba Alemania serían insuficientes. 1Golondrina. El apodo viene tanto de la figura como de la alta velocidad 6 COMPARACIÓN MODELOS 2 MOTOR DE REFERENCIA Figura 2.3: Diagrama de flujos de aire La solución a este problema se basó en sangrar aire del compresor (más de un 7 % del flujo másico) para refrigerar los aceros de baja calidad que Junkers se veía obliga- da a emplear. Aún así, la vida de estos motores solía ser de unas 50 horas. Dado que esta cifra era mayor que la vida media de un caza alemán, fue considerada aceptable para la producción en serie. Dicha refrigeración se ha modelado sangrando del compresor el 8 % del flujo admi- tido, siendo empleado el 5 % en la refrigeración de los estátores de la turbina (con lo cual dicho flujo devuelve trabajo en el rotor) y un 3 % en la refrigeración de este últi- mo, flujo que no devuelve trabajo. En cuanto a las características del motor, Hans von Ohain proporciona [2] los si- guientes datos: 2000 lb de empuje, 46.6 lb/s de flujo másico, relación de compresión de 3.14, TIT de 1427°F, 1.363 lb(lb/h) de consumo y eficiencias de: 78 % compresor, 95 % cámara de combustión y 79.5 % turbina. A la hora de calcular el ciclo, se han supuesto rendimientos del difusos y coeficiente de empuje de tobera del 98 %, rendimiento mecánico del 97 % y pérdida de presión en la cámara de combustión del 5 %. Estos valores se han ajustando mediante el software Gasturb de tal manera que el empuje calculado sea similar a los originales, que son los siguientes: E [N] Es p [m/s] TSF C [g/(kN s)] 8825 417.45 38.61 Tabla 2.1: Prestaciones del motor original 7 3 Modelo de gas perfecto Una vez establecida la metodología y el motor de referencia, el primer método de cálculo que vamos a emplear es aquél que asume que el gas es calóricamente perfec- to, esto es, que su calor específico es constante y no se ve afectado por cambios de temperatura o presión. 3.1 Análisis por estaciones La manera habitual de calcular este tipo de ciclos es realizar un análisis de cada estación del motor, empleando donde sea necesario las ecuaciones de acople de po- tencias entre ellas. 3.1.1 Difusor La temperatura de salida del difusor será la misma que la temperatura de parada del flujo libre; en este caso dado que no hay velocidad, será la misma: T02 = T00 = T0 � 1+ γ f −1 2 M 20 � = T0 (3.1) Dado que para caracterizar el difusor disponemos del factor de recuperación de presiones de parada, la presión total a la salida será pues: rd = p02 p00 ⇒ p02 = rd p00 = rd p0 � 1+ γ f −1 2 M 20 � γ f γ f −1 (3.2) Programando esta etapa en Mathematica: 1 M0 = V0/Sqrt[gf*R*T0]; 2 T02 = T0*(1 + ((gf - 1)/2)*M0^2); 3 p00 = p0*(1 + ((gf - 1)/2)*M0^2)^(gf/(gf - 1)); 4 If[etad == 0, 5 p02 = rd*p00 , 6 p02 = p0*(1 + etad *((gf - 1)/2)*M0^2)^(gf/(gf - 1)) 7 ]; 3.1.2 Compresor En el compresor se aumenta la presión del fluido, lo que lleva asociado un incre- mento de la temperatura según la ecuación de la compresión isentrópica y su rendi- miento: p03 =pic p02 T03 = T02 1+ pi γ f −1 γ f c −1 ηc  (3.3) 1 p03 = p02*pic; 2 T03 = T02*(1 + (pic^((gf - 1)/gf) - 1)/etac); 8 COMPARACIÓN MODELOS 3 MODELO DE GAS PERFECTO 3.1.3 Cámara de combustión A continuación, modelamos en primera aproximación la cámara de combustión como una adición de calor isobárica minorada con una pérdida de presión∆pc c : p04 = p03(1−∆pc c ) (3.4) El combustible necesario para alcanzar la T04 impuesta se averigua mediante un ba- lance de potencias entre la entrada y la salida, teniendo en cuenta que dicho proceso no es ideal mediante un rendimiento ηc c y que hemos sangrado para refrigeración dos flujos c1 y c2: (1− c1− c2)cp f T03+ηc c m˙ f H f = (1− c1− c2+ f )cp c T04 (3.5) 1 p04 = p03 (1 - Deltapcc); 2 f = Solve [(1 - c1 - c2) Cpf T03 + etacc f Hf == (1 - c1 - c2 + f) Cpc T04 , f][[1, 1, 2]]; 3.1.4 Turbina En primer lugar modelamos la mezcla isobárica del flujo de refrigeración del esta- tor, c1, preparando la posibilidad de que las propiedades del flujo cambien: cp m 1 = cp c (1+ f − c 1− c 2)+ cp f c1 1+ f − c2 γm 1 = cp m 1 cp m 1−R (1+ f − c1− c2)cp c T04+ c1cp f T03 = (1+ f − c2)cp m 1T041 (3.6) 1 Cpm1=(Cpc (1+f-c1-c2)+Cpf c1)/(1+f-c2);gm1=Cpm1/(Cpm1 -R); 2 T041=T041/. Flatten[Solve [(1+f-c1-c2)Cpc T04+ c1 Cpf T03 ==(1+f -c2)Cpm1 T041 ,T041 ]]; Realizando un balance de potencias en el eje que une compresor y turbina, tenien- do en cuenta un rendimiento ηm y que el flujo c1 inyectado en el estator genera tra- bajo: cp f (T03−T02) =ηm m˙a (1+ f − c2)cp m 1(T041−T045) p045 = p041 � 1− 1 ηt � 1− T045 T041 �� γc γc−1 (3.7) 1 T045 = T045 /. Flatten[Solve[Cpf*(T03 - T02) == etam *(1 + f - c2)*Cpm1*(T041 - T045), T045 ]]; 2 p045 = p04*(1 - (T041 - T045)/(T041*etat))^(gm1/(gm1 - 1)); Para acabar con la zona de la turbina, realizamos la segunda mezcla isobárica, esta vez del flujo de refrigeración del rotor, c2: cp m 2 = cp m 1(1+ f − c 2)+ cp f c2 1+ f γm 2 = cp m 2 cp m 2−R (1+ f − c2)cp m 1T045+ c2cp f T03 = (1+ f − c2)cp m 2T05 (3.8) 9 COMPARACIÓN MODELOS 3 MODELO DE GAS PERFECTO 1 Cpm2=(Cpm1 (1+f-c2)+Cpf c2)/(1+f);gm2=Cpm2/(Cpm2 -R); 2 T05=T05/. Flatten[Solve [(1+f-c2)Cpm1 T045+c2 Cpf T03 ==(1+f-c2) Cpm2 T05 ,T05]]; 3.1.5 Tobera Si la tobera está operando en diseño y expande a presión ambiente, teniendo en cuenta que la última mezcla se presume isobárica y utilizando la ecuación de la ex- pansión isentrópica ponderada con su eficiencia: p9 = p0 T9 = T05 1−ηn1−� p9 p045 � γc−1 γc  (3.9) 1 p9=p0; 2 T9=T05(1-etan(1-(p9/p045)^((gm2 -1)/gm2))); Donde se ha tenido en cuenta que a expansión completa, la eficiencia isentrópica es el cuadrado del coeficiente de empuje. Para obtener la velocidad de salida del flujo V9, debemos tener en cuenta la conversión de entalpía a energía cinética: V9 = p 2cp m 2(T05−T9) (3.10) 1 V9=Sqrt[2 Cpm2(T05 -T9)]; 3.2 Prestaciones y figuras de mérito Una vez recorrido el ciclo termodinámico, podemos calcular las prestaciones y fi- guras de mérito que nos servirán para caracterizarlo y comparar distintos modelos. 3.2.1 Empuje y empuje específico Planteando el balance de cantidad de movimiento y teniendo en cuenta que la to- bera expande hasta presión atmosférica: E = m˙a (1+ f )V9− m˙a V0 (3.11) Dividimos el empuje entre el flujo másico admitido para obtener el empuje especí- fico, que es una variable intensiva: Es = (1+ f )V9−V0 (3.12) 1 Es = (1 + f) V9 - V0; 2 Enet = Es*ma; 10 COMPARACIÓN MODELOS 3 MODELO DE GAS PERFECTO 3.2.2 Consumo específico Se trata de un cociente entre el flujo másico de combustible y el empuje. Se expresa usualmente en g/(kN s) o bien en lb/(lbf h); aquí se ha elegido lo primero: TSF C = m˙ f E = f Es (3.13) 1 TSFC=(f 1000)/(Es /1000); 3.2.3 Rendimiento térmico Es posible definirlo teniendo en cuenta la potencia desaprovechada en el escape: Pdes. escape = Ec V9 −Ec V0 = 12 (m˙a + m˙ f )(V9−V0) 2 (3.14) El rendimiento térmico asociado a esta definición será por tanto: ηt h = E V0+(Ec V9 −Ec V0 ) m˙ f H f = (1+ f )V9−V0+ 12 (1+ f )(V9−V0)2 f H f (3.15) 1 etath=(V0 Es+1/2 (1+f)(V9-V0)^2)/(f Hf); 3.3 Resultados del cálculo Empleando el modelo descrito con los datos del motor propuesto, y suponiendo las siguientes propiedades termodinámicas fijas para el gas: R = 287 J/(kg K) cp c = cp f = c t e .= 1004.5 J/kg (3.16) Se obtienen los siguientes resultados, que se discutirán más adelante junto con los del resto de modelos. Se aprecia no obstante que resultan una primera aproximación razonable a las características reales del motor: E [N] Es p [m/s] ηt h [-] TSF C [g/(kN s)] 8692.48 411.18 0.135 34.76 Tabla 3.1: Prestaciones con el modelo de propiedades constantes 11 4 Modelo de gas semiperfecto El siguiente modelo en cuanto a profundidad es aquél en que dejamos de asumir que el gas es calóricamente perfecto y pasamos a asumir que es térmicamente per- fecto o semiperfecto, de tal manera que su calor específico es una función exclusiva de la temperatura pero no de la presión. También tendremos en cuanta el cambio de R con el dosado. En este caso, la función dependiente de la temperatura la tomaremos como una regresión polinómica de séptimo grado, tanto para para el aire como para el com- bustible, y supondremos que ambas propiedades se escalan en la mezcla con el flujo másico de cada uno, así como una expresión para R basada en el dosado [2]. cp (T ) = A0+A1T +A2T 2+A3T 3+A4T 4+A5T 5+A6T 6+A7T 7 (4.1) R( f ) = 8314,17 28,97−0,946186 f (4.2) 1 R[f_] = 1.9857117/(28.97 - 0.946186 f ) 4187; Donde las constantes An son, para el caso del aire, cp a (T ), las siguientes: A0 = 2,5020051 10−1 A1 =−5,153687910−5 A2 = 6,5519486 10−8 A3 =−6,717837610−12 A4 =−1,512825910−14 A5 = 7,6215767 10−18 A6 =−1,452677010−21 A7 = 1,0115540 10−25 (4.3) 1 cpa[T_] = (A0 + A1 T + A2 T^2 - A3 T^3 + A4 T^4 + A5 T^5 + A6 T^6 - 2 A7 T^7) 4187; Mientras que para el caso del combustible, cp f (T ) usaremos las siguientes: A0 f = 7,3816638 10−2 A1 f = 1,2258630 10−3 A2 f =−1,377190110−6 A3 f = 9,9686193 10−10 A4 f =−4,205110410−13 A5 f = 1,0212913 10−16 A6 f =−1,333566810−20 A7 f = 7,2678710 10−25 (4.4) 1 cpf[T_] = (A0f + A1f T + A2f T^2 - A3f T^3 + A4f T^4 + A5f T ^5 + 2 A6f T^6 - A7f T^7) 4187; Donde hay que tener en cuenta que las constantes están expresadas en BTU/(lb ºR) y que hay que multiplicar por 4187 para expresarlas en el SI, J/(kg K). Por último, para facilitar el cálculo del ciclo, escribimos el cp del fluido en función de la proporción de aire y gas que presente: cp (T, a , f ) = a cp a (T )+ f cp f (T ) a + f (4.5) 1 cp[T_,a_,fv_]=(a cpa[T]+ fv cpf[T])/(a+fv); 12 COMPARACIÓN MODELOS 4 MODELO DE GAS SEMIPERFECTO 4.1 Análisis por estaciones Ahora ya estamos en condiciones de volver a recorrer el ciclo termodinámico del motor, tal como se hizo en el modelo anterior, usando esta vez sin embargo la expre- sión más general para el proceso isentrópico sin asumir cp constante. 4.1.1 Difusor Dado que nos encontramos en banco de ensayos, y por tanto no existe velocidad externa de entrada tenemos que: T02 = T00 = T0 (4.6) p02 = rd p00 = rd p0 (4.7) 1 T02=T0; p00=p0; p02=rd p00; 4.1.2 Compresor En este caso conocemos la relación de compresión y a través de ella hallamos la presión a la salida, que usamos con la ecuación general del cambio de entropía para averiguar la temperatura final ideal: p03 =pic p02 (4.8) ∆S = ∫ T03s T02 cp (T, 1, 0) d T T −R(0) ln � p03 p02 � = 0 (4.9) 1 p03=p02 pic; 2 DS=Integrate[cp[T,1,0] 1/T,{T,T02 ,T03s},Assumptions ->Re[T03s ]>=0]-R[0]Log[p03/p02]; 3 T03s = FindRoot[DS, {T03s , T02}][[1, 2]]; Una vez obtenida la temperatura ideal, empleamos la definición de rendimiento isentrópico del compresor para averiguar la temperatura real: ηc = W˙i W˙r = m˙a∆h i m˙a∆hr = cp (T03s , 1, 0)T03s − cp (T02, 1, 0)T02 cp (T03, 1, 0)T03− cp (T02, 1, 0)T02 (4.10) 1 T03=FindRoot [(cp[T03s ,1,0]T03s -cp[T02 ,1,0]T02)/(cp[T03 ,1,0] T03 -cp[T02 ,1,0]T02)-etac ,{T03 ,T03s }][[1 ,2]]; 4.1.3 Cámara de combustión Las ecuaciones son las mismas en este caso, con la precaución de mantener el cp en el balance de potencias y tener en cuenta el sangrado de refrigeración: p04 = p03(1−∆pc c ) (4.11) (1− c1− c2)cp (T03, 1− c1− c2, 0)T03+ηc c f H f = (1− c1− c2)cp (T04, 1− c1− c2, f )T04 (4.12) 13 COMPARACIÓN MODELOS 4 MODELO DE GAS SEMIPERFECTO 1 p04=p03(1-Dpcc); 2 f=Solve[(1-c1-c2)cp[T03 ,1-c1-c2 ,0] T03+ etacc f Hf == (1-c1- c2+f)cp[T04 ,1-c1-c2,f] T04 ,f][[5 ,1 ,2]; 4.1.4 Turbina En primer lugar, debemos añadir el flujo de refrigeración del estátor, que se mezcla a presión constante: (1+ f − c1− c2)cp (T04, 1− c1− c2, f )T04︸ ︷︷ ︸ flujo principal +c1cp (T03, c1, 0)T03︸ ︷︷ ︸ flujo refrigeración = (1+ f − c2)cp (T041, 1− c2, f )T041︸ ︷︷ ︸ flujo resultante (4.13) 1 T041=Solve [(1+f-c1-c2)cp[T04 ,1-c1-c2,f] T04+c1 cp[T03 ,c1 ,0] T03 ==(1+f-c2)cp[T041 ,1-c2,f] T041 ,T041 ][[5 ,1 ,2]]; A continuación resolvemos el rotor de la turbina, considerando la ecuación de aco- plamiento de potencia con el compresor, así como los rendimientos mecánicos y el isentrópico de la turbina, que mayoran el trabajo requerido por el compresor: ∆hc = ∫ T03 T02 cp (T, 1, 0)d T ∆h t = ∫ T041 T045 cp (T, 1− c2, f )d T ∆hc 1 ηmηt = (1+ f − c2)∆h t (4.14) 1 Dhc=Integrate[cp[T,1,0],{T,T02 ,T03}]; 2 Dht=Integrate[cp[T,1-c2,f],{T,T045 ,T041 }]; 3 T045=Solve[Dhc 1/( etam etat) == (1+f-c2)Dht ,T045 ][[5 ,1 ,2]]; Ahora necesitamos emplear la ecuación del salto entrópico para averiguar la pre- sión que queda tras esta expansión en el rotor: ∆S = ∫ T045 T041 cp (T, 1− c2, f )d T T −R( f ) ln � p045 p041 � = 0 (4.15) 1 DS = Integrate[cp[T, 1 - c2, f] 1/T, {T, T041 , T045}] - 2 R[f] Log[p045/p04]; 3 p045 = Solve[DS == 0, p045 ][[5, 1, 2]]; Por último, realizamos la segunda mezcla isóbarica del flujo de refrigeración del rotor, que por tanto no ha sido expandida por éste: (1+ f − c2)cp (T045, 1− c2, f )T045︸ ︷︷ ︸ flujo principal +c2cp (T03, c2, 0)T03︸ ︷︷ ︸ flujo refrigeración 2 = (1+ f )cp (T05, 1, f )T05︸ ︷︷ ︸ flujo resultante (4.16) 1 T05 = Solve [(1 + f - c2) cp[T045 , 1 - c2, f] T045 + c2 cp[T03 , c2, 0] T03 == (1 + f) cp[T05 , 1, f] T05 , T05][[5, 1, 2]]; 14 COMPARACIÓN MODELOS 4 MODELO DE GAS SEMIPERFECTO 4.1.5 Tobera La tobera se va a resolver de forma isentrópica ideal, debido a que su rendimiento se expresa mediante la relación de energía cinética a la salida, con lo que se aplicará al final de ésta. Sabiendo que la tobera está adaptada: p9 = p0 (4.17) ∆S = ∫ T9s T05 cp (T, 1, f ) d T T −R( f ) ln � p9 p045 � = 0 (4.18) 1 p9 = p0; 2 DS = Integrate[cp[T, 1, f] 1/T, {T, T05 , T9s}, 3 Assumptions -> Re[T9] >= 0] - R[f] Log[p9/p045]; 4 T9s=FindRoot[DS ,{T9,T045 }][[1 ,2]]; Ahora tenemos que calcular cuál es el salto entálpico asociado a este salto de pre- siones: ∆hs = ∫ T05 T9s cp (T, 1, f )d T (4.19) 1 Dhs = Integrate[cp[T, 1, f], {T, T9s , T05}]; Y por último transformar a velocidad dicho salto entálpico, usando el rendimiento isentrópico de la tobera para calcular la velocidad de salida real del flujo: V9s = p 2∆hs ηn = V 29 /2 V 29s /2 ⇒V9 = Æ ηn V 29s (4.20) 1 V9s = Sqrt[2 Dhs]; V9 = Sqrt[etan V9s ^2]; 4.2 Resultados del cálculo La manera de calcular las prestaciones del ciclo es exactamente la misma que la expuesta en el modelo anterior. Empleando las mismas fórmulas ya introducidas se obtienen los siguientes resultados: E [N] Es p [m/s] ηt h [-] TSF C [g/(kN s)] 7626.79 360.77 0.107 38.50 Tabla 4.1: Prestaciones con el modelo de propiedades variables 15 5 Modelo de equilibrio El siguiente modelo que vamos a considerar se basa en un modelo de gases reales, en base a la consideración de equilibrio químico y termodinámico entre las distintas especies presentes en el turborreactor a las condiciones en las que se encuentra en cada estación. Para ello emplearemos el programa CEA (Chemical Equilibrium with Applications) creado por la NASA [3]. De nuevo, volveremos a recorrer el ciclo termodinámico es- tación a estación, aplicando los mismos parámetros de eficiencia. El combustible se- leccionado es JP-4. 5.1 Análisis por estaciones Durante este análisis se indicará qué problemas se resuelven mediante CEA, ad- juntando el input del programa y los resultados que ofrece, y realizando los cálculos auxiliares con Mathematica. Los cálculos se realizarán en K y atm para facilitar el input. Las entalpías y entropías se van a dar en kJ/kg y kJ/(kg K) respectivamente, pero a la hora de introducirlas en el CEA se deben dividir por Ru ni v . 5.1.1 Difusor En esta primera estación se aplica el coeficiente de recuperación de presiones de parada y se resuelve un programa TP de CEA para calcular las propiedades a la salida: T02 = T0 = 288,15 p02 = rd p0 = 0,98 (5.1) 1 problem case =1234 2 tp t,k=288.15 , p,atm=0.98, 3 react 4 oxid=Air wt=100 5 end Con lo que se obtiene S02 = 6,832 y h02 =−14,379. Un ejemplo de la salida completa del programa CEA se ofrece anexa. 5.1.2 Compresor El proceso en el compresor es similar: conocemos la presión a la salida y asumimos isentropía, entrando al CEA con ambos datos: p03 =pic p02 = 3,0772 S03 =S02 = 6,832 (5.2) 16 COMPARACIÓN MODELOS 5 MODELO DE EQUILIBRIO 1 problem case =1234 2 sp p,atm =3.0772 , t,k=800 s/r=0.821746 3 react 4 oxid=Air wt=100 5 end De esta manera obtenemos la entalpía ideal h03s = 129,337 y aplicando la defini- ción de rendimiento isentrópico obtenemos la real. Entrando al CEA con dicho dato y resolviendo un programa HP: ηc = h03s −h02 h03−h02 ⇒ h03 = 129,33 (5.3) 1 problem case =1234 2 hp p,atm =3.0772 , h/r=15.5566 t,k=800 3 react 4 oxid=Air wt=100 t,k=288.15 5 end Obtenemos que la temperatura es T03 = 430,5. 5.1.3 Cámara de combustión Sabemos que la presión disminuye debido a las pérdidas, p04 = (1−∆pc c )p03 = 2,92334; después de calcularla realizamos un proceso iterativo. Se supone un dosa- do de combustible y conociendo la temperatura final T04 y la presión, se calcula un problema TP para obtener una entalpía, que se introduce en la ecuación del balance energético (donde se tiene en cuenta el calor de formación del JP-4, -22.72 kJ/kg) para obtener un dosado que se realimenta: f a r ′4 = h03−h04( f a r4) h04( f a r4)−Ho −ηc c H f (5.4) 1 problem case =1234 %fuel =0.01681544 , 2 tp t,k=1048, p,bar =2.923340 , 3 react 4 oxid=Air wt=100 5 fuel=JP -4 wt=100 6 end La iteración converge rápidamente, siendo los valores: f a r4 supuesto h04 de CEA f a r 4′ del balance ∆ f a r4 ( %) 0,02 789,93 0,016778766 19,19827852 0,016778766 791,35 0,01681544 -0,218096628 0,01681544 791,33 0,016814923 0,003071987 Tabla 5.1: Resultado de la iteración en CC 17 COMPARACIÓN MODELOS 5 MODELO DE EQUILIBRIO 5.1.4 Turbina En primer lugar, realizamos la mezcla a presión constante, p041 = p04, mediante la ecuación de balance: c1h03+(1− c2− c2+ f )h04 = (1− c2+ f )h041⇒ h041 = 757,78 (5.5) Teniendo el cuenta que el f a r en esta estación hay que recalcularlo, podemos en- trar en CEA y resolver un problema HP para hallar la temperatura, que resulta ser T041 = 1018,65 f a r41 = f 1− c2− c2 1− c1 = 0,016283 (5.6) 1 problem case =1234 %fuel =0.016283 , 2 hp p,atm =2.92334 , h/r=91.146 t,k=1000 3 react 4 oxid=Air wt=100 t,k=298.15 5 fuel=JP -4 wt=100 t,k=298.15 6 end Ahora calculamos, mediante el acoplamiento de potencia del compresor, la ental- pía necesaria: (1− c2+ f )(h014−h045) = 1 ηm (h03−h02)⇒ h045 = 612,151 (5.7) Aplicando la definición de rendimiento isentrópico dela turbina podemos hallar la entropía isentrópica, que mediante un proceso iterativo en CEA (resolviendo el pro- blema HP) nos conducen a hallar la presión. ηt = h041−h045 h041−h045s ⇒ h045s = 574,597 (5.8) 1 problem case =1234 %fuel =0.016283 , 2 hp p,atm =1.428 , h/r=69.112 t,k=1000 3 react 4 oxid=Air wt=100 t,k=298.15 5 fuel=JP -4 wt=100 t,k=298.15 6 end Para ello se ha entrado con la presión calculada mediante polinomios, p045 = 1,428, que ofrece un error en la entalpía calculada por CEA del 0.0012 % (de 7.8575 a 7.8477), considerado aceptable. El último paso en la turbina consiste en calcular la refrigeración del rotor, conside- rada isobárica (p05 = p045 = 1,428) mediante el balance, y resolver en CEA el problema HP para hallar las condiciones finales, teniendo en cuenta que reunido todo el flujo, ahora f a r = f . (1− c2+ f )h045+ c2h03 = (1+ f )h05⇒ h05 = 597,906 (5.9) 1 problem case =1234 %fuel =0.016814 , 2 hp p,atm =1.428 , h/r=71.9156 t,k=1000 3 react 4 oxid=Air wt=100 t,k=298.15 5 fuel=JP -4 wt=100 t,k=298.15 6 end Con lo que tenemos que se sale de la turbina con T05 = 877,57 y S05 = 7,8847. 18 COMPARACIÓN MODELOS 5 MODELO DE EQUILIBRIO 5.1.5 Tobera En la tobera tenemos en cuenta que se conservan en principio las entalpías totales, h09 = h05 = 597,906, las entropías S9 = S09 = S05 = 7,8847 y que la tobera está adap- tada, p9 = p0 = 1, lo que resolviendo un problema SP nos permite hallar la entalpía disponible tras la expansión: 1 problem case =1234 %fuel =0.016814 , 2 sp p,atm=1, t,k=1000 s/r=0.948364 3 react 4 oxid=Air wt=100 t,k=298.15 5 fuel=JP -4 wt=100 t,k=298.15 6 end El cálculo arroja un resultado de h9s = 512,56 y por tanto una diferencia de entalpías de h09 − h9s = 85,246 kJ/kg lo que teniendo en cuenta la ecuación del balance de energías a la salida: h9s + V 29s 2 = h09⇒V9s = 413.149 m/s (5.10) Aplicando ahora la definición de rendimiento isentrópico de la tobera, podemos calcular la velocidad de salida real: ηn = V 29 V 29s ⇒V9 = 404.886 m/s (5.11) 5.2 Resultados del cálculo Ahora ya podemos emplear las definiciones anteriores de prestaciones y figuras de mérito del ciclo termodinámico, con los datos calculados mediante este modelo de equilibrio, que arrojan los siguientes resultados: E [N] Es p [m/s] ηt h [-] TSF C [g/(kN s)] 8703.2 411.694 0.115 40.841 Tabla 5.2: Prestaciones con el modelo de equilibrio 19 6 Software de cálculo El último paso en cuanto a modelos consiste en emplear un software diseñado es- pecíficamente para el cálculo y el análisis de ciclos termodinámicos de aerorreacto- res, como por ejemplo GasTurb 11. Empleando dicho software se pueden calcular prestaciones, realizar optimizacio- nes, estudios paramétricos, modificar la geometría de los componentes, manipular mapas de compresor y turbina, realizar análisis off-design, etc. En nuestro caso vamos a introducir los datos de nuestro motor, especificando que vamos a realizar un cálculo de turborreactor de flujo único en punto de diseño, intro- duciendo además todos los datos que tenemos de nuestro motor. Figura 6.1: Ventana de datos principales de GasTurb Tal y como se ha comentado en la presentación del motor, para obtener el empuje nominal con las eficiencias dadas por von Ohain, ha sido necesario minorar los rendi- mientos mecánicos, de difusor, tobera, etc., proceso se realizará realizado en GasTurb. Cabe señalar también que GasTurb no dispone de la opción de especificar tobe- ra adaptada, así que se ha empleado la opción de tobera convergente-divergente de geometría variable y se ha lanzado un estudio de optimización del empuje variando la relación de áreas A9/A8. En este caso, dado que el flujo es subsónico, la tobera convergente expande total- mente el fluido hasta la presión ambiente, con lo que no es necesario en realidad añadir una parte convergente. El programa reacciona pues señalando que el óptimo se encuentra en A9/A8 = 1. 20 COMPARACIÓN MODELOS 6 SOFTWARE DE CÁLCULO Figura 6.2: Ventana de optimización de GasTurb Una vez establecidos los parámetros del cálculo y realizado el estudio de la situa- ción de la tobera, ya podemos ejecutar la opción de cálculo de ciclo. Lo resultados ofrecidos por GasTurb son los siguientes: Turbojet Alt= 0m ISA W T P WRstd Station kg/s K kPa kg/s FN = 8,85 kN amb 288,15 101,325 TSFC = 37,5188 g/(kN*s) 1 21,139 288,15 101,325 FN/W2 = 418,48 m/s 2 21,139 288,15 99,299 21,570 3 21,139 430,15 311,797 8,393 Prop Eff = 0,0000 31 19,448 430,15 311,797 eta core = 0,1326 4 19,779 1048,00 293,089 13,041 41 20,836 1019,23 293,089 13,548 WF = 0,33190 kg/s 49 20,836 891,24 146,271 s NOx = 0,05959 5 21,470 878,46 146,271 25,969 XM8 = 0,7565 6 21,470 878,46 146,271 A8 = 0,1158 m² 8 21,470 878,46 146,271 25,969 P8/Pamb = 1,4436 Bleed 0,000 430,15 311,797 WBld/W2 = 0,00000 -------------------------------------------- Ang8 = 0,00 ° P2/P1 = 0,9800 P4/P3 = 0,9500 P6/P5 1,0000 CD8 = 1,0000 Efficiencies: isentr polytr RNI P/P W_NGV/W2 = 0,05000 Compressor 0,7800 0,8118 0,980 3,140 WCL/W2 = 0,03000 Burner 0,9500 0,940 Loading = 100,00 % Turbine 0,7950 0,7803 0,655 2,004 e45 th = 0,78458 -------------------------------------------- far7 = 0,01570 Spool mech Eff 0,9700 Nom Spd 8700 rpm PWX = 0,00 kW -------------------------------------------- Con-Di Nozzle: A9/A8 = 1,00000 A9*(Ps9-Pamb) -1,01E-6 CFGid = 0,98000 -------------------------------------------- hum [%] war0 FHV Fuel 0,0 0,00000 43,100 Generic Tabla 6.1: Resultados del ciclo También es posible obtener el diagrama de entropía contra entalpía del ciclo que hemos modelado, en el que se observan las mezclas de la refrigeración a presión cons- tante, la expansión completa a temperatura ambiente, así como la notable ineficien- 21 COMPARACIÓN MODELOS 6 SOFTWARE DE CÁLCULO Figura 6.3: Nomenclatura y circuito de aire en GasTurb cia y baja relación de compresión del compresor, propia de un modelo de turborreac- tor ciertamente primitivo y diseñado en circunstancias extremadamente adversas. 12/11/2011 GasTurb 11 -.2 0 .2 .4 .6 .8 1 1.2 Entropy [kJ/(kg K)] -200 0 200 400 600 800 1000 E nt ha lp y [k J/ kg ] 50.0 k Pa 150 kPa 200 kPa 2 50 kP a Pso (a mbien t) = 1 01 kP a 0 2 3 441 49 5 8 s9 Figura 6.4: Diagrama del ciclo TS Resumiendo los resultados obtenidos para facilitar su comparación con los otros modelos, podemos reunir las siguientes prestaciones del ciclo termodinámico: E [N] Es p [m/s] ηt h [-] TSF C [g/(kN s)] 8850 418.48 0.1316 37.52 Tabla 6.2: Prestaciones con el modelo de GasTurb 22 7 Variaciones de los modelos Una vez planteado el funcionamiento de los distintos modelos, puede resultar útil introducir variaciones para analizar si es posible mejorar sus resultados. 7.1 Modelo de gas perfecto con 2 constantes Para refinar el modelo más simple sin llegar a realizar el cálculo con variables fun- ción de la temperatura, se pueden considerar dos valores constantes para el calor específico, extraídos de los resultados del CEA (promediando): cp f = 1010 J/(kg K) cp c = 1125 J/(kg K) (7.1) Introduciendo pues los datos anteriores en el modelo que ya se había programado, hemos obtenido los siguientes resultados: E [N] Es p [m/s] ηt h [-] TSF C [g/(kN s)] 8861.68 419.19 0.116 41.09 Tabla 7.1: Prestaciones con el modelo de propiedades constantes (2 val.) 7.2 Modelo de gas semiperfecto (correlación 2) Una manera de variar el modelo de gas semiperfecto consiste en probar diferentes correlaciones. Otra de ellas [4] es la que consiste en representar el calor específico y la constante R de la mezcla como: R( f a r ) = 287,05−0,00990 f a r +10−7 f a r 2 (7.2) cp (T, f a r ) = ‚ A0+A1 � T 1000 � +A2 � T 1000 �2 +A3 � T 1000 �3 +A4 � T 1000 �4 +A5 � T 1000 �5 +A6 � T 1000 �6 +A7 � T 1000 �7 +A8 � T 1000 �8 + f a r 1+ f a r ‚ B0+ B1 T 1000 + B2 � T 1000 �2 + B3 � T 1000 �3 + B4 � T 1000 �4 + B5 � T 1000 �5 + B6 � T 1000 �6 + B7 � T 1000 �7ŒŒ 103; (7.3) Donde las constantes empleadas son: A0 = 0,992313, A1 = 0,236688, A2 =−1,852148, A3 = 6,083152, A4 = −8,893933, A5 = 7,097112, A6 = −3,234725, A7 = 0,794571, A8 = −0,081873, A9 = 0,422178, A10 = 0,001053, B0 = −0,718874, B1 = 8,747481, B2 = −15,863157, B3 = 17,254096, B4 = −10,233795, B5 = 3,081778, B6 = −0,361112, B7 = −0,003919, B8 = 0,0555930, B9 =−0,0016079, obteniendo: E [N] Es p [m/s] ηt h [-] TSF C [g/(kN s)] 8837.05 418.02 0.106 44.86 Tabla 7.2: Prestaciones con la 2ª correlación. 23 8 Discusión Una vez realizado el cálculo según los distintos modelos termodinámicos de cálcu- lo, procede realizar una pequeña discusión de los datos obtenidos. Vamos a centrar- nos en cómo estos modelos predicen el empuje específico y el consumo específico respecto al empuje Es p [m/s] TSF C [g/(kN s)] Gas perfecto 405,060 35,286 Gas perfecto 2 val. 419,190 41,095 Gas semiperfecto 1 391,160 355,175 Gas semiperfecto 2 418,025 44,859 Equilibrio (CEA) 411,690 40,841 GasTurb 418,480 375,188 Motor real 417,450 38,610 Tabla 8.1: Resultados absolutos de los distintos modelos A la hora de realizar una comparación, y ya que disponemos de los datos del motor real que hemos intentado modelizar, lo ideal es expresar los datos anteriores como variaciones porcentuales de esta referencia: ∆Es p [%] ∆TSF C [%] Gas perfecto -2,968 -8,609 Gas perfecto 2 val. 0,417 6,436 Gas semiperfecto 1 -6,298 -8,010 Gas semiperfecto 2 0.137 16.185 Equilibrio (CEA) -1,380 5,778 GasTurb 0,247 -2,826 Tabla 8.2: Resultados relativos a los datos del motor real Perf . Perf . 2 val Semip. 1 Semip. 2 CEA GT -5 5 10 15 DEsp @%D DTSFC @%D Figura 8.1: Resultados relativos a los datos del motor real 24 COMPARACIÓN MODELOS 8 DISCUSIÓN Se puede observar que el principal problema suele ser la predicción acertada del TSF C , siendo los más exactos en este apartado el modelo de CEA y el GasTurb, con aproximadamente un 3 % de diferencia respecto al motor real. La diferencia máxima es de un 16 %, empleando la segunda correlación polinómica. En lo que respecta a la predicción del empuje específico, resulta sorprendente que el método de los dos calores específicos constantes resulte muy acertado. No lo es tanto como GasTurb, pero esto se debe a que, tal como se ha explicado, con GasTurb se han ajustado los parámetros de ineficiencias desconocidos. En este apartado también es muy preciso el CEA, con un 1 % de diferencia respeto al motor real, y la segunda correlación polinómica, con apenas un 0.13 %. Se observa sin embargo una amplia variabilidad según la correlación que estemos considerando, lo que demuestra que se trata de unos datos críticos para este modelo. Por último se grafica la diferencia en la predicción de rendimiento térmico del ciclo, en este caso respecto al obtenido mediante GasTurb, ya que no se ha hallado la figura del motor real. Pe rf . Pe rf . 2 v al Se mi p. 1 Se mi p. 2 CE A -12 -10 -8 -6 -4 -2 0 D % re sp ec to a G as Tu rb Figura 8.2: Comparación de ηt h respecto a GasTurb En este caso se observa que el método de gas perfecto se aproxima mucho a Gas- Turb, mientras que el CEA predice un resultado bastante diferente. Por supuesto, se ha de tener en cuanta que estamos comparando respecto a GasTurb, no respecto al motor real, con lo que este dato simplemente ilustra la diferencia entre los distintos modelos. 25 9 Conclusiones A la vista de los datos y la discusión anterior, se observa que el método más con- sistente en cuanto a empuje, consumo específico y rendimiento térmico es el de gas semiperfecto empleando la 2ª correlación. Se puede concluir también que la modelización de la combustión muy problemáti- ca, probablemente debido a que las ineficiencias y características del motor original son difíciles de modelar de manera analítica. También es de resaltar que el modelo de gas perfecto considerando dos calores se- gún el flujo sea frío o caliente ofrece una aproximación muy buena comparada res- pecto a la sencillez de su planteamiento y su velocidad de ejecución. En definitiva, la comparación respecto a datos reales del motor en banco subraya la importancia de la obtención de datos precisos del motor real, como las ineficiencias (mecánicas, pérdidas de presión, etc.) o los rendimientos (compresor, turbina, tobe- ra, etc.). Por muy sofisticado que sea el modelo termodinámico que empleemos, si estos da- tos no son muy exactos, no logramos obtener un aumento significativo de la precisión respecto al motor físico, porque nos limita dicha estimación, no el modelo termodi- námico en sí. A modo de ejemplo, no resultaría eficiente realizar para un programa de actuacio- nes el cálculo mediante equilibrio químico de 200 especies diferentes, si en realidad no se tienen datos del fabricante o de ensayos en banco y los parámetros se han ob- tenido de modo estimativo o mediante correlaciones de madurez tecnológica, como se hace para los problemas académicos. Referencias [1] FOSTER, J., 1945. Messerschmitt Me-262 Jet Fighter. Part II - The Power Plant., Avia- tion, Volume 44. New York: McGraw-Hill Publishing Company. [2] MATTINGLY, J.D., VON OHAIN, H. and AMERICAN INSTITUTE OF AERONAU- TICS AND ASTRONAUTICS, 2006. Elements of propulsion: gas turbines and rockets. Reston, Virginia: AIAA. [3] MCBRIDE, B.J. and GORDON, S.,1996. Computer Program for Calculation of Com- plex Chemical Equilibrium Compositions and Applications II. Users Manual and Program Description.. NASA RP 1311. [4] WALSH, P.P. and FLETCHER, P., 2004. Gas turbine performance. 2 edn. Oxford: Blackwell. 26 Apéndice A Ejemplo de salida CEA A modo de muestra, se anexa el resultado de CEA para el difusor, donde se puede apreciar las especies químicas consideradas. ******************************************************************************* NASA-GLENN CHEMICAL EQUILIBRIUM PROGRAM CEA2, MAY 21, 2004 BY BONNIE MCBRIDE AND SANFORD GORDON REFS: NASA RP-1311, PART I, 1994 AND NASA RP-1311, PART II, 1996 ******************************************************************************* problem case=1234 tp t,k=288.15, p,atm=0.98, react oxid=Air wt=100 end OPTIONS: TP=T HP=F SP=F TV=F UV=F SV=F DETN=F SHOCK=F REFL=F INCD=F RKT=F FROZ=F EQL=F IONS=F SIUNIT=T DEBUGF=F SHKDBG=F DETDBG=F TRNSPT=F T,K = 288.1500 TRACE= 0.00E+00 S/R= 0.000000E+00 H/R= 0.000000E+00 U/R= 0.000000E+00 P,BAR = 0.992985 REACTANT WT.FRAC (ENERGY/R),K TEMP,K DENSITY EXPLODED FORMULA : Air 1.000000 0.000000E+00 0.00 0.0000 N 1.56168 O 0.41959 AR 0.00937 C 0.00032 SPECIES BEING CONSIDERED IN THIS SYSTEM (CONDENSED PHASE MAY HAVE NAME LISTED SEVERAL TIMES) LAST thermo.inp UPDATE: 9/09/04 g 3/98 *Ar g 7/97 *C g 8/99 *CN g12/99 CNN tpis79 *CO g 9/99 *CO2 tpis91 *C2 g 7/00 CCN tpis91 CNC srd 01 OCCN tpis79 C2N2 g 8/00 C2O tpis79 *C3 srd 01 CNCOCN g 7/88 C3O2 g tpis *C4 g 6/01 C4N2 g 8/00 *C5 g 5/97 *N g 6/01 NCO tpis89 *NO g 4/99 NO2 j12/64 NO3 tpis78 *N2 g 6/01 NCN g 4/99 N2O g 4/99 N2O3 tpis89 N2O4 g 4/99 N2O5 tpis89 N3 g 5/97 *O tpis89 *O2 g 8/01 O3 n 4/83 C(gr) n 4/83 C(gr) n 4/83 C(gr) O/F = 0.000000 27 EFFECTIVE FUEL EFFECTIVE OXIDANT MIXTURE ENTHALPY h(2)/R h(1)/R h0/R (KG-MOL)(K)/KG 0.00000000E+00 0.00000000E+00 0.00000000E+00 KG-FORM.WT./KG bi(2) bi(1) b0i *N 0.53915890E-01 0.53915890E+01 0.53915890E-01 *O 0.14486046E-01 0.14486046E+01 0.14486046E-01 *Ar 0.32331996E-03 0.32331996E-01 0.32331996E-03 *C 0.11013248E-04 0.11013248E-02 0.11013248E-04 POINT ITN T N O AR C SINGULAR MATRIX, ITERATION 15 VARIABLE 4(EQLBRM) 1 31 288.150 -11.651 -13.123 -23.303 -171.775 THERMODYNAMIC EQUILIBRIUM PROPERTIES AT ASSIGNED TEMPERATURE AND PRESSURE CASE = 1234 REACTANT WT FRACTION ENERGY TEMP (SEE NOTE) KJ/KG-MOL K Air 1.0000000 0.000 0.000 O/F= 0.00000 %FUEL=100.000000 R,EQ.RATIO= 0.001521 PHI,EQ.RATIO= 0.000000 THERMODYNAMIC PROPERTIES P, BAR 0.99299 T, K 288.15 RHO, KG/CU M 1.2005 0 H, KJ/KG -14.379 U, KJ/KG -97.093 G, KJ/KG -1983.01 S, KJ/(KG)(K) 6.8320 M, (1/n) 28.965 (dLV/dLP)t -1.00000 (dLV/dLT)p 1.0000 Cp, KJ/(KG)(K) 1.0043 GAMMAs 1.4002 SON VEL,M/SEC 340.3 MOLE FRACTIONS *Ar 0.00937 *CO2 0.00032 *N2 0.78084 *O2 0.20948 * THERMODYNAMIC PROPERTIES FITTED TO 20000.K PRODUCTS WHICH WERE CONSIDERED BUT WHOSE MOLE FRACTIONS WERE LESS THAN 5.000000E-06 FOR ALL ASSIGNED CONDITIONS *C *CN CNN *CO *C2 CCN CNC OCCN C2N2 C2O *C3 CNCOCN C3O2 *C4 C4N2 *C5 *N NCO *NO NO2 NO3 NCN N2O N2O3 N2O4 N2O5 N3 *O O3 C(gr) NOTE. WEIGHT FRACTION OF FUEL IN TOTAL FUELS AND OF OXIDANT IN TOTAL OXIDANTS Apéndice B Ejemplo de ciclo con Mathematica A continuación se anexa el código que resuelve el ciclo termodinámico mediante el método que asume gas perfecto y por tanto calor específico constante. 1 Clear["‘*"] 3 (* Propiedades del motor *) 4 rd =0.98; (* Difusor *) 5 etad =0; 6 pic =3.14; (* Compresor *) 7 etac =0.78; 8 c1 =0.05; 9 c2 =0.03; 10 Deltapcc =0.05; (*CC*) 11 etacc =0.95; 12 Hf =43100000; 13 T04 =1048; 14 etam =0.97; (* Turbina *) 15 etat =0.795; 16 etan =0.98^2; (* Tobera *) 17 ma =21.14; 19 (* Propiedades del gas *) 20 R=287; 21 Cpf =1004.5; 22 Cpc =1004.5; 23 gf=Cpf/(Cpf -R); 24 gc=Cpc/(Cpc -R); 26 (* Propiedades del ambiente *) 27 T0 =288.15; 28 p0 =101325; 29 V0=0; 31 (*Ciclo termodinamico *) 33 (* Difusor *) 34 M0=V0/Sqrt[gf R T0]; 35 T02=T0(1+(gf -1)/2 M0^2); 36 p00=p0 (1+(gf -1)/2 M0^2)^(gf/(gf -1)); 37 If[etad==0,p02=rd p00 ,p02=p0 (1+ etad (gf -1)/2 M0^2)^(gf/(gf -1))]; 39 (* Compresor *) 40 p03=p02 pic; 41 T03s=T02 pic^((gf -1)/gf); 42 T03=T02+(T03s -T02)/etac; 43 T03=T02 (1+( pic^((gf -1)/gf) -1)/etac); 45 (* Camara de combustion *) 46 p04=p03(1-Deltapcc); 30 47 f=Solve[(1-c1-c2)Cpf T03+ etacc f Hf == (1-c1-c2+f)Cpc T04 ,f ][[1 ,1 ,2]]; 49 (* Turbina *) 50 (* Refrigeracion NGV*) 51 Cpm1=(Cpc (1+f-c1-c2)+Cpf c1)/(1+f-c2); 52 gm1=Cpm1/(Cpm1 -R); 53 T041=T041/. Flatten[Solve [(1+f-c1-c2)Cpc T04+c1 Cpf T03 ==(1+f- c2)Cpm1 T041 ,T041]] 54 1016.66 56 (*Rotor*) 57 T045=T045/. Flatten[Solve[Cpf (T03 -T02)==etam (1+f-c2) Cpm1 ( T041 -T045),T045 ]]; 58 p045=p04 (1-(T041 -T045)/(T041 etat))^(gm1/(gm1 -1)) 59 144686. 61 (* Refrigeracion rotor*) 62 Cpm2=(Cpm1 (1+f-c2)+Cpf c2)/(1+f); 63 gm2=Cpm2/(Cpm2 -R); 64 T05=T05/. Flatten[Solve [(1+f-c2)Cpm1 T045+ c2 Cpf T03 ==(1+f) Cpm2 T05 ,T05]] 65 854.14 67 (* Tobera *) 68 p9=p0; 69 T9=T05(1-etan(1-(p9/p045)^((gm2 -1)/gm2))); 70 V9=Sqrt[2 Cpm2(T05 -T9)]; 72 (* Resultados *) 73 Es=(1+f)V9-V0; 74 etath=(V0 Es+1/2 (1+f)(V9-V0)^2)/(f Hf); 75 TSFC=(f 1000)/(Es /1000); Nomenclatura Introducción Objetivo Metodología Motor de referencia Modelo de gas perfecto Análisis por estaciones Prestaciones y figuras de mérito Resultados del cálculo Modelo de gas semiperfecto Análisis por estaciones Resultados del cálculo Modelo de equilibrio Análisis por estaciones Resultados del cálculo Software de cálculo Variaciones de los modelos Modelo de gas perfecto con 2 constantes Modelo de gas semiperfecto (correlación 2) Discusión Conclusiones Referencias Ejemplo de salida CEA Ejemplo de ciclo con Mathematica