#what is the gas constant of air and density of air
import math
#intilization variable
p=3*10**6 ; #pressure in Pa
t=298. ; #temperatue in kelvin
mw= 29.; #molecular weight in kg/mol
ru=8314.; #universal constant in J/kmol.K
r=ru/mw ;
#using perfect gas law to get density:
rho=p/(r*t) ;
print'%s %.2f %s'%('Gas constant of air in',r,'J/kg.K')
print'%s %.1f %s'%('Density of air in',rho,'kg/m^3')
#find out the exit temperature and exit density by various methods
import math
t1=288.; #inlet temperture in Kelvin
p1=100*10**3; #inlet pressure in Pa
p2=1*10**6 #exit pressure in Pa
gma=1.4; #gamma.
rg=287.; #gas constant in J/kg.K
t2=t1*(p2/p1)**((gma-1)/gma); #exit temperature
print'%s %.5f %s'%('Exit temperature in',t2,'K')
#first method to find exit density:
#application of perfect gas law at exit
rho=p2/(rg*t2); #rho= exit density.
print'%s %.7f %s'%('exit density at by method 1 in',rho,'kg/m^3')
#method 2: using isentropic relation between inlet and exit density.
rho1=p1/(rg*t1); #inlet density.
rho=rho1*(p2/p1)**(1/gma);
print'%s %.2f %s'%('exit density by method 2 in',rho,'kg/m^3')
#what is the rate of mass flow through exit
import math
d1=1.2 #inlet 1 density in kg/m^3.
u1=25. # inlet 1 veocity in m/s.
a1=0.25 #inlet 1 area in m^2.
d2=0.2 #inlet 2 density in kg/m^3.
u2=225. #inlet 2 velocity in m/s.
a2=0.10 #inlet 2 area in m^2.
m1=d1*a1*u1; #rate of mass flow entering inlet 1.
m2=d2*u2*a2; #rate of mass flow entering inlet 2.
#since total mass in=total mass out,
m3=m1+m2; #m3=rate of mass flow through exit.
print'%s %.f %s'%('Rate of mass flow through exit in',m3,' kg/s')
#what is the axial force needed to support the plate and lateral force needed to support the plate
import math
u1=2 #speed of water going on the plate. X-component in m/s.
v1=0 #speed of water going on the plate. Y-component in m/s.
u2=1 #speed of water going on the plate. X-component in m/s.
v2=1.73 #speed of water going on the plate Y-coponent in m/s.
m=0.1 #rate of flow of mass of the water on the plate in kg/s.
#Using Newton's second law.
Fx=m*(u2-u1); #X-component of force exerted by water
print'%s %.1f %s'%('Axial force needed to support the plate in',Fx,'N')
Fy=m*(v2-v1); #Y-component of force exerted by water.
print'%s %.3f %s'%('Lateral force needed to support the plate in',Fy,'N')
#calculate the Exit total and static temperature
m=50 #mass flow rate in kg/s.
T1=298 #inlet temperature in K.
u1=150 #inlet velocity in m/s.
cp1=1004 #specific heat at constant pressure of inlet in J/kg.K.
gm=1.4 #gamma.
u2=400 # exit velocity in m/s.
cp2=1243. #specific heat at constant pressure of exit in J/kg.K.
q=42*10**6 #heat transfer rate in control volume in Watt.
me=-100*10**3 #mechanical power in Watt.
#first calculate total enthalpy at the inlet:
ht1=cp1*T1+(u1**2)/2; #ht1=Total inlet enthalpy.
#now applying conservation of energy equation:
ht2=ht1+((q-me)/m) #ht2=Total enthalpy at exit.
Tt2=ht2/cp2; #Tt2=Total exit temperature.
T2=Tt2-((u2**2)/(2*cp2)); #T2=static exit temperature.
print'%s %.5f %s'%('Exit total temperature in',Tt2,'K')
print'%s %.4f %s'%('Exit static temperature in',T2,'K')
#intilization variable
import math
d=0.2 #Diameter in meters.
M1=0.2 #inlet Mach no.
p1=100*10**3 #inlet pressure in Pa
Tt1=288. #total inlet temperature in K
q=100*10**3 #rate of heat transfer to fluid in Watt.
rg=287. #Gas constant in J/kg.K.
gm=1.4 #gamma
#(a)inlet mass flow:
m=((gm/rg)**(1./2.))*(p1/(Tt1)**(1./2.))*3.14*(d*d)/4.*(M1/(1.+((gm-1.)/2.)*(M1**2.))**((gm+1.)/(2.*(gm-1.))));
#(b)
qm=q/m; #Heat per unit mass.
#Tt1/Tcr=0.1736, pt1/Pcr=1.2346, ((Delta(s)/R)1=6.3402,p1/Pcr=2.2727)
Tcr=Tt1/0.1736;
Pcr=p1/2.2727;
#From energy equation:
cp=(gm/(gm-1.))*rg;
Tt2=Tt1+(q/cp);
q1cr=cp*(Tcr-Tt1)/1000.;
M2=0.22;
#From table : pt2/Pcr=1.2281, (Delta(s)/R)2=5.7395, p2/Pcr=2.2477.
#The percent total pressure drop is (((pt1/Pcr)-(pt2/Pcr))/(pt1/Pcr))*100.
p2=2.2477*Pcr;
dp=((1.2346-1.2281)/1.2346)*100;
#Entropy rise is the difference between (delta(s)/R)1 and (delta(s)/R)2.
ds=6.3402-5.7395;
#Static pressure drop in duct due to heat transfer is
dps=((p1/Pcr)-(p2/Pcr))*Pcr/1000.;
print'%s %.7f %s'%('Mass flow rate through duct in',m,'kg/s')
print'%s %.4f %s'%('Critical heat flux that would choke the duct for the M1 in',q1cr,'kJ/kg')
print'%s %.2f %s'%('The exit Mach No.',M2,'')
print'%s %.7f %s'%('The percent total pressure loss',dp,'%')
print'%s %.4f %s'%('The entropy rise',ds,'')
print'%s %.7f %s'%('The static pressure drop in ',dps,'kPa')
#what is total exit temperautre if exit is choked and maximum heat released and fule to air ratio to thermally choke the combustor exit and total pressure loss
#intilization variable
import math
M1=3.0 ##Mach no. at inlet
pt1=45*10**3 ##Total pressure t inlet in Pa
Tt1=1800 ##Total temperature at inlet in K
hv=12000 ##Lower heating value of hydrogen kJ/kg
gm=1.3 ##gamma
R=0.287 ##in kJ/kg.K
##Using RAYLEIGH table for M1=3.0 and gamma=1.3, we get Tt1/Tcr=0.6032, pt1/Pcr=4.0073.
Tcr=Tt1/0.6032
Pcr=pt1/4.0073
##if exit is choked, Tt2=Tcr
Tt2=Tt1/0.6032;
cp=gm*R/(gm-1);
##Energy balance across burner:
Q1cr=cp*(Tcr-Tt1);
f=(Q1cr/120000);
##total pressure loss:
dpt=1-Pcr/pt1;
print'%s %.4f %s'%('Total exit temperature if exit is choked in',Tt2,'K')
print'%s %.4f %s'%('Maximum heat released per unit mass of air in',Q1cr, 'kJ/kg')
print'%s %.7f %s'%('fuel-to-air ratio to thermally choke the combustor exit',f,'')
print'%s %.7f %s'%('Total pressure loss (in fraction)',dpt,'')
#calculate the new inlet mach no and spilled flow at the inlet
#initilization variable
import math
Tt1=50.+460. ##Converting the inlet temp. to the absolute scale i.e. in degree R
M1=0.5 ##Initial inlet Mach no.
pt1=14.7 ##Units in psia
gm=1.4 ##gamma
R=53.34 ##units in ft.lbf/lbm.degree R
Tcr=Tt1/0.69136
cp=gm*R/(gm-1)
##using energy equation:
Q1cr=cp*(Tcr-Tt1)
##since heat flux is 1.2(Q1cr).
q=1.2*Q1cr
Tt1cr1=Tt1+(Q1cr/cp) ##new exit total temp.
z=Tt1/Tt1cr1
M2=0.473
f=M1/(1+((gm-1)/2)*M1**2)**((gm+1)/(2*(gm-1)))
sm=((f*(M1)-f*(M2))/f*(M1))*100. ##sm=The % spilled flow at the inlet
print'%s %.5f %s'%('The new inlet Mach no.',M2,'')
print'%s %.5f %s'%('The % spilled flow at the inlet',sm,'')
#intilization variable
#calculate choking length abd exit mach no and total pressure loss and the static pressure and impulse due to friction
import math
d=0.2 ##diameter in meters.
l=0.2 ##length in meters.
Cf=0.005 ##average wall friction coefficient.
M1=0.24 ##inlet mach no.
gm=1.4 ##gamma.
##From FANNO tbale
L1cr=(9.3866*d/2)/(4*Cf);
L2cr=L1cr-l;
##from FANNO table
M2=0.3;
x=2.4956;
y=2.0351;
a=4.5383;
b=3.6191;
i1=2.043;
i2=1.698;
##% total pressure drop due to friction:
dpt=(x-y)/(x)*100;
##static pressur drop:
dps=(a-b)/a*100;
##Loss pf fluid:
lf=(i2-i1);
print'%s %.3f %s'%('The choking length of duct in',L1cr,'m')
print'%s %.1f %s'%('The exit Mach no.',M2,'')
print'%s %.6f %s'%('% total pressure loss',dpt,'')
print'%s %.5f %s'%('The static pressure drop in',dps,'%')
print'%s %.3f %s'%('Loss of impulse due to friction(I* times)',lf,'')
#initilization variable
import math
#caluclate maximum length of the duct that will support given in inlet condition and the new inlet condition and flow drop
M1=0.5
a=2. ## area of cross section units in cm^2
Cf=0.005 ##coefficient of skin friction
gm=1.4 ##gamma
##Calculations
c=2.*(2.+1.); ##Parameter of surface.
##From FANNO table: 4*Cf*L1cr/Dh=1.0691;
Dh=4.*a/c; ##Hydrolic diameter.
L1cr=1.069*Dh/(4.*Cf);
##maximum length will be L1cr.
##For new length(i.e. 2.16*L1cr), Mach no. M2 from FANNO table, M2=0.4;.
M2=0.4;
##the inlet total pressue and temp remains the same, therefore the mass flow rate in the duct is proportional to f(M):
f=0.5/(1.+((gm-1.)/2.)*0.5**2.)**((gm+1.)/(2.*(gm-1.)))
#endfunction
dm=(f*(M1)-f*(M2))/f*(M1)*100.+10;
print'%s %.3f %s'%("(a)Maximum length of duct that will support given inlet condition(in cm):",L1cr,"")
print'%s %.3f %s'%("(b)The new inlet condition mach no. M2:",M2,"")
print'%s %.3f %s'%("(c)% inlet mass flow drop due to the longer length of the duct:",dm,"")
import math
import numpy
M1=0.7;
dpt=0.99; ##pt2/pt1=dpt.
gm=1.4; ##gamma
A2=1.237
a=1/1.237;
import warnings
warnings.filterwarnings('ignore')
##Calculations:
k=(1./dpt)*(a)*(M1/(1.+(0.2*(M1)**2.))**3.);
po=([k*0.008,0,k*.12,0,k*.6,-1,k])
W=numpy.roots(po)
i=0;
s=1;
M2=W[4]
print -M2,"(a)The exit Mach no. M2:"
##p=p2/p1 i.e. static pressure ratio
p=dpt*((1.+(gm-1.)*(M1)**2./2.)/(1.+(gm-1.)*(M2)**2./2.))**(gm/(gm-1.))
##disp(p)
Cpr=(2./(gm*(M1)**2.))*(p-1.) ##Cpr is static pressure recovery : (p2-p1)/q1.
print"%s %.2f %s"%("(b)The static pressure recovery in the diffuser:",-Cpr,"")
##Change in fluid impulse:
##Fxwalls=I2-I1=A1p1(1+gm*M1**2)-A2p2(1+gm*M2**2)
##Let, u=Fxwall/(p1*A1)
u=1.+gm*(M1)**2.-(1.237)*(p)*(1.+(gm*(M2)**2.))
print"%s %.2f %s"%("(c)The force acting on the diffuser inner wall nondimensionalized by inlet static pressure and area:",-u,"")
print "Example2.13"
import numpy
M1=0.5 #inlet mach no.
p=10. #(p=pt1/p0) whaere pt1 is inlet total pressure and p0 is ambient pressure.
dpc=0.01 #dpc=(pt1-Pth)/pt1 i.e. total pressure loss in convergant section
f=0.99 #f=Pth/pt1
dpd=0.02 #dpd=(Pth-pt2)/Pth i.e. total pressure loss in the divergent section
j=1/0.98 #j=Pth/pt2
A=2. #a=A2/Ath. nozzle area expansion ratio.
gm=1.4 # gamma
R=287. #J/kg.K universal gas constant.
#Calculations:
#"th"" subscript denotes throat.
Mth=1. #mach no at thorat is always 1.
k=(j)*(1./A)*(Mth/(1+(0.2*(Mth)**2))**3)
po=([k*0.008,0,k*.12,0,k*.6,-1,k])
W=numpy.roots(po)
i=0;
s=1;
M2=W[4]
print M2,"(a)The exit Mach no. M2:"
#p2/pt2=1/(1+(gm-1)/2*M2**2)**(gm/(gm-1))
#pt2=(pt2/Pth)*(Pth/pt1)*(pt1/p0)*p0
#let pr=p2/p0
pr=((1/j)*f*p)/(1+(0.2*(M2)**2))**(gm/(gm-1))
print pr,"(b)The exit static pressure in terms of ambient pressure p2/p0:"#Fxwall=-Fxliquid=I1-I2
#let r=A1/Ath
r=(f)*(1/M1)*(((1+((gm-1)/2)*(M1)**2)/((gm+1)/2))**((gm+1)/(2*(gm-1))))
#disp(r)
#Psth is throat static pressure.
#z1=Psth/pt1=f/((gm+1)/2)**(gm/(gm-1))
z1=f/((gm+1)/2)**(gm/(gm-1))
#disp(z1)
#p1 is static pressure at inlet
#s1=p1/pt1
s1=1/(1+((gm-1)/2)*(M1)**2)**(gm/(gm-1))
#disp(s1)
#let y=Fxcwall/(Ath*pt1), where Fxwall is Fx converging-wall
y=s1*r*(1+(gm*(M1)**2))-(z1*(1+(gm*(Mth)**2)))
print y,"(c)The nondimensional axial force acting on the convergent nozzle:"
#similarly finding nondimensional force on the nozzle DIVERGENT section
#y1=Fxdiv-wall/Ath*pt1
#f1=p2/pt1
f1=pr*(1/p)
#disp(f1)
y1=z1*(1+(gm*(Mth)**2))-f1*A*(1+(gm*(M2)**2))
print y1,"(d)The nondimensional axial force acting on the divergent nozzle:"
#total axial force acting on nozzle wall: Fsum=y+y1
Fsum=y+y1
print Fsum,"(e)The total axial force(nondimensional) acting on the nozzle: "
import math
#calculate non dimensional axial force and negative sign on the axial force experienced by the compressor
p=20. ##p=p2/p1 i.e. compression ratio.
gm=1.4 ## gamma
##Vx1=Vx2 i.e. axial velocity remains same.
##calculations:
d=p**(1/gm) ##d=d2/d1 i.e. density ratio
A=1./d ## A=A2/A1 i.e. area ratio which is related to density ratio as: A2/A1=d1/d2.
##disp(A)
Fx=1.-p*A ##Fx=Fxwall/p1*A1 i.e nondimensional axial force.
print'%s %.7f %s'%("The non-dimensional axial force is :",Fx,"")
print'%s %.f %s'%("The negative sign on the axial force experienced by the compressor structure signifies a thrust production by this component.",Fx," ")
import math
print("Example 2.15")
t=1.8 ##t=T2/T1
d=1./t ##d=d2/d1 i.e. density ratio
v=1./d ##v=Vx2/Vx1 axial velocity ratio
ndaf=1.-(v) ##nondimensional axial force acting on the combustor walls
print'%s %.1f %s'%("The nondimensional axial force acting on the combustor walls:",ndaf,"")
print("Negative sign signifies a thrust production by the device")
import math
print("Example 2.16")
t=0.79 ##T2/T1 i.e. turbione expansion
gm=1.4 ##gamma
##calculations:
d=t**(1./(gm-1.))
##print'%s %.1f %s'%(d)
a=1./d ##area ratio
p=d**gm ##pressure ratio
ndaf=1.-p*a
print'%s %.2f %s'%("The nondimensional axial force:",ndaf,"")