Solve implicit equation for isentropic flow

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

Answer 1

[Star Strider]

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

.

Answer 2

[Matt J]

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]

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

Answer 3

[Torsten]

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