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])')
采纳的回答
Star Strider
2020-6-6
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])')
.
更多回答（0 个）
类别
在 帮助中心 和 File Exchange 中查找有关 Brakes and Detents 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
