Matlab Turbojet

%
% ANEXO 1 : RUTINA (MATLAB) DE CALCULO (ANALISIS PARAMETRICO DE UN TURBOJET
% REAL
%
%
%
%---------------------------------------------------------------
%Analisis parametrico de un turbojet real
%Autor: Pablo Cornejo
%Propulsion 2006%Unidades: Pa K Kgr m s
%---------------------------------------------------------------
%Entradas:
%---------------------------------------------------------------
clc
clear all
Mo=0.85; % Numero de Mach de vuelo
go=9.8; % aceleracion degravedad nivel de mar
m=613000; % Masa avion
n_mot=4; % numero de motores del avion
mfuel_min=0.357; % razon masa de combustible masa al inicio
Cl_Cd=20; % razon de planeo
o_e=0.02; % coeficiente resistencia aerodinamica motores
To=223; % Temperatura de la corrientelibre 10000 m
Po=26490 ; % Presion de la corriente libre 10000 m
ga_c=1.4; % coeficiente poliprotico en el compresor
ga_t=1.3; % coeficiente poliprotico en la turbina
Cpc=1005; % poder calorifico gas en compresor
Cpt=1240; % poder calorifico gas en turbina
hpr=4.28e7;% entalpia especifica combustible
Pidmax=0.95; % relacion de presion difusor de entrada
Pib=0.92; % relacion de presion camara de combustion
Pin=0.93; % relacion de presion tobera
ec=0.84; % eficiencia poliprotica compresor
et=0.85; % eficiencia poliprotica turbinaNb=0.91; % eficiencia camara de combustion
Nm=0.99; % eficiencia transmicion potencia turbina-compresor
Pi_out=1; % turbina perfectamente expandida (P0/P9= = 1
Tt4=1130; % temperatura total salida camara de combustion
Pic=12; % relacion de presion compresor
P0_9=1; % latobera esta perfectamente expandida
Mi=0.5; % numero de Mach a la entrada
S_ideal=0.00002312; %consumo espeficico ideal
%--------------------------------------------------------------------
%Desarrollo de las ecuaciones
%--------------------------------------------------------------------%Empuje crucero requerido

Fctotal=(m*go)/((1-o_e)*Cl_Cd);
Fc=Fctotal/4;

%Empuje especifico y consumo especifico
Rc=(ga_c-1)*(Cpc/ga_c);
Rt=(ga_t-1)*(Cpt/ga_t);

a0=sqrt(ga_c*Rc*To);
V0=a0*Mo;Taor=1 + (ga_c-1)*(Mo^2)/2;
Pir=Taor^(ga_c/(ga_c-1));

Nr=1 ; % Mo < 1
Pid=Pidmax;

Taolam=(Cpt*Tt4)/(Cpc*To);

Taoc=Pic^((ga_c-1)/(ga_c*ec));
Nc=(Pic^((ga_c-1)/ga_c)-1)/(Taoc-1);

f=(Taolam-Taor*Taoc)/((hpr*Nb/(Cpc*To))-Taolam);

Taot= 1 - [Taor*(Taoc-1)/(Nm*(1+f)*Taolam)];

Pit=Taot^(ga_t/((ga_t-1)*et));
Nt=(1-Taot)/(1-Taot^(1/et));

Pt_9_9=P0_9*Pir*Pid*Pic*Pib*Pit*Pin;

M9=sqrt((2/(ga_t-1))*[(Pt_9_9)^((ga_t-1)/ga_t)-1]);...