Solve implicit equation for isentropic flow

27 views (last 30 days)
I have to resolve the equation of Isentropic flow that links Area ratio and Mach number.
I have solved in this way but the result is different from the true result that you can find online with an isentropic calculator or with tables.
fcn = @(M) (1./M)*((2/(g+1)).*(1+(((g-1)/2)*M.^2)).^((g+1)/(2*(g-1)))) - Aratio;
M = fzero(fcn, 1);
For example Aratio=4, g=1.4, the true result is M=2.94, but with the coding I get 2.557 and this value doesn't depend on the initial value.
I could use also this code
[mach,T,P,rho,area] = flowisentropic(gamma,4,'sup') and the result of Mach number is correct
but I would like to know why the previous code is incorrect

Accepted Answer

Star Strider
Star Strider on 13 Sep 2022
Edited: Star Strider on 13 Sep 2022
When in doubt, plot the function —
Aratio=4;
g=1.4;
% fcn = @(M) (1./M).*((2/(g+1)).*(1+(((g-1)/2).*M.^2)).^((g+1)/(2*(g-1)))) - Aratio;
fcn = @(M) (1./M).*( (2/(g+1)) .* (1+(g-1)/2*M.^2) ) .^ ((g+1)/(2*(g-1))) - Aratio;
Mv = linspace(0, 4);
idx = find(diff(sign(fcn(Mv))))
idx = 1×2
4 73
for k = 1:numel(idx)
[Mc(k),fv(k)] = fzero(fcn, Mv(idx(k)));
end
Mc
Mc = 1×2
0.1465 2.9402
fv
fv = 1×2
1.0e-15 * 0.8882 0.8882
figure
plot(Mv, fcn(Mv), '-b', Mc, zeros(size(Mc)),'sr')
grid
text(Mc,zeros(size(Mc)),compose(' \\leftarrow %.4f',Mc))
There are two roots. There appears to be a slight coding error in ‘fcn’ since when I re-coded it, I get the desired root.
EDIT — (13 Sep 2022 at 11:45)
Recoded ‘fcn’.
.

More Answers (2)

Matt J
Matt J on 13 Sep 2022
Edited: Matt J on 13 Sep 2022
For example Aratio=4, g=1.4, the true result is M=2.94
We can see below that M=2.94 is not a solution, so either you have atypo in fcn, or your expectations are wrong.
Aratio=4; g=1.4;
fcn = @(M) (1./M)*((2/(g+1)).*(1+(((g-1)/2)*M.^2)).^((g+1)/(2*(g-1)))) - Aratio;
fcn(2.94)
ans = 1.7590
  1 Comment
Matt J
Matt J on 13 Sep 2022
Edited: Matt J on 13 Sep 2022
so either you have atypo in fcn, or your expectations are wrong.
Apparently, the former:
fcn=@(M)isentropic(M,4,1.4);
[M,fval]=fzero(fcn,[1,4])
M = 2.9402
fval = 8.8818e-16
function out=isentropic(M, Aratio, gamma)
gp1=gamma+1; gm1=gamma-1;
tmp=1+gm1/2*M^2;
tmp=tmp*2/gp1;
tmp=tmp.^(gp1/2/gm1);
out=tmp/M-Aratio;
end

Sign in to comment.


Torsten
Torsten on 13 Sep 2022
Edited: Torsten on 13 Sep 2022
gamma = 1.4;
[mach,T,P,rho,area] = flowisentropic(gamma,4,'sup')
mach = 2.9402
T = 0.3664
P = 0.0298
rho = 0.0813
area = 4
Aratio = 4.0;
gamma = 1.4;
fcn = @(M) (1/M)*(2/(gamma+1)*(1+(gamma-1)/2*M.^2))^((gamma+1)/(2*(gamma-1))) - Aratio;
sol = fzero(fcn,2)
sol = 2.9402

Categories

Find more on Gas Dynamics in Help Center and File Exchange

Products


Release

R2022a

Community Treasure Hunt

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

Start Hunting!