summaryrefslogtreecommitdiff
path: root/3685/CH17/EX17.1/Ex17_1.sce
blob: 8b7f4636d830a359d54b44a5df93a9b1f6f54897 (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
47
48
49
50
51
52
53
clc
t0 = 37 // Stangnation temperature in degree Celsius
P = 40 // Duct static pressure in kPa
g = 1.4 // Heat capacity ratio
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);
T0 = t0+273;
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)));

printf("\n Example 17.1\n")
printf("\n Mach number is %f ",M)
printf("\n Velocity is %f m/s",V)
printf("\n Stagnation pressure is %f kPa",P0)
//The answers vary due to round off error