MATLAB: Trying to plot equation that changes at a certain point on x-axis

engineeringMATLABplots

In my code below, there is an equation for Tpmax that uses Mc^2. When Mc is less than M0, this is correct, but I need to plot the same equation using M0^2 on the same curve when Mc is greater than M0. The two equations can be seen below as well as an image of how the graph for fuel to air ratio should look. I need this setup to be possible on the graphs for fuel-to-air ratio and TSFC. Thank you

clear
clc
%Constants
Po = 5529.31; %N/m^2
T0 = 217;   %K
gam = 1.4;  
Cp = 1.004; %kJ/kg*K
T4 = 1600;
HPR = [43190 43500 119600];
M0 = 0:0.01:14;
Mc = 4; 
%Figures
f1 = figure();
ax1 = axes();
hold(ax1); 
f2 = figure();
ax2 = axes();
hold(ax2); 
f3 = figure();
ax3 = axes();
hold(ax3); 
f4 = figure();
ax4 = axes();
hold(ax4);
f5=figure();
ax5 = axes();
hold(ax5);
f6=figure();
ax6 = axes();
hold(ax6);
f7=figure();
ax7 = axes();
hold(ax7);
%Equations over Mc
for i=1:numel(HPR)
    hpr = HPR(i); 
    % Gas Constant
    R = ((gam-1)/ gam)*Cp;
    %Velocity of Sound
    a0 = (1000*gam*R*T0).^(1/2);
    %Total to Static Temperature Ratio
    Taur = 1 + (((gam-1)/2)*(M0.^2));
    %Burner Exit Total Temp, T'max
    Tpmax = T4*(1+((gam-1)/2)*Mc.^2);
    
    %Burner Exit Enthalpy to Ambient Enthalpy Ratio
    Taugam = Tpmax/T0;
    %Velocity Ratio
    VR = M0.*((Taugam./Taur).^(1/2));
    %Exit Velocity
    v9 = a0*VR;
    %Temperature at Nozzle Exit
    T9 = (((v9/a0).^2)/M0.^2)*T0;
    %Burner Exit Temp to Burner Inlet Temp
    Taub = T9/T0;
    %Temp at station t2
    Tt2 = Tpmax/Taub;
    %Specific Thrust
    ST = (a0.*M0).*((sqrt(Taugam./Taur)-1));
    %Fuel to Air Ratio
    f = ((Cp*T0)/hpr)*(Taugam-Taur);
    %Thrust Specific Fuel Consumption
    TSFC = ((Cp*T0*(Taugam-Taur))./(a0.*M0.*hpr.*(sqrt(Taugam./Taur)-1)))*10^6;
    %Thermal Efficiency
    Efft = (1-(1./Taur))*100;
    
    %Propulsive Efficiency
    Effp = (2./(((Taugam./Taur).^(1/2))+1))*100;
    
    %Overall Efficiency
    Effo = (Efft.*Effp)/100;
    %Total to Static Temperature Ratio for Maximum Specific Thrust
    Taurmax = Taugam.^(1/3);
    %Mach Number of Flight for Maximum Specific Thrust
    M0max = ((2/(gam-1))*(Taugam.^(1/3)-1)).^1/2;
    %Mach Number of Flight for Minimum TSFC
    A = ((-sqrt(Taugam)/2)+sqrt(((-sqrt(Taugam).^2)/4)+((2^3)/27))).^1/3;
    B = ((-sqrt(Taugam)/2)-sqrt(((-sqrt(Taugam).^2)/4)+((2^3)/27))).^1/3;
    Z = A + B;
    M0min = sqrt((2/(gam-1))*(Z.^2-1));
    %Temperature at Station Tt0
    Tt0 = Taur*T0;
    
    %Thrust Flux Calculations
        %Mass flux
            g = sqrt((gam/R)*(2/(gam+1)^((gam+1)/gam-1)));
            Astar = ((1/Mc.^2)*((2/(gam+1))*(1+((gam-1)/2)*Mc.^2))^((gam+1)/(gam-1))).^(-1/2);
            mflux = g*Astar*(Po/sqrt(T0))*Taugam.^3;
         %Thrust flux
            Tflux = (ST.*mflux)/1000000;
        
    
    %Plots
    plot(ax1, M0,Efft)
    idx2 = Effp < 90;
    plot(ax2, M0(idx2),Effp(idx2))
    idx3 = Effo < 90;
    plot(ax3, M0(idx3),Effo(idx3))
    plot(ax4, M0,TSFC)
    idx5 = M0 > 4;
    plot(ax5, M0,f)
    plot(ax6, M0,ST)
    plot(ax7, M0,Tflux)
end
%Plot Configurations
title(ax1, 'M_{o} vs Thermal Efficiency');
xlabel(ax1,'M_{o}')
ylabel(ax1,'\eta_{t} (%)')
xlim(ax1,[0,14])
ylim(ax1,[0,100])
title(ax2, 'M_{o} vs Propulsive Efficiency');
xlabel(ax2,'M_{o}')
ylabel(ax2,'\eta_{p} (%)')
xlim(ax2,[0,14])
ylim(ax2,[0,100])
legend(ax2,'JP-8: h_{pr}=43,190 kJ/kg','JP-7: h_{pr}=43,500 kJ/kg','Liquid Hydrogen: h_{pr}=119,600 kJ/kg','Location','SouthEast')
title(ax3, 'M_{o} vs Overall Efficiency');
xlabel(ax3,'M_{o}')
ylabel(ax3,'\eta_{o} (%)')
xlim(ax3,[0,14])
ylim(ax3,[0,100])
legend(ax3,'JP-8: h_{pr}=43,190 kJ/kg','JP-7: h_{pr}=43,500 kJ/kg','Liquid Hydrogen: h_{pr}=119,600 kJ/kg','Location','SouthEast')
title(ax4,'M_{o} vs Thrust Specific Fuel Consumption');
xlabel(ax4,'M_{o}')
ylabel(ax4,'TSFC')
xlim(ax4,[0 7])
ylim(ax4,[0 160])
legend(ax4,'JP-8: h_{pr}=43,190 kJ/kg','JP-7: h_{pr}=43,500 kJ/kg','Liquid Hydrogen: h_{pr}=119,600 kJ/kg','Location','NorthEast')
title(ax5, 'M_{o} vs Fuel to Air Ratio');
xlabel(ax5,'M_{o}')
ylabel(ax5,'f')
ylim(ax5,[0 0.2])
legend(ax5,'JP-8: h_{pr}=43,190 kJ/kg','JP-7: h_{pr}=43,500 kJ/kg','Liquid Hydrogen: h_{pr}=119,600 kJ/kg','Location','NorthEast')
title(ax6, 'M_{o} vs Specific Thrust');
xlabel(ax6,'M_{o}')
ylabel(ax6,'Specific Thrust (N/(kg/s)')
xlim(ax6,[0 14])
ylim(ax6,[0 2500])
legend(ax6,'JP-8: h_{pr}=43,190 kJ/kg','JP-7: h_{pr}=43,500 kJ/kg','Liquid Hydrogen: h_{pr}=119,600 kJ/kg','Location','NorthEast')
title(ax7, 'M_{o} vs Thrust Flux');
xlabel(ax7,'M_{o}')
ylabel(ax7,'Thrust Flux (N/(cm^{2})')
xlim(ax7, [0 11])
ylim(ax7, [0 5000])
legend(ax7,'JP-8: h_{pr}=43,190 kJ/kg','JP-7: h_{pr}=43,500 kJ/kg','Liquid Hydrogen: h_{pr}=119,600 kJ/kg','Location','NorthEast')
fprintf('The Mach Number of Flight for Maximum Specific Thrust is %.2f.\n', M0max);
fprintf('The Mach Number of Flight for Minimum TSFC is %.2f.\n', M0min);

Best Answer

Looks like Mc is a constant and M0 is an array. Not sure if this is what you are looking for and how it will affect your other calculations.

%Burner Exit Total Temp, T'max
Tpmax = T4*(1+((gam-1)/2)*M0.^2);
Tpmax(M0>Mc)= T4*(1+((gam-1)/2)*Mc.^2);

Related Solutions

MATLAB: Plotting Curve of Function over 5 parameters

See this code

%Constants
T0 = 217;
gam = 1.4;
Cp = 1.004;
T4 = 1600;
hpr = 42800;
M0 = 0:1:14;
MC = 0:4; % <--- define as vector
%Equations
f = figure();
ax = gca;
hold(ax); % to plot all lines on same axes
for i=1:numel(MC)
    Mc = MC(i); % select one value of Mc
    % Gas Constant
    R = ((gam-1)/ gam)*Cp;
    %Velocity of Sound
    a0 = (1000*gam*R*T0).^(1/2);
    %Total to Static Temperature Ratio
    Taur = 1 + (((gam-1)/2)*(M0.^2));
    %Burner Exit Total Temp, T'max
    Tpmax = T4*(1+((gam-1)/2)*Mc.^2);
    %Burner Exit Enthalpy to Ambient Enthalpy Ratio
    Taugam = Tpmax/T0;
    %Velocity Ratio
    VR = M0.*((Taugam./Taur).^(1/2));
    %Exit Velocity
    v9 = a0*VR;
    %Temperature at Nozzle Exit
    T9 = (((v9/a0).^2)/M0.^2)*T0;
    %Burner Exit Temp to Burner Inlet Temp
    Taub = T9/T0;
    %Temp at station t2
    Tt2 = Tpmax/Taub;
    %Specific Thrust
    ST = (a0.*M0).*((sqrt(Taugam./Taur)-1));
    %Fuel to Air Ratio
    f = ((Cp*Tt2)/hpr)*((Tpmax/Tt2)-1);
    %Thrust Specific Fuel Consumption
    TSFC = (Cp*T0*(Taugam-Taur))./(a0.*M0.*hpr.*(sqrt(Taugam./Taur)-1));
    %Thermal Efficiency
    Efft = (1-(1./Taur))*100;
    %Propulsive Efficiency
    Effp = (2./(((Taugam./Taur).^(1/2))+1))*100;
    %Overall Efficiency
    Effo = (Efft.*Effp)/100;
    %Total to Static Temperature Ratio for Maximum Specific Thrust
    Taurmax = Taugam.^(1/3);
    %Mach Number of Flight for Maximum Specific Thrust
    M0max = ((2/(gam-1))*((T4/T0)*(1+((gam-1)/2)*Mc.^2)).^1/3).^1/2;
    %Temperature at Station Tt0
    Tt0 = Taur*T0;
    plot(M0,Effo)
end

MATLAB: Matlab freezing up when plotting a piecewise function

It seems that fplot has problems evaluating the piecewise function. The only way to get around that is to create your own anonymous function from it, and use that. That was initially straightforward, just lightly editing the piecewise expression:

ST = @(M0) (M0 <= 1.0) .* (- 297.2421235*M0 - 297.2421235*(20.96007134/(0.632149*M0^2 + 0.90307)^(3/13))^(1/2)*(0.02171962353*(0.7*M0^2 + 1.0)^(3/13) - 1.292400727)*(3.2114*(0.632149*M0^2 + 0.90307)^(3/13) - 1.0)^(1/2)) + (1.0 < M0) .* (- 297.2421235*M0 - 297.2421235*(3.2114*(-1.0*(0.7*M0^2 + 1.0)*(0.06773025*(M0 - 1.0)^(27/20) - 0.90307))^(3/13) - 1.0)^(1/2)*(20.96007134/(-1.0*(0.7*M0^2 + 1.0)*(0.06773025*(M0 - 1.0)^(27/20) - 0.90307))^(3/13))^(1/2)*(0.02171962353*(0.7*M0^2 + 1.0)^(3/13) - 1.292400727));

however fplot then threw:

Warning: Function behaves unexpectedly on array inputs.

Using str2func and vectorize to eliminiate that problem:

ST = str2func(['@(M0) ' vectorize('(- 297.2421235*M0 - 297.2421235*(20.96007134/(0.632149*M0^2 + 0.90307)^(3/13))^(1/2)*(0.02171962353*(0.7*M0^2 + 1.0)^(3/13) - 1.292400727)*(3.2114*(0.632149*M0^2 + 0.90307)^(3/13) - 1.0)^(1/2)) + (1.0 < M0) .* (- 297.2421235*M0 - 297.2421235*(3.2114*(-1.0*(0.7*M0^2 + 1.0)*(0.06773025*(M0 - 1.0)^(27/20) - 0.90307))^(3/13) - 1.0)^(1/2)*(20.96007134/(-1.0*(0.7*M0^2 + 1.0)*(0.06773025*(M0 - 1.0)^(27/20) - 0.90307))^(3/13))^(1/2)*(0.02171962353*(0.7*M0^2 + 1.0)^(3/13) - 1.292400727));')]);

the fplot call works:

fplot(ST,[0.1 4.5]) %Specifc thrust vs flight mach
title('specific thrust - freestream mach relationship at 99% of tauc')
xlabel('M0 - freestream mach number')
ylabel('ST - specific thrust (N/[kg/s])')

Best Answer

Related Solutions

MATLAB: Plotting Curve of Function over 5 parameters

MATLAB: Matlab freezing up when plotting a piecewise function

Related Question