How to run script iteratively

6 views (last 30 days)
Vickson Leji
Vickson Leji on 7 Nov 2018
Commented: madhan ravi on 7 Nov 2018
So I just finished a code for a class of mine that works fine. However, I need to be able to change the value of just one variable.
I know I could make that variable equal 2:3:200, but I already have a variable designed in such a way. And if I were to define the variable that I want to change in such a way it wouldn't work because they contain a different number of values. Any Ideas? This is my script:
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=4
Va=Ma*sqrt(siga*Ra*Ta)
Cpoa= (siga*Ra)./(siga - 1)
TINT = 1100:100:1800
PRC = sqrt(OPR)
% 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= (f./Fs).*3600
OPR is the variable I'm interested in

Accepted Answer

madhan ravi
madhan ravi on 7 Nov 2018
Edited: madhan ravi on 7 Nov 2018
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
  5 Comments
Vickson Leji
Vickson Leji on 7 Nov 2018
Awesome, exactly what I needed.
madhan ravi
madhan ravi on 7 Nov 2018
Anytime :) you can also vote for the answer

Sign in to comment.

More Answers (0)

Categories

Find more on MATLAB in Help Center and File Exchange

Tags

Products

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!