Function LAMBD(T0;T;a)
//function to calculate the air factor knowing a (4 for CH4), T0 recatant initial temperature and T adiabatic temperature
//CH4  +  lambda (1+a/4) (O2 + 3.76 N2) <--> CO2  +  2 H2O  +  3.76 lambda (1+a/4) N2 + (lambda-1) (1+a/4)O2
	num= enthalpy(CH4;T=T0)-enthalpy(CO2;T=T)-a/2*enthalpy(H2O;T=T)+(1+a/4)*enthalpy(O2;T=T)
	denom=(1+a/4)*3,76*(enThalpy(N2;T=T)-enthalpy(N2;T=T0)) +(1+a/4)*(enthalpy(O2;T=T)-enthalpy(O2;T=T0))
	LAMBD= num/denom
END
Function h_products(T;a;lambda)
//function to calculate the enthalpy of reactants knowing a (4 for CH4),  air factor lambda and temperature T
//CH4  +  lambda (1+a/4) (O2 + 3.76 N2) <--> CO2  +  (1+a/4) H2O  +  3.76 lambda (1+a/4) N2 + (lambda-1) (1+a/4) O2
	nb_moles=(1+(1+a/4)+3,76*lambda*(1+a/4)+(lambda-1)*(1+a/4))
	molar_mass=MolarMass(CO2)/nb_moles+MolarMass(H2O)*(1+a/4)/nb_moles+MolarMass(N2)*3,76*lambda*(1+a/4)/nb_moles+MolarMass(O2)*(lambda-1)*(1+a/4)/nb_moles
	fract_CO2=1/nb_moles*MolarMass(CO2)/molar_mass
	fract_H2O=(1+a/4)/nb_moles*MolarMass(H2O)/molar_mass
	fract_N2=3,76*lambda*(1+a/4)/nb_moles*MolarMass(N2)/molar_mass
	fract_O2=(lambda-1)*(1+a/4)/nb_moles*MolarMass(O2)/molar_mass
	h_products=enThalpy(CO2;T=T)*fracT_CO2+enthalpy(H2O;T=T)*fract_H2O+enthalpy(N2;T=T)*fract_N2+enthalpy(O2;T=T)*fract_O2
END
Function s_products(T;P;a;lambda)
//function to calculate the enthalpy of reactants knowing a (4 for CH4),  air factor lambda and temperature T
//CH4  +  lambda (1+a/4) (O2 + 3.76 N2) <--> CO2  +  (1+a/4) H2O  +  3.76 lambda (1+a/4) N2 + (lambda-1) (1+a/4) O2
	nb_moles=(1+(1+a/4)+3,76*lambda*(1+a/4)+(lambda-1)*(1+a/4))
	molar_mass=MolarMass(CO2)/nb_moles+MolarMass(H2O)*(1+a/4)/nb_moles+MolarMass(N2)*3,76*lambda*(1+a/4)/nb_moles+MolarMass(O2)*(lambda-1)*(1+a/4)/nb_moles
	fract_CO2=1/nb_moles*MolarMass(CO2)/molar_mass
	fract_H2O=(1+a/4)/nb_moles*MolarMass(H2O)/molar_mass
	fract_N2=3,76*lambda*(1+a/4)/nb_moles*MolarMass(N2)/molar_mass
	fract_O2=(lambda-1)*(1+a/4)/nb_moles*MolarMass(O2)/molar_mass
	mixEntropy=-8,314/molar_mass*(1/nb_moles*ln(1/nb_moles)+(a/2)/nb_moles*ln((a/2)/nb_moles)+3,76*lambda*(1+a/4)/nb_moles*ln(3,76*lambda*(1+a/4)/nb_moles)+(lambda-1)*(1+a/4)/nb_moles*ln((lambda-1)*(1+a/4)/nb_moles))
	s0=entropy(CO2;T=T;P=P)*fract_CO2+entropy(H2O;T=T;P=P)*fract_H2O+entropy(N2;T=T;P=P)*fract_N2+entropy(O2;T=T;P=P)*fract_O2
	s_products= s0 +mixEntropy
END
//EQUATIONS
//Units: SI, Temperatures in Celsius, pressures in bar
//Project file: D:\_classement\_Thopt\THERMOPTIM_Pro_282\proj\GT_comb.prj
//Date and Time: 2024-08-24 19:08:57

//Flow rate unit: t/s

//GAS COMPOSITIONS
//burnt gases
//CO2	0.042254601489356794
//H2O	0.03459358954165343
//O2	0.1670338655210087
//N2	0.7438928324480963
//Ar	0.012225110999884659
//air
//N2	0.7555302216468832
//Ar	0.012416359476160373
//O2	0.2320534188769565
//CH4 ` methane
//CH4 ` methane	1.0

//PROCESSES

//Process: gas outlet
//Equation: 1
m_dot_gasoutlet = m_dot_turbine // Upstream process - turbine
//Equation: 2
m_dot_gasoutlet = 1.01827 // Given value

//Process: fuel
//Equation: 3
m_dot_fuel = 0.015643905 // Given value

//Process: air inlet
//Equation: 4
m_dot_airinlet = 1.0 // Given value
//Equation: 5
T_airinlet = 25.0// Given value (Celsius)
//Equation: 6
p_airinlet = 1.0// Given value (bar)
//Equation: 7
h_airinlet = calcH_TP("air";T = T_airinlet ;P = p_airinlet) // Downstream point - air inlet
//Equation: 8
m_dot_compressor = m_dot_airinlet //Flow propagation

//Process: compressor
//Equation: 9
m_dot_compressor = m_dot_airinlet // Upstream process - air inlet
//Equation: 10
s_airinlet = calcS_PH("air";P = p_airinlet;H = h_airinlet) // Upstream point - air inlet - Downstream point - 2
// Comment = Isentropic reference
//Equation: 11
hs_2 = calcH_PS("air";P = p_2;S = s_airinlet) // Downstream point - 2
//Equation: 12
etaT_compressor = 0.85// Isentropic efficiency
//Equation: 13
h_2 = h_airinlet + (hs_2 - h_airinlet)/etaT_compressor // Upstream point - air inlet - Downstream point - 2
//Equation: 14
T_2 = calcT_PH("air";P = p_2 ;H = h_2) // Downstream point - 2
// Comment = Given outlet pressure
//Equation: 15
p_2 = 16.0// Outlet pressure
//Equation: 16
W_dot_compressor = m_dot_compressor*(h_2 - h_airinlet) // DeltaH 

//Process: turbine
//Equation: 17
m_dot_turbine = m_dot_combustionchamber // Upstream process - combustion chamber
//Equation: 18
s_3 = calcS_PH("burnt gases";P = p_3;H = h_3) // Upstream point - 3 - Downstream point - 4
// Comment = Isentropic reference
//Equation: 19
hs_4 = calcH_PS("burnt gases";P = p_4;S = s_3) // Downstream point - 4
//Equation: 20
etaT_turbine = 0.85// Isentropic efficiency
//Equation: 21
h_4 = h_3 - etaT_turbine*(h_3 - hs_4) // Upstream point - 3 - Downstream point - 4
//Equation: 22
T_4 = calcT_PH("burnt gases";P = p_4 ;H = h_4) // Downstream point - 4
//Equation: 23
s_4 = calcS_PH("burnt gases";P = p_4;H = h_4) // Entropy
// Comment = Given outlet pressure
//Equation: 24
p_4 = 1.0// Outlet pressure
//Equation: 25
W_dot_turbine = m_dot_turbine*(h_4 - h_3) // DeltaH 

//Process: combustion chamber
// Comment = Calculate lambda simplified model oxidizer air, fuel CH4
//Equation: 26
T_3 = 1065.0// Given value (Celsius)
//Equation: 27
a_combustionchamber = 4// for CH4
//Equation: 28
lambda_combustionchamber = LAMBD(T_2;T_3;a_combustionchamber)// air factor lambda
//Equation: 29
h_3 = h_products(T_3;a_combustionchamber;lambda_combustionchamber)// enthalpy of the reactants
//Equation: 30
hfict_2 = h_products(T_2;a_combustionchamber;lambda_combustionchamber)// enthalpy of a fictitious inlet point for calculating the heat released
//Equation: 31
m_dot_combustionchamber = m_dot_compressor + m_dot_fuel // Upstream process - compressor - Fuel process fuel - Downstream process - combustion chamber
//Equation: 32
Q_dot_combustionchamber = (h_3 - hfict_2)*m_dot_combustionchamber // DeltaH
//Equation: 33
DeltaHr_combustionchamber = (-(-74850) +(-393520)+a_combustionchamber/2*(-242000))/16 // DeltaHr (kJ/kg) = (-(-74850) +(-393520) + a/2* (-242000))/16 for methane
//Equation: 34
m_dot_fuel = abs(Q_dot_combustionchamber/DeltaHr_combustionchamber) // fuel flow rate
// Comment = Isobaric process
//Equation: 35
p_3 = p_2// Isopressure
//Equation: 36
T_fuel = 15.0// Given value (Celsius)
//Equation: 37
p_fuel = 20.0// Given value (bar)
//Equation: 38
h_fuel = calcH_TP("CH4 ` methane";T = T_fuel ;P = p_fuel) // Fuel point - fuel

//NODES

//HEAT EXCHANGERS

//Number of equations: 38

//POINTS WITH SATURATION TEMPERATURE SET


//OTHER POINTS WITH PRESSURE SET


//SET FLOW RATES


//OVERALL BALANCE

//Equation: 39
useful_Energy = W_dot_compressor + W_dot_turbine
//Equation: 40
purchased_Energy = Q_dot_combustionchamber
//Equation: 41
eta_global = abs(useful_Energy/purchased_Energy)