Simple python function utilizing dictonary to calculate thermal cycle analysis of a turbojet engine. Simple but accurate turbojet thermodynamic model. visit www.sukhbindersingh.com for more info.
# Lets define all imports and matplotlib figure defaults
import numpy as np
import matplotlib.pyplot as plt
from __future__ import division
import math
import matplotlib
%matplotlib inline
matplotlib.rc('figure', figsize=(14, 7))
matplotlib.rc('font', weight='bold', size=12)
# Thermal cyclic analysis of turbojet engine
# define constants
prgcon={"llp":10,"llt":4,"cp":1005,"r":287,"g":9.81,"tlr":6.5,
"cfuel":45868.0,"prc":2.0,"tet":800,"pa":1.1013,"ta":291.15,
"masc":0.1289,"effc":0.95,"effcom":0.99,"plcom":0.05,
"bleed":0.01,"efmech":0.98,"effct":0.84,"effn":0.95,"powreq":1.0
}
def toyturbojet(prgcon):
result={"temps":[],"prs":[],"thrust":0.0,"spthrust":0.0,"sfc":0.0}
prc=prgcon["prc"]
t03=prgcon["tet"]
gam=(prgcon["cp"]/prgcon["r"])/(prgcon["cp"]/prgcon["r"]-1)
gamx=prgcon["cp"]/prgcon["r"]
gamy=(gam-1.0)/2.0
macfl=0.0
r1=1.0
pa=prgcon["pa"]
ta=prgcon["ta"]
p01=pa*r1*((1.0+gamy*macfl**2)**gamx)
t01=ta*(1.0+gamy*macfl**2)
den1=(p01*10000.0*prgcon["g"])/(prgcon["r"]*t01)
#print p01,t01,den1
trais=prc**(1.0/gamx)
tdrpis=t01*(trais-1.0)
tdrp=tdrpis/prgcon["effc"]
powerc=prgcon["masc"]*prgcon["cp"]*tdrp/1000.0
p02=p01*prc
t02=t01+tdrp
mascom=prgcon["masc"]-prgcon["bleed"]*prgcon["masc"]
tracom=t03-t02
fara1=1.0
if(tracom > 000 and tracom < 400):
fara1=0.99 *((t03-t02-10)*(1.+t02/3250.))/(prgcon["cfuel"]*prgcon["effcom"])
elif(tracom > 400 and tracom < 9000):
fara1=1.10*((t03-t02-50.)*(1.+t02/3250.))/(prgcon["cfuel"]*prgcon["effcom"])
masfuel=fara1*mascom
p03=p02*(1.0-prgcon["plcom"])
powert=powerc/prgcon["efmech"]
masct=mascom+masfuel
cp34=950.0+0.21*t03
tdrpct=powert/(masct*cp34/1000.0)
gamxc=cp34/prgcon["r"]
tdctis=tdrpct/prgcon["effct"]
t04is=t03-tdctis
pract=(t03/t04is)**gamxc
p04=p03/pract
t04=t03-tdrpct
#print p01,p02,p03,p04,t01,t02,t03,t04
rthst=prgcon["powreq"]
p05i=p04
t05=t04
cp45=950.0+0.21*t04
gamxp=cp45/prgcon["r"]
gamn=gamxp/(gamxp-1.0)
popan=p05i/pa
if p05i < 1.1*pa:
result["temps"]=([t01,t02,t03,t04,t05])
result["prs"]=([p01,p02,p03,p04,0.0])
result["thrust"]=0.0
result["sfc"]=0.0
result["spthrust"]=0.0
return result
emach=math.sqrt(2*(popan**(1/gamxp)-1)/(gamn-1))
t5i=t05/(1+0.5*(gamn-1)*emach*emach)
t5=t05-prgcon["effn"]*(t05-t5i)
vjet=math.sqrt(2*cp45*(t05-t5))
p05=pa*((t05/t5)**(gamn/(gamn-1)))
massj=masct+prgcon["bleed"]+prgcon["masc"]
thst=massj*vjet/prgcon["g"]
ssp=thst/prgcon["masc"]
sfc=masfuel/thst*3600.
result["temps"]=([t01,t02,t03,t04,t05])
result["prs"]=([p01,p02,p03,p04,p05])
result["thrust"]=thst
result["sfc"]=sfc
result["spthrust"]=ssp
return result
#print p01,p02,p03,p04,p05,t01,t02,t03,t04,t05,thst,ssp,sfc
Run and display Pressures and Temperature at each stage.
res=toyturbojet(prgcon)
print res
plt.plot(res["prs"])
plt.xticks([0,1,2,3,4])
plt.title("Pressure at each stage ")
plt.ylabel("Pressure (bar)")
{'temps': [291.15, 358.2351135385468, 800, 738.4911985475384, 738.4911985475384], 'thrust': 8.814671292959224, 'spthrust': 68.38379591124301, 'sfc': 0.5491223915272165, 'prs': [1.1013, 2.2026, 2.0924699999999996, 1.439652255801095, 1.4198539164535902]}
<matplotlib.text.Text at 0x10ee72a50>
plt.plot(res["temps"])
plt.xticks([0,1,2,3,4])
plt.title("Temperature at each stage ")
plt.ylabel("Teperature (K)")
<matplotlib.text.Text at 0x10eeae190>
pr1 = np.linspace(2.0,24.0,20)
tet1=np.linspace(800.0,1400.0,6)
pr= np.zeros((pr1.shape[0],tet1.shape[0]))
#print pr.shape[0]
sp=np.zeros((pr.shape[0],tet1.shape[0]))
sfc=np.zeros((pr.shape[0],tet1.shape[0]))
thr=np.zeros((pr.shape[0],tet1.shape[0]))
for j, tet in enumerate(tet1):
prgcon["tet"]=tet
pr[:,j]=pr1
for i,pp in enumerate(pr1):
prgcon["prc"]=pp
res =toyturbojet(prgcon)
#print i,pp
sp[i,j]=res["spthrust"]
sfc[i,j]=res["sfc"]
thr[i,j]=res["thrust"]
plt.subplot(3,1,1)
plt.plot(pr,sfc)
plt.title("sfc vs pr ")
plt.ylabel("sfc")
plt.xlabel("pr")
plt.legend="best"
plt.subplot(3,1,2)
plt.plot(pr,sp)
plt.title("Sp. Thrust vs pr ")
plt.ylabel("Sp. Thrust")
plt.xlabel("pr")
plt.subplot(3,1,3)
plt.legend="best"
plt.plot(pr,thr)
plt.title("Thrust vs pr ")
plt.ylabel("Thrust")
plt.xlabel("pr")
<matplotlib.text.Text at 0x10f89b390>
visit www.sukhbindersingh.com for more info.