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-14 14:55:03
 
//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 = enthalpy(Air;T = T_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 = entropy(Air;P = p_airinlet;H = h_airinlet) // Upstream point - air inlet - Downstream point - 2
// Comment = Isentropic reference
//Equation: 11
hs_2 = enthalpy(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 = temperature(Air;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 = entropy(burnt gases;P = p_3;H = h_3) // Upstream point - 3 - Downstream point - 4
s_3 = s_products(T_3;P_3;a_combustionchamber;lambda_combustionchamber)
// Comment = Isentropic reference
//Equation: 19
//hs_4 = enthalpy(burnt gases;P = p_4;S = s_3) // Downstream point - 4
hs_4 = h_products(Tis;a_combustionchamber;lambda_combustionchamber)
s_3 = s_products(Tis;P_4;a_combustionchamber;lambda_combustionchamber)
//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 = temperature(burnt gases;H = h_4) // Downstream point - 4
h_4 = h_products(T_4;a_combustionchamber;lambda_combustionchamber)
//Equation: 23
//s_4 = entropy(burnt gases;P = p_4;H = h_4) // Entropy
s4 = s_products(T_4;P_4;a_combustionchamber;lambda_combustionchamber)
// 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 = enthalpy(CH4;T = T_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)