-
Notifications
You must be signed in to change notification settings - Fork 1
/
combustionchambera.m
63 lines (56 loc) · 1.46 KB
/
combustionchambera.m
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
54
55
56
57
58
59
60
61
62
63
function [M6,P6,T6,f,pstar,Tstar] = combustionchambera(M5,T5,P5,qhv,Tcomb)
g = 1.4;
R = 287;
Cp = 1004.5;
[T05,~] = isentropic(M5,T5,P5);
Tstar = T5./M5.^2.*((1 + g.*M5.^2)./(g + 1)).^2;
display(Tstar);
pstar = P5.*(1 + g.*M5.^2)./(1 + g);
if Tcomb > Tstar
% M6 = 1;
% P6 = pstar.*(1+g)./(1+g.*M6.^2);
% T6 = Tstar;
syms M6;
T6 = 0.95.*Tstar;
eqn1 = T6./Tstar == M6.^2.*((1 + g)./(1 + g.*M6.^2)).^2;
soln = solve(eqn1,M6);
a = size(soln);
a = max(a);
M6s = [];
for k = 1:a
M6s(k) = double(soln(k,1));
t = isreal(M6s(k));
if M6s(k) > 0 && t==1
if M5 < 1 && M6s(k) < 1
M6 = M6s(k);
% elseif M5 >1 && M6s(k) >1
% M6 = M6s(k);
end
end
end
P6 = pstar.*(1+g)./(1+g.*M6.^2);
elseif Tcomb < Tstar
syms M6;
eqn1 = Tcomb./Tstar == M6.^2.*((1 + g)./(1 + g.*M6.^2)).^2;
soln = solve(eqn1,M6);
a = size(soln);
a = max(a);
M6s = [];
for k = 1:a
M6s(k) = double(soln(k,1));
t = isreal(M6s(k));
if M6s(k) > 0 && t==1
if M5 < 1 && M6s(k) < 1
M6 = M6s(k);
elseif M5 >1 && M6s(k) >1
M6 = M6s(k);
end
end
end
P6 = pstar.*(1+g)./(1+g.*M6.^2);
T6 = Tcomb;
end
[T06,p06] = isentropic(M6,T6,P6);
q = Cp.*(T06 - T05);
f = q./qhv;
end