# Solve implicit equation for isentropic flow

32 views (last 30 days)
MATTIA FIORETTO on 13 Sep 2022
Edited: Torsten on 13 Sep 2022
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

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’.
.

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
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

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

R2022a

### Community Treasure Hunt

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

Start Hunting!