blob: 4ada5e0ed393e069cccf87066ec2131581aa7484 (
plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
|
T0 = 37+273; P = 40; g = 1.4;
function [x] = speed(a,b,f)
N = 100;
eps = 1e-5;
if((f(a)*f(b))>0) then
error('no root possible f(a)*f(b)>0');
abort;
end;
if(abs(f(a))<eps) then
error('solution at a');
abort;
end
if(abs(f(b))<eps) then
error('solution at b');
abort;
end
while(N>0)
c = (a+b)/2
if(abs(f(c))<eps) then
x = c ;
x;
return;
end;
if((f(a)*f(c))<0 ) then
b = c ;
else
a = c ;
end
N = N-1;
end
error('no convergence');
abort;
endfunction
deff('[y]=p(x)',['y = x^4 + (5*(x^2)) - 3.225 '])
x = speed(0.5,1,p);
M = x; // Mach number
g = 1.4; // gamma
R = 0.287;
T = T0/(1+((g-1)/2)*M^2);
c = sqrt(g*R*T*1000);
V = c*M;
P0 = P*((T0/T)^(g/(g-1)));
disp(M,"Mach number is")
disp("m/s",V,"Velocity is")
disp("kPa",P0,"Stagnation pressure is")
|
