Programming – Thermodynamics – Compressible flow with Nozzles and Diffusers

Dr. Stephen W. Crown

1) Excel Macro

Public Function VC(T As Double) As Double

K = 1.4

VC = (K * 287 * T) ^ 0.5

End Function

Public Function TT0(M As Double) As Double

K = 1.4

TT0 = 1 / (1 + (K - 1) / 2 * M ^ 2)

End Function

Public Function PP0(M As Double) As Double

K = 1.4

PP0 = 1 / (1 + (K - 1) / 2 * M ^ 2) ^ (K / (K - 1))

End Function

Public Function MPP0(PP01 As Double) As Double

K = 1.4

MPP0 = ((1 / PP01 ^ ((K - 1) / K) - 1) / ((K - 1) / 2)) ^ 0.5

End Function

Public Function MTT0(TT01 As Double) As Double

K = 1.4

MTT0 = ((1 / TT01 - 1) / ((K - 1) / 2)) ^ 0.5

End Function

Public Function AAS(M As Double) As Double

K = 1.4

AAS = 1 / M * (2 / (K + 1) * (1 + (K - 1) / 2 * M ^ 2)) ^ ((K + 1) / (2 * (K - 1)))

End Function

Public Function MY(MX As Double) As Double

K = 1.4

MY = ((MX ^ 2 + 2 / (K - 1)) / ((2 * K / (K - 1) * MX ^ 2) - 1)) ^ 0.5

End Function

Public Function P0YP0X(MX As Double) As Double

K = 1.4

MY1 = ((MX ^ 2 + 2 / (K - 1)) / ((2 * K / (K - 1) * MX ^ 2) - 1)) ^ 0.5

P0YP0X = MX / MY1 * ((1 + (K - 1) / 2 * MY1 ^ 2) / (1 + (K - 1) / 2 * MX ^ 2)) ^ ((K + 1) / (2 * (K - 1)))

End Function

Public Function MAAS(AAS1 As Double) As Double

K = 1.4

M = 1

S = 0.99

For I = 1 To 100

AAS2 = 1 / M * (2 / (K + 1) * (1 + (K - 1) / 2 * M ^ 2)) ^ ((K + 1) / (2 * (K - 1)))

DL = D

D = Abs(AAS1 - AAS2)

If DL < D Then S = -S / 2

If Abs(D) < AAS1 / 100000 Then Exit For

M = M + S

Next I

MAAS = M

End Function

Excel Spreadsheet Using Functions Defined in Macro

2) TI-83

ClrHome

1.4-->K: 1-->U: 0.287-->R: 0-->Q

Menu("--Main Menu--","Isentropic",M1,"Normal Shock",M2,"Stagnation Props",M3,"T-V-P-A-m",M4,"Quit",M5)

Lbl M5: Stop

Lbl M1

1-->Q

Menu("--Isentropic--","M Given",O1,"TTo Given",O2,"P/Po Given",O3,"A/A*(M<1) Given",O4,"A/A*(M>1) Given",O5)

Lbl M2

2-->Q

Menu("--Shock--","Mx Given",O6,"My Given",O7,"Py/Px Given",O8,"Ty/Tx Given",O9,"Poy/Pox Given",O0)

Lbl O1

Input "Mach(M)=",M

1/(1+(K-1)/2*M^2)-->T

1/(1+(K-1)/2*M^2)^(K/(K-1))-->P

1/M*(2/(K+1)*(1+(K-1)/2*M^2))^((K+1)/(2*(K-1)))-->A

Goto E1

Lbl O2

Input "T/To=",T

((1/T-1)*2/(K-1))^0.5-->M

1/(1+(K-1)/2*M^2)^(K/(K-1))-->P

1/M*(2/(K+1)*(1+(K-1)/2*M^2))^((K+1)/(2*(K-1)))-->A

Goto E1

Lbl O3

Input "P/Po(201/H)=",P

P^((K-1)/K)-->T

((1/T-1)*2/(K-1))^0.5-->M

1/M*(2/(K+1)*(1+(K-1)/2*M^2))^((K+1)/(2*(K-1)))-->A

Goto E1

Lbl O4

Input "A/A*(15/N)=",A

1-->M

0.99-->S

0-->D

For(I,1,100,1)

1/M*(2/(K+1)*(1+(K-1)/2*M^2))^((K+1)/(2*(K-1)))-->B

D-->F

abs(A-B)-->D

If F<D

-S/2-->S

If abs(D)<(A/100000)

Goto F1

M+S-->M

End

Lbl F1

1/(1+(K-1)/2*M^2)-->T

1/(1+(K-1)/2*M^2)^(K/(K-1))-->P

Goto E1

Lbl O5

Input "A/A*=",A

1-->M

-0.99-->S

0-->D

For(I,1,100,1)

1/M*(2/(K+1)*(1+(K-1)/2*M^2))^((K+1)/(2*(K-1)))-->B

D-->F

abs(A-B)-->D

If F<D

-S/2-->S

If abs(D)<(A/100000)

Goto F2

M+S-->M

End

Lbl F2

1/(1+(K-1)/2*M^2)-->T

1/(1+(K-1)/2*M^2)^(K/(K-1))-->P

Goto E1

Lbl O6

Input "Mx(M)=",M

((M^2+2/(K-1))/(2*K/(K-1)*M^2-1))^0.5-->Y

(1+(K-1)/2*M^2)/(1+(K-1)/2*Y^2)-->X

(1+K*M^2)/(1+K*Y^2)-->Z

M/Y*(X)^((K+1)/(2*(1-K)))-->W

O/W-->N

J*W-->H

Goto E1

Lbl M3

Input "Inlet T=",D

Input "Inlet P=",B

Input "Inlet V=",F

Input "A*=",N

N-->O

(K*R*D*U*1000)^0.5-->C

F/C-->M

D*(1+(K-1)/2*M^2)-->G

B*(1+(K-1)/2*M^2)^(K/(K-1))-->H

H-->J

Disp "C=",C,"M=",M

Pause

Disp "To=",G,"Po=",H

Pause

Goto E1

Lbl M4

T*G-->D

M*(K*R*D*U*1000)^0.5-->B

P*H-->S

N*A-->F

F*B*S/(R*U*D)/10000-->I

Disp "Given A*=",N

Disp "Given Po=",H

Disp "T=",D,"V=",B

Pause

Disp "P=",S,"A=",F,"m_dot=",I

Pause

Goto E1

Lbl O7: Lbl O8: Lbl O9: Lbl O0

Disp "Under-Constr.","","You-are-free","to-fix-this."

Pause

Lbl E1

ClrHome

Disp ""

If Q=1

Then

Disp "M=",M,"T/To=",T

Pause

Disp "P/Po=",P,"A/A*=",A

Pause

End

If Q=2

Then

Disp "Mx=",M,"My=",Y,"Py/Px=",Z

Pause

Disp "Ty/Tx=",X,"Poy/Pox=",W,"A*y=",N

Pause

End

3) TI-89

()

Prgm

Local qqm,qqq,k,aas2,i,s,d,dl

setMode("Exact/Approx","Approximate")

ClrIO

Disp ""

1.4-->k: 0.287-->r: 1-->uc

PopUp {"Stagnation Props","Isentropic","Normal Shock","Props (TVPAm)","Quit"},qqq

If qqq=5: Return

If qqq=1 Then

Input "Inlet T=",t1

Input "Inlet P=",p1

Input "Inlet V=",v1

Input "A*=",as

as-->asx

(k*r*t1*uc*1000)^0.5-->cv1

v1/cv1-->m1

1/(1+(k-1)/2*m1^2)-->tto

1/(1+(k-1)/2*m1^2)^(k/(k-1))-->ppo

1/m1*(2/(k+1)*(1+(k-1)/2*m1^2))^((k+1)/(2*(k-1)))-->aas

t1/tto-->to

p1/ppo-->po

po-->pox

Disp "C="&string(cv1),"M="&string(m1),"To="&string(to),"Po="&string(po)

Goto end1

EndIf

If qqq=2

PopUp {"Given M","Given TTo","Given P/Po","Given A/A*(M<1)","Given A/A*(M>1)"},qqm

If qqq=3 Then

PopUp {"Given Mx","Given My","Given Py/Px","Given Ty/Tx","Given Poy/Pox"},qqm

qqm+5-->qqm

EndIf

If qqq=4 Then

Disp "A*="&string(as),"Po="&string(po)

Disp "Calculated w/above Props"

Pause

1/(1+(k-1)/2*m^2)*to-->tt

(k*r*tt*uc*1000)^0.5*m-->vt

1/(1+(k-1)/2*m^2)^(k/(k-1))*po-->pt

aas*as-->at

(at*vt*pt)/(r*tt*10000)-->md

Disp "T="&string(tt),"V="&string(vt),"P="&string(pt),"A="&string(at),"m_dot(md)="&string(md)

Goto end1

EndIf

If qqm=1 Then

Input "Mach Number=",m

1/(1+(k-1)/2*m^2)-->tto

1/(1+(k-1)/2*m^2)^(k/(k-1))-->ppo

1/m*(2/(k+1)*(1+(k-1)/2*m^2))^((k+1)/(2*(k-1)))-->aas

EndIf

If qqm=2 Then

Input "T/To=",tto

(((1/tto-1)*2)/(k-1))^0.5-->m

1/(1+(k-1)/2*m^2)^(k/(k-1))-->ppo

1/m*(2/(k+1)*(1+(k-1)/2*m^2))^((k+1)/(2*(k-1)))-->aas

EndIf

If qqm=3 Then

Input "P/Po=",ppo

ppo^((k-1)/k)-->tto

(((1/tto-1)*2)/(k-1))^0.5-->m

1/m*(2/(k+1)*(1+(k-1)/2*m^2))^((k+1)/(2*(k-1)))-->aas

EndIf

If qqm=4 Then

Input "A/A*=",aas

1-->m: 0.99-->s: 0-->d

For i,1,100,1

1/m*(2/(k+1)*(1+(k-1)/2*m^2))^((k+1)/(2*(k-1)))-->aas2

d-->dl

abs(aas-aas2)-->d

If dl<d

-s/2-->s

If abs(d)<aas/100000

exitfor

m+s-->m

EndFor

1/(1+(k-1)/2*m^2)-->tto

1/(1+(k-1)/2*m^2)^(k/(k-1))-->ppo

EndIf

If qqm=5 Then

Input "A/A*=",aas

1-->m: -0.99-->s: 0-->d

For i,1,100,1

1/m*(2/(k+1)*(1+(k-1)/2*m^2))^((k+1)/(2*(k-1)))-->aas2

d-->dl

abs(aas-aas2)-->d

If dl<d

-s/2-->s

If abs(d)<aas/100000

exitfor

m+s-->m

EndFor

1/(1+(k-1)/2*m^2)-->tto

1/(1+(k-1)/2*m^2)^(k/(k-1))-->ppo

EndIf

If qqm=6 Then

Input "Mach Number(Mx)=",mx

((mx^2+2/(k-1))/((2k)/(k-1)*mx^2-1))^0.5-->my

(1+(k-1)/2*mx^2)/(1+(k-1)/2*my^2)-->tytx

(1+k*mx^2)/(1+k*my^2)-->pypx

mx/my*tytx^((k+1)/(2*(1-k)))-->poypox

pox*poypox-->po

as/poypox-->as

EndIf

If qqm=7 Then

Disp "Under Construction", "Given My"

Pause

EndIf

If qqm=8 Then

Disp "Under Construction", "Given Py/Px"

Pause

EndIf

If qqm=9 Then

Disp "Under Construction", "Given Ty/Tx"

Pause

EndIf

If qqm=10 Then

Disp "Under Construction", "Poy/Pox"

Pause

EndIf

ClrIO

Disp ""

If qqm<6 Then

m-->x1: tto-->x2: ppo-->x3: aas-->x4

Disp "X1) M="&string(x1),"X2) T/To="&string(x2)

Disp "X2) P/Po="&string(x3),"X4) A/A*="&string(x4)

Else

mx-->x5: my-->x6: pypx-->x7: tytx-->x8:poypox-->x9

as-->x10: po-->x11

Disp "X5) Mx="&string(x5),"X6) My="&string(x6),"X7) Py/Px="&string(x7)

Disp "X8) Ty/Tx="&string(x8),"X9) Poy/Pox="&string(x9)

Pause

Disp "X10) A*y="&string(x10),"X11) Poy="&string(x11)

EndIf

Lbl end1

Pause

t2()

EndPrgm