Loop to create matrix with 2 inputs and 2 outputs
Show older comments
Hello,
I recently had an assignment to create a code and got it right, but afterwards it was recommended to me that it would have been easier to do it a different way. The way the code below is set up was to be more like a calculator. You input the 2 values and it displays 2 outputs. This was very tedious to do and a loop would have been much quicker.
For my situation I have 2 inputs and 2 outputs that are a function of the 2 inputs. I want to have the inputs change through a loop and save the output variables to a matrix.
The value of TINT needs to be set to range from 1100 to 1800 at intervals of 100 while the OPR needs to be 4 to 30 increments of 2. I have a basic understanding of how to do this if it was only one variable changing. I am not sure how to do 2 variables.
%------------------------------------------------------
%Given Values |
Ra = 0.287; %kJ/kg/k |
ga = 1.4; %Gamma Air |
gg = 1.3333; %Gamma gas |
Nd = 0.93; %Efficiency of Diffuser |
Nc = 0.87; %Efficiency of Compressor |
Nb = 0.98; %Efficiency of Combustion Chamber |
Nm = 0.99; %Efficiency of Mechanical |
Nt = 0.90; %Efficiency of Turbines |
Nn = 0.95; %Efficiency of Nozzle |
Ma = 0.84; %Mach |
Alt = 5; %km |
%OPR = 9; %Overal Pressure Ratio |
prompt = 'Overall Pressure Ratio? '; % |
OPR = input(prompt); % |
PRClp = sqrt(OPR); %Pressure Ratio over LP Compressor |
PRChp = sqrt(OPR); %Pressure Ratio over HP Compressor |
deltap0 = 0.04; % 4% of delievery a loss |
%-----------------------------------------------------|
%From Handout 5 based on Altitude of 5km--------------|
Ta = 255.7; %K |
Pa = 54.05; %kPa |
rho = 0.736; %kg/m^3 |
u = 1.63; %N*s/m^2 |
%-----------------------------------------------------|
Va = Ma * sqrt(ga*Ra*Ta*1000); %in m/s
Cpoa = (ga*Ra)/(ga-1); %in kJ/kg/K
T01 = Ta*(1 + ((ga-1)/2)*Ma^2);
P01 = Pa*(1 + Nd*((ga-1)/2)*(Ma^2))^(ga/(ga-1));
P02 = P01*PRClp;
T02 = T01 + ((T01*((PRClp)^((ga-1)/ga)-1)/Nc));
Wclp = Cpoa*(T01-T02);
P03 = P02*PRChp;
T03 = T02 + ((T02*((PRChp)^((ga-1)/ga)-1)/Nc));
Wchp = -Cpoa*(T03-T02);
Tr = T03;
% Tp
prompt = 'Temperature of Products? ';
Tp = input(prompt);
T04 = Tp;
if (Tp <= 1600 && Tp>=400)
COAp = 299180; COBp = 37.85; COCp = -4571.9; CO2Ap = 56835;
CO2Bp = 66.27; CO2Cp = -11634; H2OAp = 88923; H2OBp = 49.36;
H2OCp = -7940.8; N2Ap = 31317; N2Bp = 37.46; N2Cp = -4559.3;
O2Ap = 43388; O2Bp = 42.27; O2Cp = -6635.4;
%Tp > 1600K
elseif (Tp > 1600);
COAp = 309070; COBp = 39.29; COCp = -6201; CO2Ap = 93048;
CO2Bp = 68.58; CO2Cp = -16979; H2OAp = 154670; H2OBp = 60.43;
H2OCp = -19212; N2Ap = 44639; N2Bp = 39.32; N2Cp = -6753.4;
O2Ap = 127010; O2Bp = 46.25; O2Cp = -18798;
end
MC = 14.4;
MH = 24.9;
MO = 0;
Hrp = -8561991.6;
Ycc = MC + (MH/4) + (MO/2);
Ymin = Ycc - (MC/2);
Mfuel = MC*12 + MH + MO*16;
Mair = 28.97;
%Find f
%based on Tr < 1600K
COAr = 299180; COBr = 37.85; COCr = -4571.9; CO2Ar = 56835;
CO2Br = 66.27; CO2Cr = -11634; H2OAr = 88923; H2OBr = 49.36;
H2OCr = -7940.8; N2Ar = 31317; N2Br = 37.46; N2Cr = -4559.3;
O2Ar = 43388; O2Br = 42.27; O2Cr = -6635.4;
%based on Excess Air
N1 = 0;
N2 = MC;
N3 = (MH/2);
HbarCO = 0;
HbarCO2 = (CO2Ap + CO2Bp*Tp + CO2Cp*log(Tp)) - (CO2Ar + CO2Br*Tr + CO2Cr*log(Tr));
HbarH2O = (H2OAp + H2OBp*Tp + H2OCp*log(Tp)) - (H2OAr + H2OBr*Tr + H2OCr*log(Tr));
HbarO2 = (O2Ap + O2Bp*Tp + O2Cp*log(Tp)) - (O2Ar + O2Br*Tr + O2Cr*log(Tr));
HbarN2 = (N2Ap + N2Bp*Tp + N2Cp*log(Tp)) - (N2Ar + N2Br*Tr + N2Cr*log(Tr));
syms Y
eqn = Hrp + (N1 * 282800) + (N2 * HbarCO2) + (N3 * HbarH2O) + ((Y - Ycc)* HbarO2) + ((3.76 * Y)* HbarN2)==0;
Y = solve(eqn,Y);
%Finding f
fideal= ((1/(4.76*Y)))*((Mfuel)/(Mair));
f = fideal/0.98; %Getting f actual from ideal and efficiency
P04 = (1- deltap0)*P03;
Cp0g = (gg*Ra)/(gg-1); %in kJ/kg/K
T05 = T04 +((Wchp)/(Nm*(1+0.0168170)*Cp0g)); %K
P05 = (P04)*(1-((1-(T05/T04))/Nt))^(gg/(gg-1)); %kPa
T06 = T05 + (Wclp/(Nm*(1+f)*Cp0g));
%double(T06);
P06 = P05*(1-((1-(T06/T05))/Nt))^(gg/(gg-1)); %Page 26!!!!
double(P06);
%Return to the particular converging nozzle of the turbojet problem
P0i = P06;
T0i = T06;
T0e = T06;
PstarovP06 = (1-(1/Nn)*(1-(2/(gg+1))))^(gg/(gg-1));
PaovP06 = Pa/P06;
%Choke Test
if(PstarovP06 >= PaovP06)
choked = 1;
elseif (PstarovP06 > PaovP06)
choked =2;
end
%If Choked
if(choked == 1)
Me = 1;
Pe = P06*(PstarovP06);
Te = T06*(2/(gg+1));
rhoe = Pe/(Ra*Te); %kg/m^3
Ve = sqrt(gg*Ra*1000*Te);
double(Ve);
%If not choked
elseif (choked == 2)
Pe = Pa;
Te = T06*(1-(Nn(1-(Pe/P06)^((gg-1)/gg))));
rhoe = Pe/(Ra*Te); %kg/m^3
Me = sqrt(((T06/Te)-1)*(2/(gg-1))); %Me<1 if not choked
Ve = Me*sqrt(gg*Ra*1000*Te);
end
%Compute Fs (Specific Thrust)
Fs = ((1+f)*Ve)-Va +(Pe*1000-Pa*1000)*((1+f)/(rhoe*Ve)); %N/(kg/sec)
double(Fs)
%Compute tsfc
tsfc = f/Fs*3600;
double(tsfc) %kgfuel/hour/N
Below is the way that I tried to edit the code to accomplish this and it has been running for about 5 minutes with no results.
clc clear format long g
%------------------------------------------------------
%Given Values |
Ra = 0.287; %kJ/kg/k |
ga = 1.4; %Gamma Air |
gg = 1.3333; %Gamma gas |
Nd = 0.93; %Efficiency of Diffuser |
Nc = 0.87; %Efficiency of Compressor |
Nb = 0.98; %Efficiency of Combustion Chamber |
Nm = 0.99; %Efficiency of Mechanical |
Nt = 0.90; %Efficiency of Turbines |
Nn = 0.95; %Efficiency of Nozzle |
Ma = 0.84; %Mach |
Alt = 5; %km |
%OPR = 9; %Overal Pressure Ratio |
%prompt = 'Overall Pressure Ratio? '; % |
%OPR = input(prompt); % |
%PRClp = sqrt(OPR); %Pressure Ratio over LP Compressor |
%PRChp = sqrt(OPR); %Pressure Ratio over HP Compressor |
deltap0 = 0.04; % 4% of delievery a loss |
%-----------------------------------------------------|
%From Handout 5 based on Altitude of 5km--------------|
Ta = 255.7; %K |
Pa = 54.05; %kPa |
rho = 0.736; %kg/m^3 |
u = 1.63; %N*s/m^2 |
%-----------------------------------------------------|
% Tp
%prompt = 'Temperature of Products? ';
%Tp = input(prompt);
%T04 = Tp;
for TINT = 1000:1700
TINT = TINT + 100;
for OPR = 2:28
OPR = OPR + 2;
PRClp = sqrt(OPR); %Pressure Ratio over LP Compressor
PRChp = sqrt(OPR); %Pressure Ratio over HP Compressor
Va = Ma * sqrt(ga*Ra*Ta*1000); %in m/s
Cpoa = (ga*Ra)/(ga-1); %in kJ/kg/K
T01 = Ta*(1 + ((ga-1)/2)*Ma^2);
P01 = Pa*(1 + Nd*((ga-1)/2)*(Ma^2))^(ga/(ga-1));
P02 = P01*PRClp;
T02 = T01 + ((T01*((PRClp)^((ga-1)/ga)-1)/Nc));
Wclp = Cpoa*(T01-T02);
P03 = P02*PRChp;
T03 = T02 + ((T02*((PRChp)^((ga-1)/ga)-1)/Nc));
Wchp = -Cpoa*(T03-T02);
Tr = T03;
Tp = TINT;
T04 = Tp;
if (Tp <= 1600 && Tp>=400)
COAp = 299180; COBp = 37.85; COCp = -4571.9; CO2Ap = 56835;
CO2Bp = 66.27; CO2Cp = -11634; H2OAp = 88923; H2OBp = 49.36;
H2OCp = -7940.8; N2Ap = 31317; N2Bp = 37.46; N2Cp = -4559.3;
O2Ap = 43388; O2Bp = 42.27; O2Cp = -6635.4;
%Tp > 1600K
elseif (Tp > 1600);
COAp = 309070; COBp = 39.29; COCp = -6201; CO2Ap = 93048;
CO2Bp = 68.58; CO2Cp = -16979; H2OAp = 154670; H2OBp = 60.43;
H2OCp = -19212; N2Ap = 44639; N2Bp = 39.32; N2Cp = -6753.4;
O2Ap = 127010; O2Bp = 46.25; O2Cp = -18798;
end
MC = 14.4;
MH = 24.9;
MO = 0;
Hrp = -8561991.6;
Ycc = MC + (MH/4) + (MO/2);
Ymin = Ycc - (MC/2);
Mfuel = MC*12 + MH + MO*16;
Mair = 28.97;
%Find f
%based on Tr < 1600K
COAr = 299180; COBr = 37.85; COCr = -4571.9; CO2Ar = 56835;
CO2Br = 66.27; CO2Cr = -11634; H2OAr = 88923; H2OBr = 49.36;
H2OCr = -7940.8; N2Ar = 31317; N2Br = 37.46; N2Cr = -4559.3;
O2Ar = 43388; O2Br = 42.27; O2Cr = -6635.4;
%based on Excess Air
N1 = 0;
N2 = MC;
N3 = (MH/2);
HbarCO = 0;
HbarCO2 = (CO2Ap + CO2Bp*Tp + CO2Cp*log(Tp)) - (CO2Ar + CO2Br*Tr + CO2Cr*log(Tr));
HbarH2O = (H2OAp + H2OBp*Tp + H2OCp*log(Tp)) - (H2OAr + H2OBr*Tr + H2OCr*log(Tr));
HbarO2 = (O2Ap + O2Bp*Tp + O2Cp*log(Tp)) - (O2Ar + O2Br*Tr + O2Cr*log(Tr));
HbarN2 = (N2Ap + N2Bp*Tp + N2Cp*log(Tp)) - (N2Ar + N2Br*Tr + N2Cr*log(Tr));
syms Y
eqn = Hrp + (N1 * 282800) + (N2 * HbarCO2) + (N3 * HbarH2O) + ((Y - Ycc)* HbarO2) + ((3.76 * Y)* HbarN2)==0;
Y = solve(eqn,Y);
%Finding f
fideal= ((1/(4.76*Y)))*((Mfuel)/(Mair));
f = fideal/0.98; %Getting f actual from ideal and efficiency
P04 = (1- deltap0)*P03;
Cp0g = (gg*Ra)/(gg-1); %in kJ/kg/K
T05 = T04 +((Wchp)/(Nm*(1+0.0168170)*Cp0g)); %K
P05 = (P04)*(1-((1-(T05/T04))/Nt))^(gg/(gg-1)); %kPa
T06 = T05 + (Wclp/(Nm*(1+f)*Cp0g));
%double(T06);
P06 = P05*(1-((1-(T06/T05))/Nt))^(gg/(gg-1)); %Page 26!!!!
double(P06);
%Return to the particular converging nozzle of the turbojet problem
P0i = P06;
T0i = T06;
T0e = T06;
PstarovP06 = (1-(1/Nn)*(1-(2/(gg+1))))^(gg/(gg-1));
PaovP06 = Pa/P06;
%Choke Test
if(PstarovP06 >= PaovP06)
choked = 1;
elseif (PstarovP06 > PaovP06)
choked =2;
end
%If Choked
if(choked == 1)
Me = 1;
Pe = P06*(PstarovP06);
Te = T06*(2/(gg+1));
rhoe = Pe/(Ra*Te); %kg/m^3
Ve = sqrt(gg*Ra*1000*Te);
double(Ve);
%If not choked
elseif (choked == 2)
Pe = Pa;
Te = T06*(1-(Nn(1-(Pe/P06)^((gg-1)/gg))));
rhoe = Pe/(Ra*Te); %kg/m^3
Me = sqrt(((T06/Te)-1)*(2/(gg-1))); %Me<1 if not choked
Ve = Me*sqrt(gg*Ra*1000*Te);
end
%Compute Fs (Specific Thrust)
Fs = ((1+f)*Ve)-Va +(Pe*1000-Pa*1000)*((1+f)/(rhoe*Ve)); %N/(kg/sec)
double(Fs);
%Compute tsfc
tsfc = f/Fs*3600;
double(tsfc); %kgfuel/hour/N
Fs(TINT) = Fs;
end
end
Fs
Thank you for any help you may be able to provide,
Colby
Accepted Answer
Geoff Hayes
on 7 Nov 2015
Colby - you almost have it. Your code for the two for loops is as follows
for TINT = 1000:1700
TINT = TINT + 100;
for OPR = 2:28
OPR = OPR + 2;
Since TINT and OPR are your indexing variables, you don't need to increment either through code. This will happen automatically so long as you set your step value correctly in the for loop. The above can be rewritten as
for TINT = 1100:100:1800
for OPR = 4:2:30
% etc.
end
end
So on the first iteration of the outer for loop, TINT will be assigned a value on 1100, and subsequent iterations of this loop will assign values of 1200, 1300, 1400,...,1800 to TINT.
The inner for loop will initialize *OPR to 4 with subsequent iterations assigning it values of 6,8,10,...,30.
3 Comments
Geoff Hayes
on 7 Nov 2015
You probably don't want to use
Fs(TINT)
since TINT is at least 1100 and so your output matrix Fs will be larger than it should be. Also, by doing the above, you have created a matrix named Fs which is the same name as your local variable (scalar) Fs. So on each iteration of your inner for loop, you will be overwriting the matrix with the scalar. (In the end, once all iterating has completed, you would have a matrix but it would only have one non-zero value which would be the last Fs).
You will want to create an mxn array where m is the number of TINT values and n is the number of OPR values (so an 8x14 matrix) as
FsArray = zeros(8,14);
TsfcArray = zeros(8,14);
On each iteration of your inner for loop, you will need to update the appropriate element in each array. Something like
FsArray(u,v) = Fs;
TsfcArray(u,v) = tsfc;
The trick then is to keep track of u and v as they range from 1 to 8 and 1 to 14 respectively.
colbyjh
on 8 Nov 2015
More Answers (0)
Categories
Find more on Number Theory in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
