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\regen_GT_CH4.prj
//Date and Time: 2024-08-24 19:26:25

//Flow rate unit: -

//GAS COMPOSITIONS
//burnt gases
//burnt_gas
//CO2	0.04419006337631066
//H2O	0.03617814048485812
//O2	0.1640556604215191
//N2	0.7433597848138069
//Ar	0.012216350903505186
//air
//N2	0.7555302216468832
//Ar	0.012416359476160373
//O2	0.2320534188769565
//CH4 ` methane
//CH4 ` methane	1.0

//PROCESSES

//Process: regen gas
//Equation: 1
m_dot_regengas = m_dot_turbine // Upstream process - turbine
// Comment = isobaricExchange
//Equation: 2
p_5 = p_4 // Upstream point - 4 - Downstream point - 5
// Comment = //Exchange process connected to a heat exchanger

//Process: regen air
//Equation: 3
m_dot_regenair = m_dot_compressor // Upstream process - compressor
// Comment = isobaricExchange
//Equation: 4
p_2bis = p_2 // Upstream point - 2 - Downstream point - 2 bis
// Comment = //Exchange process connected to a heat exchanger

//Process: gas outlet
//Equation: 5
m_dot_gasoutlet = m_dot_regengas // Upstream process - regen gas
//Equation: 6
m_dot_gasoutlet = 1.0163722 // Given value

//Process: air inlet
//Equation: 7
m_dot_airinlet = 1.0 // Given value
//Equation: 8
T_airinlet = 15.0// Given value (Celsius)
//Equation: 9
p_airinlet = 1.0// Given value (bar)
//Equation: 10
h_airinlet = calcH_TP("air";T = T_airinlet ;P = p_airinlet) // Downstream point - air inlet
//Equation: 11
m_dot_compressor = m_dot_airinlet //Flow propagation

//Process: fuel
//Equation: 12
m_dot_fuel = 0.016372202 // Given value
//Equation: 13
T_fuel = 15.0// Given value (Celsius)
//Equation: 14
p_fuel = 20.0// Given value (bar)
//Equation: 15
h_fuel = calcH_TP("CH4 ` methane";T = T_fuel ;P = p_fuel) // Downstream point - fuel

//Process: compressor
//Equation: 16
m_dot_compressor = m_dot_airinlet // Upstream process - air inlet
//Equation: 17
s_airinlet = calcS_PH("air";P = p_airinlet;H = h_airinlet) // Upstream point - air inlet - Downstream point - 2
// Comment = Polytropic reference
//Equation: 18
ds_2 = (1 - etaT_compressor)/etaT_compressor*8.314/M_airinlet*ln(p_2/p_airinlet) // Upstream point - air inlet - Downstream point - 2
//Equation: 19
s_2 = s_airinlet + ds_2 // Entropy - 2
//Equation: 20
M_airinlet = 28.957763399999997 // Molar mass - 2
//Equation: 21
etaT_compressor = 0.85// Polytropic efficiency
// Comment = Polytropic coefficient: k = -Math.log(aval.p/amont.p)/Math.log(aval.V/amont.V)
//Equation: 22
h_2 = calcH_PS("air";P = p_2;S = s_2) // Enthalpy
//Equation: 23
T_2 = calcT_PH("air";P = p_2 ;H = h_2) // Downstream point - 2
// Comment = Given outlet pressure
//Equation: 24
p_2 = 16.0// Outlet pressure
//Equation: 25
W_dot_compressor = m_dot_compressor*(h_2 - h_airinlet) // DeltaH 

//Process: turbine
//Equation: 26
m_dot_turbine = m_dot_combustionchamber // Upstream process - combustion chamber
//Equation: 27
s_3 = calcS_PH("burnt_gas";P = p_3;H = h_3) // Upstream point - 3 - Downstream point - 4
// Comment = Polytropic reference
//Equation: 28
ds_4 = -(1 - etaT_turbine)*8.314/M_3*ln(p_4/p_3) // Upstream point - 3 - Downstream point - 4
//Equation: 29
s_4 = s_3 + ds_4 // Entropy - 4
//Equation: 30
M_3 = 28.587049318748424 // Molar mass - 4
//Equation: 31
etaT_turbine = 0.85// Polytropic efficiency
//Equation: 32
h_4 = calcH_PS("burnt_gas";P = p_4;S = s_4) // Enthalpy
// Comment = Polytropic coefficient: k = -Math.log(aval.p/amont.p)/Math.log(aval.V/amont.V)
//Equation: 33
h_4 = calcH_PS("burnt_gas";P = p_4;S = s_4) // Enthalpy
//Equation: 34
T_4 = calcT_PH("burnt_gas";P = p_4 ;H = h_4) // Downstream point - 4
// Comment = Given outlet pressure
//Equation: 35
p_4 = 1.0// Outlet pressure
//Equation: 36
W_dot_turbine = m_dot_turbine*(h_4 - h_3) // DeltaH 

//Process: combustion chamber
// Comment = Calculate lambda simplified model oxidizer air, fuel CH4
//Equation: 37
T_3 = 1150.0// Given value (Celsius)
//Equation: 38
a_combustionchamber = 4// for CH4
//Equation: 39
lambda_combustionchamber = LAMBD(T_2bis;T_3;a_combustionchamber)// air factor lambda
//Equation: 40
h_3 = h_products(T_3;a_combustionchamber;lambda_combustionchamber)// enthalpy of the reactants
//Equation: 41
hfict_2bis = h_products(T_2bis;a_combustionchamber;lambda_combustionchamber)// enthalpy of a fictitious inlet point for calculating the heat released
//Equation: 42
m_dot_combustionchamber = m_dot_regenair + m_dot_fuel // Upstream process - regen air - Fuel process fuel - Downstream process - combustion chamber
//Equation: 43
Q_dot_combustionchamber = (h_3 - hfict_2bis)*m_dot_combustionchamber // DeltaH
//Equation: 44
DeltaHr_combustionchamber = (-(-74850) +(-393520)+a_combustionchamber/2*(-242000))/16 // DeltaHr (kJ/kg) = (-(-74850) +(-393520) + a/2* (-242000))/16 for methane
//Equation: 45
m_dot_fuel = abs(Q_dot_combustionchamber/DeltaHr_combustionchamber) // fuel flow rate
// Comment = Isobaric process
//Equation: 46
p_3 = p_2bis// Isopressure
//Equation: 47
T_fuel = 15.0// Given value (Celsius)
//Equation: 48
p_fuel = 20.0// Given value (bar)
//Equation: 49
h_fuel = calcH_TP("CH4 ` methane";T = T_fuel ;P = p_fuel) // Fuel point - fuel

//NODES

//HEAT EXCHANGERS

//Heat exchanger: regenerator
//Equation: 50
mCp_regenair = (h_2bis - h_2)/(T_2bis - T_2)*m_dot_regenair // mCpf =deltaH/deltaT - regen air
//Equation: 51
mCp_regengas = (h_5 - h_4)*m_dot_regengas/(T_5 - T_4) // mCpc =-deltaH/deltaT - regen gas
// Comment = mCpc>mCpf
//Equation: 52
UA_regenerator = NTU_regenerator *mCp_regenair // Cold fluid - regen air
//Equation: 53
R_regenerator = mCp_regenair /mCp_regengas // Hot fluid - regen gas - Cold fluid - regen air
//Equation: 54
T_2bis = T_2 + epsilon_regenerator*(T_4 - T_2) // Hot fluid outlet temperature
//Equation: 55
h_2bis =  calcH_TP("air";T = T_2bis;P = p_2bis)// Enthalpy
//Equation: 56
h_5 = - m_dot_regenair /m_dot_regengas*(h_2bis - h_2) + h_4 // Hot fluid - regen gas - Cold fluid - regen air
//Equation: 57
T_5 =  calcT_PH("burnt_gas";P = p_5;H = h_5)// Hot fluid outlet temperature
// Comment = epsilon given value
//Equation: 58
epsilon_regenerator = 0.8499764599907744 // Given value
//Equation: 59
argLn_regenerator = (1 - epsilon_regenerator*R_regenerator)/(1 - epsilon_regenerator)
//Equation: 60
NTU_regenerator = 1/(1 - R_regenerator)*ln(argLn_regenerator) // Counterflow heat exchanger
//Equation: 61
Q_dot_regengas = m_dot_regengas*(h_5 - h_4) // DeltaH hot fluid
//Equation: 62
Q_dot_regenair = m_dot_regenair*(h_2bis - h_2) // DeltaH cold fluid

//Number of equations: 62

//POINTS WITH SATURATION TEMPERATURE SET


//OTHER POINTS WITH PRESSURE SET


//SET FLOW RATES


//OVERALL BALANCE

//Equation: 63
useful_Energy = W_dot_compressor + W_dot_turbine
//Equation: 64
purchased_Energy = Q_dot_combustionchamber
//Equation: 65
eta_global = abs(useful_Energy/purchased_Energy)