clc; clear; close all;
format long
Ma=.84
Ra=.287
Rg=.287
siga=1.4
sigg=1.3333
Pa=54.05
Ta=225.7
OPR=2:3:200
tsfc=cell(1,67) ; %PRE-ALLOCATION
for i = 1:numel(OPR)
Va=Ma*sqrt(siga*Ra*Ta)
Cpoa= (siga*Ra)./(siga - 1)
TINT = 1100:100:1800
PRC = sqrt(OPR(i))
% Diffuser
efd=.93
To1 = Ta*(1+((siga-1)/2)*(Ma^2))
Toa=To1
Po1 = Pa*(1 + efd*((siga - 1)/2)*(Ma^2))^(siga/siga-1)
Poa=Pa*(Toa/Ta)^(siga/siga-1)
%Compressors
efc=.87
%Low Pressure
Po2=Po1.*(PRC)
To2 = To1 + ((To1*((PRC.^((siga-1)/siga)-1)))/efc)
WC_LP = Cpoa*(To1-To2)
%High Pressure
To3 = To2 + ((To2.*((PRC.^((siga-1)/siga)-1)))/efc)
Po3=Po2.*(PRC)
WC_HP = Cpoa*(To2-To3)
%Combustion
efb =.98
deltPob = .4
Po4= (1-deltPob)*Po3
To4=TINT
%Combustion Calculations for f
MC=14.4
MH=24.9
MO=0
Ycc=MC + (MH/4)- (MO/2)
Tp=To4
Tr=To3
To= 298
HRP_to = -8561991.6
if Tp <= 1600
hco2p = 56835 + 66.27*Tp + -11634.0*log(Tp)
hh2op = 88923 + 49.36*Tp + -7940.8*log(Tp)
hn2p= 31317 + 37.46*Tp + -4559.3*log(Tp)
ho2p = 43388 + 42.27*Tp + -6635.4*log(Tp)
else
hco2p = 93048 + 68.58*Tp + -16979.0*log(Tp)
hh2op = 154670 + 60.43*Tp + -19212.0*log(Tp)
hn2p = 44639 + 39.32*Tp + -6753.4*log(Tp)
ho2p =127010 + 46.25*Tp + -18798*log(Tp)
end
if Tr <= 1600
hco2r = 56835 + 66.27*Tr + -11634.0*log(Tr)
hh2or = 88923 + 49.36*Tr + -7940.8*log(Tr)
hn2r= 31317 + 37.46*Tr + -4559.3*log(Tr)
ho2r = 43388 + 42.27*Tr + -6635.4*log(Tr)
else
hco2r = 93048 + 68.58*Tr + -16979.0*log(Tr)
hh2or = 154670 + 60.43*Tr + -19212.0*log(Tr)
hn2r = 44639 + 39.32*Tr + -6753.4*log(Tr)
ho2r =127010 + 46.25*Tr + -18798*log(Tr)
end
hco2 = hco2p-hco2r
hh2o = hh2op-hh2or
hn2 = hn2p-hn2r
ho2 =ho2p-ho2r
h2oa=30.5
h2ob=.0103
co2a=28.8
co2b=.028
o2a=27
o2b=.0079
fuela= 38.1
fuelb=.656
Delt_HRP = (MH/2)*(h2oa*Tr +.5*h2ob*Tr.^2) - (h2oa*To +.5*h2ob*To.^2)+MC*(co2a*Tr +.5*co2b*Tr.^2) - (co2a*To +.5*co2b*To.^2)- Ycc*(o2a*Tr +.5*o2b*Tr.^2) - (o2a*To +.5*o2b*To.^2)-(fuela*Tr +.5*fuelb*Tr.^2) - (fuela*To +.5*fuelb*To.^2)
HRP =HRP_to + Delt_HRP
y = (-HRP - MC*hco2 - (MH/2).*hh2o + Ycc.*ho2)./(3.76.*hn2+ho2)
fideal=197.7./(4.76*y*28.97)
f=fideal*efb
%Turbines | HIgh and Low have equal efficiencies
eft=.90
efm=.99
Cpog = (sigg*Rg)/(sigg-1)
% High Pressure
To5= To4 + (WC_HP./(efm*(1+f)*Cpog))
Po5 = Po4*(1-((1-(To5/To4))/eft))^(sigg/(sigg-1))
% Low Pressure
To6= To5 + (WC_LP./(efm*(1+f)*Cpog))
Po6 = Po5*(1-((1-(To6/To5))/eft))^(sigg/(sigg-1))
%Nozzle
%Choke Test
efn = .95
P_Po6 = (1-((1/efn)*(1-(2/(sigg+1)))))^(sigg/(sigg-1))
R=Pa./Po6
Me=1
Pe=Po6*(P_Po6)
Te=To6*(2/(sigg+1))
roe= Pe./(Rg*Te)
Ve= Me*sqrt(sigg*Rg*1000*Te)
% Fs and T.S.F.C Calculations
Fs = ((1+f).*Ve-Va) + (Pe-Pa).*((1+f)./(roe.*Ve))
tsfc{i}= (f./Fs).*3600;
end
celldisp(tsfc) %each subcell represents each value of OPR
tsfc = vertcat(tsfc{:}) %double array
