#calcualte overspeed mach no
print("Example 5.1")
Md=1.5
##From isentropic table,
gm=1.4 ##gamma
A=1.176 ##A=A1/Ath=A1/Acr
##for same A, from isentropic table for M<1
My=0.61
##for My=0.61, from normal shock table
Mx=1.8
Mos=Mx
print'%s %.1f %s'%("Overspeed Mach no.",Mos,"")
import math
#calculate contractio ratio and the maximum pressure recovery
print("Example 5.2")
Md=2.65
Mx=Md
##for Mx=2.65, from normal shock table
My=0.4996
M1=My
##from isentropic table for M1=0.5,
A=1.34
##for Md=2.65, from isentropic table (A=A1/Acr)
A1=3.036
Af=A1/A
##from isentropic table Af,
Mth=2.35
##for Mth=2.35, from normal shock table
p=0.5615 ##p=pty/ptx
print'%s %.2f %s'%("Maximum total pressure recovery:",p,"")
#calculate inlet design contraction ratio and throat mach no
print("Example 5.3")
Md=3.3 ##from isentropioc table
A=5.629 ## A=A1/Acr=A1/Ath
Mx=Md ##from normal shock table
My=0.4596
M1=My
##from isentropic table
A11=1.425
pt=((1./A11-1./A)/(1./A))*100.
Af=A/A11
##for Af=3.95, from isentropic table for M>1
M1th=2.95
print'%s %.2f %s'%("Inlet design contraction ratio A1/Ath:",A," ")
print'%s %.2f %s'%("The % opening of the throat:",pt," ")
print'%s %.3f %s'%("Throat Mach no.:",M1th,"")
import math
#calculate inlet pressure recovery with the shock at the lip
print("Example 5.4")
M0=1.4
##from normal shock table
p=0.9582 ##p=pt2/pt0
M1=M0
##from isentropic table:
A=1.115 ##A=A1/Acr
A11=1.1 ##A11=Ax/A1
Af=A11*A
##from normal shock table for M>1
Mx=1.56
##from normal table
p1=0.91 ##p=pt2/pt0
p2=p
print'%s %.2f %s'%("(a)The best backpressure :",p,"")
print'%s %.3f %s'%("(b)The supercritical mode inlet total pressure recovery:",p1,"")
print'%s %.2f %s'%("(c)Inlet pressure recovery in subcritical mode with 10% spillage:",p2,"")
#calculate total pressure recovery of this inlet
import math
print("Example 5.5")
##th=theta and b=beta.
gm=1.4 ##gamma
##OBLIQUE SHOCK 1
M0=2.
th=8. ##degree
##from theta-beta-M chart,
b1=37. ##degree
Mn1=M0*math.sin(b1/57.3)
p1=0.993 ##p=pt2/pt1
Mn2=((2.+(gm-1.)*Mn1**2.)/(2.*gm*Mn1**2.-(gm-1.)))**(1/2.)
M2=Mn2/math.sin(b1-th/57.3)
##OBLIQUE SHOCK 2
M0=M2
th=12.
##from oblique shock chart,
b2=48.7
Mn1=M0*math.sin(b2/57.3)
p2=0.978
Mn2=((2.+(gm-1.)*Mn1**2)/(2.*gm*Mn1**2.-(gm-1.)))**(1/2.)
M3=Mn2/math.sin(b1-th/57.3)
##NORMAL SHOCK
M0=M3
b3=90.
pNS=0.977
Po=p1*p2*pNS
print'%s %.3f %s'%("Total pressure recovery:",Po,"")
import math
#calculate percent increase in gross thurst
print("Example 5.6")
M9=1. ## Mach no.
p=1/8. ##p=p0/pt7
gm=1.3 ##gamma
V9cd=(2.*(1.-p**((gm-1.)/gm)))**(1/2.)
px=p*((gm+1.)/2.)**(gm/(gm-1.))
V9c=(2.*(gm-1.)/(gm+1.))**(1/2.)
FR=(V9cd/V9c)/(1.+(1.-px)/gm)
pr=(FR-1.)*100./1.
print'%s %.3f %s'%("% increase in gross thrust:",pr," ")
#calculate velocites at various point and coefficent
import math
print("Example 5.7")
p98=0.95 ##p98=pt9/pt8
p87=0.98 ##p98=pt8/pt7
p70=8. ##p70=pt7/pt0
p97=8. ##p97=pt9/pt7
Cp=1243.7 ##specific heat in J/kg.K
gm=1.3 ##gamma
Tt9=900. ##Total temp. of the gas entering a C-D nozzle
Tt7=Tt9
p90=1. ##p90=p9/p0
p99=p98*p87*p70*p90 ##p99=pt9/p9
M9=(2./(gm-1.)*(p99**((gm-1.)/gm)-1.))**(1/2.) ##exit mach no.
T9=Tt9/(1.+(gm-1.)*M9**2/2.) ##The nozzle exit static temp.
a9=((gm-1.)*Cp*T9)**(1/2.) ##speed of sound in exit plane
V9=a9*M9 ##exit velocity
V9s=(2.*Cp*Tt7*(1.-p97**-((gm-1.)/gm)))**(1/2.)
p89=p87*p70*p90 ##p89=pt8/p9
V9i=(2.*Cp*Tt7*(1.-p89**-((gm-1.)/gm)))**(1/2.)
Cv=V9/V9i
print'%s %.1f %s'%("(a)V9 in",V9," m/s:")
print'%s %.1f %s'%("(b)V9s in",V9s," m/s:")
print'%s %.1f %s'%("(c)V9i in ",V9i,"m/s:")
print'%s %.3f %s'%("(d)The velocity coefficient Cv:",Cv,"")
%matplotlib inline
#plot the graphs
print "Example 5.8"
#calculate and graph the divergnece correction factor
import math
import numpy
import matplotlib
from matplotlib import pyplot
alfa=0 #alfa=cone half angle
dx=numpy.linspace(0,44,146)
x=numpy.zeros(146)
count=0;
for alfa in dx:
Ca=(1+math.cos(alfa*math.pi/180.))/2.; #Flow angularity loss coefficient
x[count]=Ca;
count=count+1;
#disp(Ca,"Divergence correction factor Ca:")
pyplot.plot(dx,x)
pyplot.title("Flow convergence loss in a conical nozzle")
pyplot.xlabel("Cone half-angle")
pyplot.ylabel("Flow angularity loss coefficient")
print "Example 5.9"
#plot the graphs
import math
import numpy
import matplotlib
from matplotlib import pyplot
%matplotlib inline
alfa=0.1
dx=numpy.linspace(0.1,44,88)
x=numpy.zeros(88)
g1=numpy.zeros(88)
count=0
g2=numpy.zeros(88)
gc1=0;
for alfa in dx:
Ca=(math.sin(alfa*math.pi/180.))/(alfa*math.pi/180.)
Cac=(1+math.cos(alfa*math.pi/180.))/2.
x[count]=Ca
count=count+1;
g1[gc1]=Cac;
gc1=gc1+1;
pyplot.plot(dx,g1)
pyplot.plot(dx,x)
pyplot.legend(["Conical","2D-CD"])
pyplot.xlabel("Divegent flap angle or Cone half-angle(degree)")
pyplot.ylabel("Flow angularity loss coefficient")
pyplot.title("Divergent loss of a conical nozzle and a 2D-CD nozzle")
print ("Example 5.10")
%matplotlib inline
#plot the graphs
import numpy
import matplotlib
from matplotlib import pyplot
p=0.96 #p=p't8/pt8
f=0.02
fAB=0.04
z0=numpy.linspace(0.45,0.63,7)
gmr=1.3/1.33 #gm=gm/gm' gm=gamma
gm=1.33
gm1=1.3
tlAB=7.
tl=6.
i=0;
z1=numpy.linspace(7,9,3)
for tlAB in z1:
tt=6.5
g1=numpy.zeros(7)
gc1=0;
for tt in z0:
A=(1+f+fAB)/(1+f)*((gmr)**(1./2))*1/p*((tlAB/(tl*tt))**(1./2))*((((gm1+1)/2.)**((gm1+1)/(2*(gm1-1))))/(((gm+1)/2.)**((gm+1)/(2.*(gm-1)))))
g1[gc1]=A
gc1=gc1+1;
number=0;
pyplot.plot(z0,g1)
i=i+1;
pyplot.xlabel("Turbine expansion parameter")
pyplot.ylabel("A8-AB-ON/A8-AB-OFF")
pyplot.title("Nozzle throat area variation with ")
pyplot.legend(["tau(AB)=7","tau(AB)=8","tau(AB)=9"])
print ("Example 5.12")
#plot the graphs
%matplotlib inline
import numpy
import matplotlib
from matplotlib import pyplot
import warnings
warnings.filterwarnings('ignore')
gm=1.1
M0=2.5
g1=numpy.zeros(40)
z0=numpy.linspace(0,4,40)
i=1;
z1=numpy.linspace(1.1,1.4,4)
for gm in z1:
gc1=0;
for M in z0:
p0=(1+(gm-1)/2*(M**2))**(gm/(gm-1))
p20=.4*p0-.5*p0
M=3
p42=0.37
NPR=p20*p42
g1[gc1]=p0
gc1=gc1+1;
pyplot.plot(z0,g1)
pyplot.title("Total-to-static pressure ratio")
pyplot.xlabel("Flight Mach no. (M0)")
pyplot.ylabel("pt0/p")
pyplot.legend(["gamma=1.1","gamma=1.2","gamma=1.3","gamma=1.4"])
i=i+1;
print "Example 5.13"
#plot the graphs
%matplotlib inline
import numpy
from numpy import linspace
import matplotlib
from matplotlib import pyplot
#T=Th/Tc
z0=numpy.linspace(0,8,160)
i=0
z1=numpy.linspace(1,4.5,7)
for T in z1:
g1=numpy.zeros(160);
gc1=0;
for alfa in z0:
FR=((1+alfa)**(1./2)*(T+alfa)**(1./2))/(T**(1./2)+alfa)
g1[gc1]=FR
gc1=gc1+1;
number=0;
pyplot.plot(z0,g1)
i=i+1;
pyplot.xlabel("Bypass ratio(alfa)")
pyplot.ylabel("Ratio of mixed to seperate-flow turbofan engines gross thrust")
pyplot.legend("T(hot)/T(cold)=1.5","T(hot)/T(cold)=2","T(hot)/T(cold)=2.5 so on")
pyplot.title("Ideal gross thrust gain with a perfect mixer")