from __future__ import division
import math
import sympy
from sympy import solve,symbols
#Initializing the variables
g = 9.81;
rho = 1000;
rhoHg = 13.6*rho;
d1 = 0.075;
d2 = 0.025;
Pi = 0.250;
Pt = 0.150;
P_Hg = 0.760;
rho1 = 1.6;
gma = 1.4;
#Calculations
P1 = (Pi+P_Hg)*rhoHg*g;
P2 = (Pt+P_Hg)*rhoHg*g;
rho2 = rho1*(P2/P1)**(1/gma);
V0=symbols('V0')
V1=symbols('V1')
Velo = solve([d2**2*V1*rho2-d1**2*V0*rho1,0.5*(V1**2 - V0**2)*((gma-1)/gma)*(rho2*rho1/(rho2*P1-rho1*P2))-1],[V0,V1])
s=(Velo[1])[0]
Flow = math.pi*d1**2/4*s;
print "Volume of flow (m3/s):",round(Flow,3)
from __future__ import division
import math
#Initializing the variables
Ma = 4;
gma = 1.4;
At = 500; # in mm
#Calculations
N = 1 + (gma-1)*Ma**2/2;
D = (gma+1)/2 ;
#ratio of A/At ==R
R = round( (N/D)**((gma+1)/(2*(gma-1)))/Ma,2);
A=At*R
print "Area of test section (mm^2):",A
from __future__ import division
import math
#Initializing the variables
Ma1 = 2;
gma = 1.4;
T1 = 15+273; # In kelvin
P1 = 105;
#Calculations
Ma2 = (((gma-1)*Ma1**2 +2)/(2*gma*Ma1**2-gma+1))**0.5;
P2 = P1*(1+gma*Ma1**2)/(1+gma*Ma2**2);
T2 = T1*(1 +(gma-1)/2*Ma1**2)/(1 +(gma-1)/2*Ma2**2);
print "Mach No downstream shock wave :",round(Ma2,3)
print "Pressure (bar) of downstream shock wave :",round(P2)
print "Temperature (Degree C) of downstream shock wave :",T2 - 273
from __future__ import division
import math
#Initializing the variables
gma = 1.4;
f = 0.00375;
d = 0.05;
#Calculations
m = d/4;
def x(Ma):
A =(1 -Ma**2 )/(gma*Ma**2);
B = (gma+1)*Ma**2/(2+(gma-1)*Ma**2);
y = m/f*(A+ (gma+1)*math.log(B)/(2*gma));
return y
X1 = x(0.2); # At entrance Ma = 0.2;
X06_X1 =x(0.6); # Section(b) Ma = 0.6;
X06 = X1-X06_X1;
print "The Distance X1 at which the Mach number is unity (m) :",round(X1,2)
print "Distance from the entrance (m) :",round(X06,2)
from __future__ import division
import math
from scipy.optimize import fsolve
#Initializing the variables
gma = 1.4;
Q = 28/60; # m3/s
d = 0.1;
p1 = 200*10**3;
f = 0.004;
x_x1 = 60;
R = 287;
T = 15+273;
#Calculations
A = math.pi*d**2/4;
m = d/4;
v1 = Q/A;
pa = p1*(1-f*(x_x1)*v1**2/(m*R*T))**0.5;
def g(p):
A = (v1*p1)**2/(R*T)
B = f*A*x_x1/(2*m);
y = (p**2 - p1**2)/2 -A*math.log(p/p1) +B;
return y
pb=fsolve(g,pa) # Guessing solution around pa
pb=pb[0]
print "Pressure at the outlet, neglecting velocity changes (kN/m2) :",round(pa/1000,1)
print "Pressure at the outlet, allowing for velocity changes (kN/m2) :",round(pb/1000,1)