import math
import sys
# Variables
T0 = 37+273
P = 40.
g = 1.4;
def speed(a,b,f):
N = 100.;
eps = 1e-5;
if((f(a)*f(b))>0):
print('no root possible f(a)*f(b)>0');
sys.exit(0)
if(abs(f(a))<eps):
print('solution at a');
sys.exit(0)
if(abs(f(b))<eps):
print('solution at b');
sys.exit(0)
while(N>0):
c = (a+b)/2.
if(abs(f(c))<eps):
x = c ;
return x;
if((f(a)*f(c))<0 ):
b = c ;
else:
a = c ;
N = N-1;
print('no convergence');
sys.exit(0)
def p(x):
return x**4 + (5*(x**2)) - 3.225
x = speed(0.5,1,p);
M = x; # Mach number
g = 1.4; # gamma
R = 0.287;
# Calculation
T = T0/(1+((g-1)/2)*M**2);
c = math.sqrt(g*R*T*1000);
V = c*M;
P0 = P*((T0/T)**(g/(g-1)));
# Results
print "Mach number is",round(M,2)
print "Velocity is",round(V,1),"m/s"
print "Stagnation pressure is",round(P0,2),"kPa"
# note: rounding off error is there.
import math
P1 = 0.18e03; # in Kpa
R = 0.287
T1 = 310.;
P0 = 0.1e03;
A1 = 0.11;
V1 = 267.;
# Calculation
w = (P1/(R*T1))*A1*V1;
g = 1.4;
c1 = math.sqrt(g*R*T1*1000);
M1 = V1/c1;
A1A_ = 1.0570; # A1./A* A* = A_
P1P01 = 0.68207;
T1T01 = 0.89644;
F1F_ = 1.0284;
A2A1 = 0.44/0.11 ; # A2A1 = A2/A1
A2A_ = A2A1*A1A_;
M2 = 0.135;
P2P02 = 0.987;
T2T02 = 0.996;
F2F_ = 3.46;
P2P1 = P2P02/P1P01;
T2T1 = T2T02/T1T01;
F2F1 = F2F_/F1F_;
P2 = P2P1*P1;
T2 = T2T1*T1;
A2 = A2A1*A1;
F1 = P1*A1*(1+g*M1**2);
F2 = F2F1*F1;
Tint = F2-F1;
Text = P0*(A2-A1);
NT = Tint - Text ;
# Results
print "Net thrust is",round(NT,2),"kN"
# rounding off error is there.
import math
M2 = 2.197;
P2P0 = 0.0939;
T2T0 = 0.5089;
P0 = 1000.
T0 = 360.;
g = 1.4;
R = 0.287;
# Calculation and Results
P2 = P2P0*P0;
T2 = T2T0*T0;
c2 = math.sqrt(g*R*T2*1000);
V2 = c2*M2;
# for air
P_P0 = 0.528; T_T0 = 0.833; # T_ == T*
P_ = P_P0*P0; T_ = T_T0*T0;
rho_ = P_/(R*T_);
V_ = math.sqrt(g*R*T_*1000);
At = 500e-06; # throat area
w = At*V_*rho_;
print ("When divergent section act as a nozzle")
print "Maximum flow rate of air is",round(w,3),"kg/s"
print "Static temperature is",round(T2,1),"K"
print "Static Pressure is",P2,"kPa"
print "Velocity at the exit from the nozzle is",round(V2,0),"m/s"
# Part (b)
Mb = 0.308;
P2P0b = 0.936;
T2T0b = 0.9812;
P2b = P2P0b*P0;
T2b = T2T0b*T0;
c2b = math.sqrt(g*R*T2b*1000);
V2b = c2b*Mb;
print ("\nWhen divergent section act as a diffuser")
print "Maximum flow rate of air is",round(w,3),"kg/s"
print "Static temperature is",round(T2b,1),"K"
print "Static Pressure is",P2b,"kPa"
print "Velocity at the exit from the nozzle is",round(V2b,0),"m/s"
# Variables
Px = 16. #kPa;
Poy = 70. #kPa;
Mx = 1.735;
Pyx = 3.34; # Pyx = Py/Px
rho_yx = 2.25;
Tyx = 1.483; Poyox = 0.84; My = 0.631;
Tox = 573; Toy = Tox;
# Calculation
Tx = Tox/(1+((g-1)/2.)*Mx**2);
Ty = Tyx*Tx;
Pox = Poy/Poyox ;
# From table
Mx = 1.735;
# Results
print "Mach number of the tunnel is",Mx
import math
# Variables
Ax = 18.75;
A_ = 12.50; # A_= A*
AA_ = 1.5; # A/A*
Mx = 1.86;
Pxox = 0.159;
R = 0.287;
Pox = 0.21e03; # in kPa
# Calculation
Px = Pxox*Pox;
# from the gas table on normal shock
Mx = 1.86; My = 0.604; Pyx = 3.87; Poyx = 4.95; Poyox = 0.786;
Py = Pyx*Px;
Poy = Poyx*Px;
My = 0.604;
Ay_ = 1.183;
A2 = 25.; Ay = 18.75;
A2_ = (A2/Ay)*Ay_;
# From isentropic table
M2 = 0.402;
P2oy = 0.895;
P2 = P2oy*Poy;
syx = -R*math.log(Poy/Pox); # sy-sx
# Results
print "Exit mach number is M2",M2
print "Exit pressure is",round(P2,2),"kPa"
print "Exit Stagnation pressure is",round((Pox-Poy),1),"kPa"
print "Entropy increase is",round(syx,4),"kJ/kg K"
import math
# Variables
g = 1.4
R = 0.287
d = 1.4; # del
P0 = 1.4; # in bar
T0 = 280.
T1 = T0;
cp = 1.005
A2 = 0.0013
# Calculation
P_ = P0/((g+1)/2.)**(d/(d-1)) ; # P_ = P*
P1 = P0; Pb = 1.; P2 = Pb;
T2 = T1*(P2/P1)**((d-1)/d);
V2 = math.sqrt(2*cp*(T1-T2)*1000);
m_dot = (A2*V2*P2*100)/(R*T2);
# Results
print "Mass flow rate is",round(m_dot,4),"kg/s"
print "The mass flow rate can be increased by raising the supply pressure"
import math
# Variables
Mx = 1.8
Pyx = 3.6133;
Px = 0.5; Tx = 280.; Ty = 429.;
R = 0.287
# Calculation and Results
Py = Pyx*Px; cp = 1.005;
print "Pressure Py is",round(Py,4),"bar"
Pxox = 0.17404;
Pox = Px/Pxox;
print "Stagnation pressure is",round(Pox,2),"bar"
Txox = 0.60680;
Tox = Tx/Txox;
print "Stagnation temperature is",round(Tox,1),"K"
sysx = cp*math.log(Ty/Tx)-R*math.log(Py/Px);
print "The change in specific entropy is",round(sysx,5),"kJ/kg K"
# note : rounding off error.
# variables
M1 = 0.39
T1ox = 0.97032
Tox = 546.2
P1ox = 0.9
pox = 3.778 # bar
Mx = 2
Pxox = 0.1278
Px = 0.4828 # bar
My = 0.57735
Pyoy = 0.79737
Poy = 2.735 # bar
Toy = 546.2 # K
Tox = Toy
# calculations and Results
A1 = 0.5*0.287*530/(180*3.4*100) # m**2
print "A1 = %.4e m**2"%A1
M1 = 0.39
A1 = 1.6346
Astar = 7.602*10**-4
print "A* = %.3e m**2"%(Astar)
M = 0.57735
A2 = 1.25*Astar
print "A2 = Exit plane area = %.3e m**2"%A2