Matlab freezing up when plotting a piecewise function

Hi there,

I'm running into some difficulties with matlab when I try to plot the function ST (y) against M0 (x) in the range 0.1<M0<4.5. When I try to run the script matlab seems to lock up and after 30 mins of ticking over, still hasn't plotted the function. I tried plotting over a smaller range (4.0<M0<4.5) and had no luck but did manage to get it to plot when I set the piecewise function (Pii) to a regular function. This leads me to think that there is some interaction with the piecewise function that is causing the issues.

It would be great if anyone knew why this is happening or had a solution to the problem. I am using MATLAB R2019b.

Regards,

Matt.

close all

clear all

syms a0 gamt gamc gam0 Rt Rc R0 tauLam tauc tauR Pin Pib Pii Pi0 M0 T0 Cpc Cpt Cpf Cp02 CpR H PhiR PhiB fst f02

tauc = 3.2114

Cpc = 1004; %J/(kg*K)

Cpt = 1239; %J/(kg*K)

Cp02 = 918; %J/(kg*K)

Cpf = 2010; %J/(kg*K)

CpR = 1844; %J/(kg*K)

fst = 0.067; %stoichiometric fuel-air ratio

f02 = 0.22; %Atmospheric mass fraction of oxygen

H = 44e6; %J/kg fuel heating value

gamt = 1.3;

gamc = 1.4;

gam0 = gamc;

Rc = 286.86; %J/(kg*K)

Rt = 285.92; %J/(kg*K)

R0 = Rc; %J/(kg*K)

T0 = 220; %K

a0 = sqrt(gam0*R0*T0); %m/s

tauLam = 2400/T0; %Tt4/T0

tauR = 2000/T0; %TtR0/T0

Pii = piecewise(M0<=1,0.98,M0>1,0.98*(1-0.075*(M0-1)^1.35)) %Stagnation pressure ratio for frictional losses in the inlet and inlet shocks

Pib = 0.95; %Burner stagnation pressure ratio

Pin = 0.97; %Nozzle stagnation pressure ratio

Pi0 = 1+100000/((200000*R0*T0)/(a0^2*M0^2)); %Inlet stagnation pressure ratio

tau0 = Pi0^((gamt-1)/gamt);

PhiB0 = ((Cpc*tau0)/(Cpt*tauLam))-1;

PhiB1 = fst*((CpR*tauR)/Cpf)-1;

PhiB2 = fst+(f02*Cp02)/Cpf+(fst*H)/(Cpf*T0)-(fst+f02)*((CpR*tauR)/Cpf);

PhiB3 = (fst/(Cpt*tauLam))*(Cpf+H/T0);

PhiB4 = PhiB3+(Cp02*f02)/(Cpt*tauLam)-(fst+f02);

PhiB = PhiB0/(fst-PhiB3-(PhiB1/PhiB2)*PhiB4);

PhiR = (PhiB*(fst*CpR*tauR-fst*Cpf))/(fst*Cpf+fst*(H/T0)-fst*CpR*tauR+f02*Cp02-f02*CpR*tauR);

beta0 = Cpc/(Cpt*tauLam);

beta1 = tau0+(PhiR+PhiB)*(fst/Cpc)*(Cpf+(H/T0))+PhiR*f02*(Cp02/Cpc);

beta = beta0*beta1-1

ST0 = (1+beta);

ST1 = tauc*(Pin*Pib*Pii*Pi0)^((gamt-1)/gamt);

ST2 = sqrt(((gamt*Rt)/(gamc*Rc))*(tauLam/ST1)*(2/(gamt-1)));

ST3 = sqrt(ST1-1);

ST = a0*(ST0*ST2*ST3-M0)

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])')

Answer 1

[Star Strider]

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])')

.

[Matthew Browning]

Worked great, thanks for your help.

[Star Strider]

As always, my pleasure!