MATLAB Answers

Trying to plot specific values from a vector

4 views (last 30 days)
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?

Accepted Answer

Mathieu NOE
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
  4 Comments
Mathieu NOE
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

Sign in to comment.

More Answers (0)

Tags

Products


Release

R2020a

Community Treasure Hunt

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

Start Hunting!