Solve implicit equation for isentropic flow

17 次查看（过去 30 天）
MATTIA FIORETTO
MATTIA FIORETTO 2022-9-13
编辑： Torsten 2022-9-13
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
Star Strider 2022-9-13
编辑：Star Strider 2022-9-13
Ran in:
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’.
.

更多回答（2 个）

Matt J
Matt J 2022-9-13
编辑：Matt J 2022-9-13
Ran in:
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 个评论
Matt J
Matt J 2022-9-13
编辑：Matt J 2022-9-13
Ran in:
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
Torsten 2022-9-13
编辑：Torsten 2022-9-13
Ran in:
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

类别

在 Help Center 和 File Exchange 中查找有关 Linear Algebra 的更多信息

标签

产品


版本

R2022a

Community Treasure Hunt

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

Start Hunting!

Translated by