"Fuel Air Cycle Analysis - Stoichiometric Only at this Point.  Otherwise some problems with EES
arise"

"Fuel is n-octane.  C8H18"

"This calculation uses the chemical equilibrium module to find the composition of
  the gas right after combustion."

y = 18/8
epsilon = 4/(4+y)
phi = 1.0
FA_s = (12.011+1.008*y)/(34.56*(4+y))
FA = phi*FA_s
{x_b = 0.08}
M_f = MOLARMASS(C8H18)
psi = 3.773

"Composition Based on 1 mole of O2 - Please see table 4.4 in Text lean column"

n_f = 4*(1-x_b)*(1+2*epsilon)*phi/M_f
n_O2 = 1-x_b*phi
n_N2 = psi
n_CO2 = x_b*epsilon*phi
n_H2O = 2*x_b*(1-epsilon)*phi
n_u= n_f+n_O2+n_N2+n_CO2+n_H2O
y_f = n_f/n_u
y_O2 = n_O2 /n_u
y_N2 = n_N2/n_u
y_CO2 = n_CO2 / n_u
y_H2O = n_H2O / n_u
m_RP = 32+4*phi*(1+2*epsilon)+28.16*psi
m_air = 32+28.16*psi
M_u = m_RP/n_u
R_u = R# / M_u

"Engine Data"

r_c = 8

"Reference internal energy and entropy"

T_0 = 298.15									"[K]"
P_0 = Po#

ubar_0 = y_f*INTENERGY(C8H18,T=T_0)+y_O2*INTENERGY(O2,T=T_0)+y_N2*INTENERGY(N2,T=T_0)+y_CO2*INTENERGY(CO2,T=T_0)+y_H2O*INTENERGY(H2O,T=T_0)
{sbar_0 = y_f*ENTROPY(C8H18,T=T_0,P=y_f*P_0)+y_O2*ENTROPY(O2,T=T_0,P=y_O2*P_0)+y_N2*ENTROPY(N2,T=T_0,P=y_N2*P_0)+y_CO2*ENTROPY(CO2,T=T_0,P=y_CO2*P_0)+y_H2O*ENTROPY(H2O,T=T_0,P=y_H2O*P_0)}

"Property calculations at starting point"

T_1 = 350
P_1 = Po#
P_1 * v_1 = R_u * T_1

"Sensible Internal Energy at 1"

ubar_1 = y_f*INTENERGY(C8H18,T=T_1)+y_O2*INTENERGY(O2,T=T_1)+y_N2*INTENERGY(N2,T=T_1)+y_CO2*INTENERGY(CO2,T=T_1)+y_H2O*INTENERGY(H2O,T=T_1)-ubar_0

u_1 = ubar_1/M_u
u_1air = u_1*m_RP/m_air

"Sensible Entropy at 1"

sbar_1 = y_f*ENTROPY(C8H18,T=T_1,P=y_f*P_1)+y_O2*ENTROPY(O2,T=T_1,P=y_O2*P_1)+y_N2*ENTROPY(N2,T=T_1,P=y_N2*P_1)+y_CO2*ENTROPY(CO2,T=T_1,P=y_CO2*P_1)+y_H2O*ENTROPY(H2O,T=T_1,P=y_H2O*P_1)

s_1= sbar_1 / M_u

"Properties at 2"

v_2 = v_1 / r_c
P_2 * v_2 = R_u * T_2

"Entropy at 2"

sbar_2 = y_f*ENTROPY(C8H18,T=T_2,P=y_f*P_2)+y_O2*ENTROPY(O2,T=T_2,P=y_O2*P_2)+y_N2*ENTROPY(N2,T=T_2,P=y_N2*P_2)+y_CO2*ENTROPY(CO2,T=T_2,P=y_CO2*P_2)+y_H2O*ENTROPY(H2O,T=T_2,P=y_H2O*P_2)

sbar_2 = sbar_1
s_2 = sbar_2/M_u

"Sensible Internal Energy at 2"

ubar_2 = y_f*INTENERGY(C8H18,T=T_2)+y_O2*INTENERGY(O2,T=T_2)+y_N2*INTENERGY(N2,T=T_2)+y_CO2*INTENERGY(CO2,T=T_2)+y_H2O*INTENERGY(H2O,T=T_2)-ubar_0

u_2 = ubar_2/M_u
u_2air = u_2*m_RP/m_air

"Specific work done in compressing from 1 to 2"

w_12 = u_2 - u_1
w_12air = w_12 * m_RP/m_air

"Now do constant volume combustion.  Start by calculating the enthalpy of formation"

Delu_CO2 = -393500
Delu_H2O = -240600
Delu_C8H18 = -204300

Delubar = y_f*Delu_C8H18 + y_CO2*Delu_CO2+y_H2O*Delu_H2O
Delu = Delubar/M_u
Delu_air = Delu*m_RP/m_air
Delu_test = -118.2 - 2956*x_b


alpha = 8
beta = 18
as = 12.5

AO = 1.e-5
CO = alpha/(2*as/phi)
HO = beta/(2*as/phi)
NO=3.773

"Chemical Equilibrium Burned Gas"
{T_3 = 2700
P_3 = 7000}

CALL CHEM_EQUIL(P_3,T_3,AO,CO,HO,NO:x_H2,x_O2,x_H2O,x_CO,x_CO2,x_OH,x_H,x_O,x_N2,x_N,x_NO,x_NO2,x_CH4,x_A)

"Molecular weight of burned gas residual."

Mw=x_CO*MOLARMASS(CO)+x_CO2*MOLARMASS(CO2)+x_H*1+x_H2*MOLARMASS(H2)+x_H2O*MOLARMASS(H2O)+x_N*7+x_N2*MOLARMASS(N2)+x_NO*MOLARMASS(NO)+x_O*16+x_O2*MOLARMASS(O2)+x_OH*17+x_NO2*MOLARMASS(NO2)+x_A*MOLARMASS(Argon)+x_CH4*MOLARMASS(CH4)

"Properties of  gas components at T_3 from the JANAF tables"

  	CALL JANAF('CO',T_3:CP_CO,H_CO,S0_CO)
	CALL JANAF('CO2',T_3:CP_CO2,H_CO2,S0_CO2)
	CALL JANAF('H',T_3:CP_H,H_H,S0_H)
	CALL JANAF('H2',T_3:CP_H2,H_H2,S0_H2)
	CALL JANAF('H2O',T_3:CP_H2O,H_H2O,S0_H2O)
	CALL JANAF('N',T_3:CP_N,H_N,S0_N)
	CALL JANAF('N2',T_3:CP_N2,H_N2,S0_N2)
	CALL JANAF('NO',T_3:CP_NO,H_NO,S0_NO)
	CALL JANAF('O',T_3:CP_O,H_O,S0_O)
	CALL JANAF('O2',T_3:CP_O2,H_O2,S0_O2)
	CALL JANAF('OH',T_3:CP_OH,H_OH,S0_OH)
	CALL JANAF('NO2',T_3:CP_NO2,H_NO2,S0_NO2)
	CALL JANAF('Ar',T_3:CP_Ar,H_Ar,S0_Ar)
	CALL JANAF('CH4',T_3:CP_CH4,H_CH4,S0_CH4)

"Properties of Burned Gas Mixture.  Enthalpy and Internal Energy"

hbar_3=x_CO*h_CO+x_CO2*h_CO2+x_H*h_H+x_H2*h_H2+x_H2O*h_H2O+x_N*h_N+x_N2*h_N2+x_NO*h_NO+x_O*h_O+x_O2*h_O2+x_OH*h_OH+x_NO2*h_NO2+x_A*h_Ar+x_CH4*h_CH4

ubar_3 = hbar_3 - R# * T_3
R_b = R#/Mw
u_3 = ubar_3/M_u
u_3air = u_3*m_RP/m_air

u_3 = u_2 + Delu

v_3 = v_2
P_3 * v_3 = R_b * T_3

"Entropy of Burned Gas Mixture at 3"

S_CO = S0_CO - R# * ln( x_CO * P_3 / Po#)
S_CO2 = S0_CO2 - R# * ln( x_CO2 * P_3 / Po#)
S_H = S0_H - R# * ln( x_H * P_3 / Po#)
S_H2 = S0_H2 - R# * ln( x_H2 * P_3 / Po#)
S_H2O = S0_H2O - R# * ln( x_H2O * P_3 / Po#)
S_N = S0_N - R# * ln( x_N * P_3 / Po#)
S_N2 = S0_N2 - R# * ln( x_N2 * P_3 / Po#)
S_NO = S0_NO - R# * ln( x_NO * P_3 / Po#)
S_O = S0_O - R# * ln( x_O * P_3 / Po#)
S_O2 = S0_O2 - R# * ln( x_O2 * P_3 / Po#)
S_OH = S0_OH - R# * ln( x_OH * P_3 / Po#)
S_NO2 = S0_NO2 - R# * ln( x_NO2 * P_3 / Po#)
S_Ar = S0_Ar - R# * ln( x_A * P_3 / Po#)
trick=x_CH4 * P_3 / Po#
S_CH4 = S0_CH4 - R# * ln( trick)

sbar_3=x_CO*S_CO+x_CO2*S_CO2+x_H*S_H+x_H2*S_H2+x_H2O*S_H2O+x_N*S_N+x_N2*S_N2+x_NO*S_NO+x_O*S_O+x_O2*S_O2+x_OH*S_OH+x_NO2*S_NO2+x_A*S_Ar+x_CH4*S_CH4
s3 = sbar_3/Mw
s3_air = s3* m_RP/m_air

"We now have a picture of what's happened after combustion.  Last step.  Isentropic expansion
  to Point 4.  We start with the burned mixture composition in Table 4-3."

n_bCO2 = epsilon*phi
n_bH2O = 2*(1-epsilon)*phi
n_bO2 = 1-phi
n_bN2 = psi
n_b = (1-epsilon)*phi+1+psi

y_bCO2 = n_bCO2 /n_b
y_bH2O = n_bH2O/n_b
y_bO2 = n_bO2/n_b
y_bN2 = n_bN2/n_b

"Entropy at 4"

sbar_4 = y_bN2*ENTROPY(N2,T=T_4,P=y_bN2*P_4)+y_bCO2*ENTROPY(CO2,T=T_4,P=y_bCO2*P_4)+y_bH2O*ENTROPY(H2O,T=T_4,P=y_bH2O*P_4)
s4 = sbar_4/M_b
s4_air = s4* m_RP/m_air

s4=s3

v_4 = v_1
M_b = y_bO2*MOLARMASS(O2)+y_bN2*MOLARMASS(N2)+y_bH2O*MOLARMASS(H2O)+y_bCO2*MOLARMASS(CO2)

R_4 = R# / M_b
P_4 * v_4 = R_4 * T_4

ubar_4 = y_bN2*INTENERGY(N2,T=T_4)+y_bCO2*INTENERGY(CO2,T=T_4)+y_bH2O*INTENERGY(H2O,T=T_4)
u_4 = ubar_4/M_b
u_4air = u_4*m_RP/m_air

w_34 = u_3-u_4
w_34air = w_34 * m_RP/m_air


"Entropy at 5"

P_5 =Po#
sbar_5 = y_bN2*ENTROPY(N2,T=T_5,P=y_bN2*P_5)+y_bCO2*ENTROPY(CO2,T=T_5,P=y_bCO2*P_5)+y_bH2O*ENTROPY(H2O,T=T_5,P=y_bH2O*P_5)
s5 = sbar_5/M_b
s5_air = s5* m_RP/m_air

s5=s4
P_5 * v_5 = R_4 * T_5

x_b = v_2/v_5					"This starts iteration to get correct residual fraction"

"Calculation of the Goodies!!"

w_c = w_34 - w_12
w_cair = w_c* m_RP/m_air

Q_LHV = 44400
mfm = (1-x_b)/(1+1/FA)
q_add = mfm*Q_LHV

eta_f = w_c / q_add

eta_f1 = 1 - 1/r_c^(0.3)

imep = w_c / v_1