Trying to plot specific values from a vector

4 views (last 30 days)
Caroline Kenemer on 3 May 2021
Commented: Mathieu NOE on 4 May 2021
Hi, I am trying to plot the TSFC values and Pcr values for specfic values of To4 which are 1500,1600, and 1700. So that I will end up with three lines. How could I do this?

Mathieu NOE on 3 May 2021
hello
this is my suggestion - I put the equations in a sub function to make it easy after to do a for loop for multiple To4 values
also changes the Pcr from linspace to logspace vector as it seems to better render in the semilogx plot (my opinion)
clc
clearvars
%
%Turbojet - Rolls Royce Olympus 593
%Range = 6440 km
%Engine Weight = 3386 kg;
%Ramjet is based off the Blackbird SR-71
%% Takeoff TJ
Pa = 101325; % in Pa
Pe = Pa;
rho = 1.225; %kg/m^3
Ai = 0.4902; %inlet area in meters
ma = 64.4; %mass flow rate in kg/s
u = ma/(rho*Ai);
h = 0; %km
T = 169000; % Thrust in N
R = 287; %J/kg*K
Ta = 288.15; % ambient temperature at Sea Level
gamma = 1.4;
Pcr = logspace(0,2,25);
To4 = [1500;1600;1700];
for ci = 1:length(To4)
TSFC(ci,:) = mysubfn1(Pcr,u,gamma,R,Ta,Pa,ma,T,To4(ci))
leg_str{ci} = ['To4 = ' num2str(To4(ci))];
end
semilogx(Pcr,TSFC)
legend(leg_str);
function [TSFC] = mysubfn1(Pcr,u,gamma,R,Ta,Pa,ma,T,To4)
M = u/sqrt(gamma*R*Ta);
%compressor inlet conditions
To2 = Ta*(1+((gamma-1)/2)*M^2);
yd = 1.4;
nd = 0.97;
Po2 = Pa*(1+nd*((To2/Ta)-1))^(yd/(yd-1));
% Pcr = linspace(0,100,10);
% To4 = linspace(1500,1700,10);
Po3 = Po2.*Pcr;
yc = 1.37;
nc = 0.85;
To3 = To2.*(1+(1./nc).*(Pcr.^((yc-1)./yc))-1);
Q = 45000000; %heating value in K/kg
Cp = 1000; %J/kgK
f = ((To4./To3)-1)./(((Q.*To3)./Cp)-To4./To3);
Po4 = Po3;
To5 = To4-(To3-To2);
nt = 0.90;
yt =1.33;
Po5 = Po4.*(1-(1/nt).*(1-(To5./To4))).^(yt/(yt-1));
%nozzle conditions with no afterburner
To6 = To5;
Po6 = Po5;
yn = 1.36;
nn = 0.98;
ue = sqrt(2.*nn.*(yn./(yn-1)).*R.*To6.*(1-(Pa./Po6).^(yn-1)./yn));
mf = f.*ma;
TSFC = mf./T;
end
Mathieu NOE on 4 May 2021
hello again
1/ if you have multiple functions in one code , the functions must always be at the very end of the code, all "main" code must be written in the upper section; you can have multiple functions , simply make sure they have individual names
2/ xlabel, ylabel, and title : it's exactly this way to code it :
semilogx(Pcr,TSFC)
title('TSFC vs. Pcr vs. To4')
legend(leg_str);
xlabel('Pcr')
ylabel('TSFC')
3/ sure, you can always get back to linspace
or this way : Pcr = (0:10:100);
did you notice that Pcr = 0 will issue a NaN result because you make a division by zero in this line :
f = ((To4./To3)-1)./(((Q.*To3)./Cp)-To4./To3);
I send back the full code :
clc
clearvars
%
%Turbojet - Rolls Royce Olympus 593
%Range = 6440 km
%Engine Weight = 3386 kg;
%Ramjet is based off the Blackbird SR-71
%% Takeoff TJ
Pa = 101325; % in Pa
Pe = Pa;
rho = 1.225; %kg/m^3
Ai = 0.4902; %inlet area in meters
ma = 64.4; %mass flow rate in kg/s
u = ma/(rho*Ai);
h = 0; %km
T = 169000; % Thrust in N
R = 287; %J/kg*K
Ta = 288.15; % ambient temperature at Sea Level
gamma = 1.4;
% Pcr = logspace(0,2,25);
Pcr = (0:10:100);
To4 = [1500;1600;1700];
for ci = 1:length(To4)
TSFC(ci,:) = mysubfn1(Pcr,u,gamma,R,Ta,Pa,ma,T,To4(ci))
leg_str{ci} = ['To4 = ' num2str(To4(ci))];
end
figure(1),
plot(Pcr,TSFC)
title('TSFC vs. Pcr vs. To4')
legend(leg_str);
xlabel('Pcr')
ylabel('TSFC')
function [TSFC] = mysubfn1(Pcr,u,gamma,R,Ta,Pa,ma,T,To4)
M = u/sqrt(gamma*R*Ta);
%compressor inlet conditions
To2 = Ta*(1+((gamma-1)/2)*M^2);
yd = 1.4;
nd = 0.97;
Po2 = Pa*(1+nd*((To2/Ta)-1))^(yd/(yd-1));
% Pcr = linspace(0,100,10);
% To4 = linspace(1500,1700,10);
Po3 = Po2.*Pcr;
yc = 1.37;
nc = 0.85;
To3 = To2.*(1+(1./nc).*(Pcr.^((yc-1)./yc))-1);
Q = 45000000; %heating value in K/kg
Cp = 1000; %J/kgK
f = ((To4./To3)-1)./(((Q.*To3)./Cp)-To4./To3);
Po4 = Po3;
To5 = To4-(To3-To2);
nt = 0.90;
yt =1.33;
Po5 = Po4.*(1-(1/nt).*(1-(To5./To4))).^(yt/(yt-1));
%nozzle conditions with no afterburner
To6 = To5;
Po6 = Po5;
yn = 1.36;
nn = 0.98;
ue = sqrt(2.*nn.*(yn./(yn-1)).*R.*To6.*(1-(Pa./Po6).^(yn-1)./yn));
mf = f.*ma;
TSFC = mf./T;
end

R2020a

Community Treasure Hunt

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

Start Hunting!