#*---------------------------------------------------------------------
* EMSO Model Library (EML) Copyright (C) 2004 - 2009 ALSOC.
*
* This LIBRARY is free software; you can distribute it and/or modify
* it under the therms of the ALSOC FREE LICENSE as available at
* http://www.enq.ufrgs.br/alsoc.
*
* EMSO Copyright (C) 2004 - 2009 ALSOC, original code
* from http://www.rps.eng.br Copyright (C) 2002-2004.
* All rights reserved.
*
* EMSO is distributed under the therms of the ALSOC LICENSE as
* available at http://www.enq.ufrgs.br/alsoc.
*
*----------------------------------------------------------------------
* Equilibrium modeling of a biomass gasifier
*----------------------------------------------------------------------
*
*   Description:
*       Thermodynamic equilibrium modeling of a biomass gasifier.
*
*   Assumptions:
*       * thermodynamic equilibrium
*		* steady-state
*		* ideal gas relations
*		* ambient and atmospheric input conditions (To,Po)
*		* global gasification reaction:
*			CHONS + H2O + Air -> CO + CO2 + CH4 + H2 
*				+ H2O + O2 + N2 + SO2
*		* all oxygen is consumed in the process
*
*	Specify:
*		* ultimate biomass analysis (dry.massfrac)
*		* moisture biomass ratio (moisture.massfrac)
*		* equivalence ratio (phi)
*		* relative air humidity (air.rh)
*		* equilibrium temperature (Teq)
*
*----------------------------------------------------------------------
* Author: Rodolfo Rodrigues and Argimiro R. Secchi
* GIMSCOP/UFRGS - Group of Integration, Modeling, Simulation, 
*					Control, and Optimization of Processes
*   LASIM/UFRGS - Simulation Laboratory
*   LPR/UFRGS - Residues Processing Laboratory
*		Federal University of Rio Grande do Sul
*		Porto Alegre (RS), Brazil
* $Id$
*--------------------------------------------------------------------*#

using "types";


#*---------------------------------------------------------------------
* Melhorias a fazer:
*	* adicionar o 'char' aos balanços de massa e energia;
*	* considerar a cinza nos balanços de massa e energia;
*	* adicionar cromo e cloro ao modelo;
*--------------------------------------------------------------------*#	



#*---------------------------------------------------------------------
*	Model of dry biomass
*--------------------------------------------------------------------*#	
Model dry_biomass
	ATTRIBUTES
	Brief = "Model of dry biomass";
	
	
	PARAMETERS
	NElem		as Integer 	(Brief="Number of elements", Default=5); # CHONS...
	Mw_(NElem)	as molweight(Brief="Molecular weight of elements");
	
	HHVcalc 	as Switcher (Brief="High heat value calculation", 
				Valid=["(Boie,1952)","(Zainal,2001)","(Higman,2003)","(Souza-Santos,2004)",
				"(Basu,2006)","Known data"], Default="(Souza-Santos,2004)");


	SET
	NElem = 5; # C,H,O,N,S
	Mw_ = [12.011,1.0079,15.99994,14.0067,32.06]*'kg/kmol'; # C,H,O,N,S

	
	VARIABLES
	massfrac(NElem)	as fraction	(Brief="Mass fraction (Ultimate analysis)", Symbol="frac_{mass}", Unit='kg/kg');
	molfrac(NElem)	as fraction	(Brief="Molar fraction", Symbol="frac_{mol}", Unit='kmol/kmol');
	ash			as fraction	(Brief="Mass ash fraction", Unit='kg/kg'); # only used to HHV calculation
	Mw			as molweight(Brief="Molecular weight");
	
	F			as flow_mol	(Brief="Molar flow rate");
	Fmass		as flow_mass(Brief="Mass flow rate", Symbol="F_{mass}");
	na(NElem)	as positive	(Brief="Matrix of elements");
	H			as enth_mol	(Brief="Molar enthalpy");
	Hmass		as enth_mass(Brief="Mass enthalpy", Upper=1e9, Symbol="H_{mass}");
	HHV			as enth_mass(Brief="Mass high heat value", Symbol="HHV_{mass}", Upper=1e10);
	HHV_		as enth_mass(Brief="Known data for high heat value", Upper=1e10);
	HHVmol		as enth_mol (Brief="Molar high heat value", Symbol="HHV_{mol}");
	LHV			as enth_mass(Brief="Mass low heat value", Upper=1e9, Symbol="LHV_{mass}");
	LHVmol		as enth_mol (Brief="Molar low heat value", Symbol="LHV_{mol}");
	
	
	EQUATIONS	
	"Molar fraction of fuel formula"
	na = massfrac*Mw_(1)/massfrac(1)/Mw_;
	
	"Mass fraction normalisation"
	sum(massfrac) = 1;

	"Molecular weight of fuel formula"
	Mw = sum(Mw_*na);
	
	"Molar fraction of biomass"
	molfrac = na/sum(na);	
	
	"Mass flow rate"
	Fmass = F*Mw;
	
	"Molar enthalpy"
	Hmass*Mw = H;
	
	switch HHVcalc # CHONSA
		case "(Boie,1952)":
		"Equation of Boie"
		HHV = sum([35.16,116.225,11.09,6.28,10.465]*massfrac)*'kJ/kg';
		
		case "(Zainal,2001)":
		"Reed and Levie (1985)"
		HHVmol = (sum([34.0945,13.23,-11.986,0,0]*massfrac/(1+ash)) - 1.53*ash + 6.85)*'kJ/kmol';
		
		case "(Higman,2003)":
		"Dulong formula"
		HHV = (sum([34.91,117.83,-10.34,-1.51,10.05]*massfrac/(1+ash)) - 2.11*ash)*'MJ/kg';
		
		case "(Souza-Santos,2004)":
		"Souza-Santos (2004)"
		HHV = (sum([34.245,110.198,-11.985,-11.985,0]*massfrac/(1+ash)) - 1.53*ash + 0.0685)*'MJ/kg';
		
		case "(Basu,2006)":
		"Dulong and Petit formula"
		HHV = sum([33.823,144.25,-14.28,0,9.418]*massfrac)*'MJ/kg';

		case "Known data":
		"Known data"
		HHV = HHV_;
	end
	
	
	"Molar high heat value"
	HHVmol = HHV*Mw;
	
	"Molar low heat value"
	LHVmol = LHV*Mw;
end



#*---------------------------------------------------------------------
*	Model of moisture biomass
*--------------------------------------------------------------------*#	
Model moisture_biomass
	ATTRIBUTES
	Brief = "Model of moisture biomass";
	
	
	PARAMETERS
#	NComp	as Integer 	(Brief="Number of components", Default=8);
	Mw		as molweight(Brief="Molecular weight");
	
	Ho		as enth_mol (Brief="Molar standard enthalpy of formation");
	Hv		as enth_mol (Brief="Molar enthalpy of vaporization");
	
#	Ho(NComp) 	as enth_mol (Brief="Molar standard enthalpy of formation");
#	Hv(NComp)	as enth_mol (Brief="Molar enthalpy of vaporization");
#	Tb(NComp)	as temperature (Brief="Boiling point temperature");
	
	
	SET
#	NComp = 8; #PP.NumberOfComponents;
	Mw = 18.0152*'kg/kmol'; # H2O
	
	Ho = -2.42e5*'kJ/kmol'; #PP.IdealGasEnthalpyOfFormationAt25C();
	Hv = 4.065e4*'kJ/kmol'; #PP.IdealGasEnthalpyOfFormation(Tb);
#	Tb = PP.NormalBoilingPoint();
	
	
	VARIABLES
	F			as flow_mol	(Brief="Molar flow rate");
	Fmass		as flow_mass(Brief="Mass flow rate", Symbol="F_{mass}");
	molfrac		as positive	(Brief="Molar relative moisture of biomass", Symbol="frac_{mol}", Unit='kmol/kmol');
	massfrac	as positive	(Brief="Mass relative moisture of biomass", Symbol="frac_{mass}", Unit='kg/kg', Default=0.1);
	H			as enth_mol	(Brief="Molar enthalpy");
	Hmass		as enth_mass(Brief="Mass enthalpy", Lower=-1e9, Symbol="H_{mass}");
	
	
	EQUATIONS
	"Mass flow rate"
	Fmass = F*Mw;
	
	"Mole enthalpy"
	H = Ho + Hv;	
	
	"Mass enthalpy"
	Hmass*Mw = H;
end



#*---------------------------------------------------------------------
*	Model of a biomass feed
*--------------------------------------------------------------------*#	
Model raw_biomass
	ATTRIBUTES
	Brief = "Model of a biomass feed";
	
	
	PARAMETERS
	To	as temperature	(Brief="Ambient temperature", Default=298);
	Po	as pressure		(Brief="Atmospheric pressure", Default=1);
	
	
	VARIABLES
	dry			as dry_biomass	(Brief="Dry biomass", Symbol="_{dry}");
	moisture	as moisture_biomass	(Brief="Moisture biomass", Symbol="_{moist}");
	
	F			as flow_mol		(Brief="Molar flow rate");
	Fmass		as flow_mass	(Brief="Mass flow rate", Symbol="F_{mass}");
	Mw			as molweight	(Brief="Molecular weight");
	
	h			as enth_mol 	(Brief="Molar enthalpy");
	hmass		as enth_mass	(Brief="Mass enthalpy", Upper=1e9, Symbol="H_{mass}");
	T			as temperature 	(Brief="Temperature");


	EQUATIONS
	"Molecular weight"
	Mw = (1 - moisture.molfrac)*dry.Mw + moisture.molfrac*moisture.Mw;
	
	"Molar fraction of moisture" # kmol of M by kmol of B
	moisture.molfrac*(moisture.massfrac+(1 - moisture.massfrac)*moisture.Mw/dry.Mw) = 
		moisture.massfrac;

	"Dry biomass flow rate"
	dry.Fmass = (1 - moisture.massfrac)*Fmass;
	
	"Moisture biomass flow rate"
	moisture.Fmass = moisture.massfrac*Fmass;

	"Raw biomass flow rate"
	Fmass = F*Mw;
	
	"Mass enthalpy"
	hmass = (1 - moisture.massfrac)*dry.Hmass + moisture.massfrac*moisture.Hmass;
	
	"Mole enthalpy"
	h = hmass*Mw;
	
	"Mass low heat value" ## adicionado de 'stream_feed'
	dry.LHV = dry.HHV - moisture.Hv*dry.massfrac(2)/2/dry.Mw_(2); # (Souza-Santos,2004)
#	dry.LHV = (dry.HHV*(1-moisture.massfrac)
#		- moisture.Hv*moisture.massfrac - (1-.moisture.massfrac))
#		*(18*dry.massfrac(2)/200); # (Mansaray,1998) p.26 -> ?????

#	dry.LHV = dry.HHV - 22604*'kJ/kg'*dry.massfrac(2) 
#		- 2581*'kJ/kg'*moisture.massfrac; # (Basu,2006) p.449	
	
	"Temperature"
	T = To;
end



#*---------------------------------------------------------------------
*	Model of an air stream
*--------------------------------------------------------------------*#	
Model air_stream
	ATTRIBUTES
	Brief = "Model of an air stream";
	
	
	PARAMETERS
outer PP 		as Plugin	(Brief="External physical properties", Type="PP");
	NComp		as Integer 	(Brief="Number of components", Default=8);
	Mw(NComp)	as molweight;
	
	To	as temperature	(Brief="Ambient temperature", Default=298);
	Po	as pressure		(Brief="Atmospheric pressure", Default=1);
	
	
	VARIABLES
	F			as flow_mol	(Brief="Molar flow rate");
	Fmass		as flow_mass(Brief="Mass flow rate", Symbol="F_{mass}");
	Fvol		as flow_vol	(Brief="Volumetric flow rate", Symbol="F_{vol}");
	
	z(NComp)	as fraction	(Brief="Molar fraction", Unit='kmol/kmol');
	zmass(NComp)as fraction	(Brief="Mass fraction", Unit='kg/kg', Symbol="z_{mass}");

	sh			as positive	(Brief="Specific humidity", Unit='g/kg');
	rh			as percent 	(Brief="Relative humidity");

	Mws			as molweight;
	vm			as volume_mol;
	T			as temperature;
	P			as pressure;
	Pv(NComp)	as pressure;
	h			as enth_mol;
	
	
	SET
	NComp = 8; # PP.NumberOfComponents;
	Mw = PP.MolecularWeight();
	
	
	EQUATIONS
	"Mixture Molecular weight"
	Mws = sum(z*Mw);
	
	"Mass flow rate"
	Fmass = F*Mws;
	
	"Volumetric flow rate"
	Fvol = F*vm;
	
	"Temperature"
	T = To;
	
	"Pressure"
	P = Po;
	

	"Molar Volume"
	vm = PP.VapourVolume(T,P,z);
	
	
	"Mass water fraction"
	zmass(1) = sh;
	
	"Mass oxygen fraction"
	zmass(2) = 0.2316*(1-zmass(1)); # (Basu,2006) p.446: 23.16% O2, 76.8% N2, 0.04% inerts
	
	"Mass nitrogen fraction"
	sum(zmass) = 1;
	
	"Air composition" # <<---
	z(4:NComp) = 0;
#	zmass(4:8) = 0;
	
	"Molar fraction"
	z*sum(zmass/Mw) = zmass/Mw;

	"Relative humidity"
	rh = z(1)*P/Pv(1)*100;
	
	"Vapour pressure"
	Pv = PP.VapourPressure(T);
	
	"Enthalpy"
	h = sum(PP.IdealGasEnthalpyOfFormationAt25C()*z) + PP.VapourEnthalpy(T,P,z);
end



Model stream_gasifier
	ATTRIBUTES
	Brief = "Model of a gasifier stream";
	
	
	PARAMETERS
outer PP 		as Plugin	(Brief="External physical properties", Type="PP");
	NComp		as Integer 	(Brief="Number of components", Default=8);
	NElem		as Integer 	(Brief="Number of elements", Default=5);
	na(NElem,NComp)as positive(Brief="Matrix of elements per component");
	Mw(NComp)	as molweight(Brief="Molecular weight of components");
	Ho(NComp) 	as enth_mol	(Brief="Molar component enthalpy");
	
	
	VARIABLES
	F 			as flow_mol (Brief="Molar Flow Rate");
	Fmass		as flow_mass(Brief="Mass flow rate", Symbol="F_{mass}");
	z(NComp) 	as fraction	(Brief="Molar Fraction");
	zmass(NComp)as fraction	(Brief="Mass fraction", Unit='kg/kg', Symbol="z_{mass}");
	
	N(NComp)	as positive	(Brief="Mole fraction of component by initial biomass", Unit='kmol/kmol', Symbol="N_{mol}");
	Nmass(NComp)as positive	(Brief="Mass fraction of component by initial biomass", Unit='kmol/kmol', Symbol="N_{mass}");

	Mws			as molweight(Brief="Molecular weight of stream");

	T 			as temperature(Brief="Temperature");
	P 			as pressure (Brief="Pressure");
	
	h 			as enth_mol (Brief="Molar stream enthalpy");
	hmass		as enth_mass(Brief="Mass stream enthalpy", Lower=-1e10, Upper=1e10, Symbol="h_{mass}");
	E			as enth_mol	(Brief="Specific energy");
	v 			as fraction	(Brief="Vapourization fraction");
	

	SET
	NComp = 8; # PP.NumberOfComponents;
	NElem = 5; # C,H,O,N,S
	
	#		   C H O N S
	na(:,1) = [0,2,1,0,0]; # H2O
	na(:,2) = [0,0,2,0,0]; # O2
	na(:,3) = [0,0,0,2,0]; # N2
	na(:,4) = [1,0,1,0,0]; # CO
	na(:,5) = [1,0,2,0,0]; # CO2
	na(:,6) = [1,4,0,0,0]; # CH4
	na(:,7) = [0,2,0,0,0]; # H2
	na(:,8) = [0,0,2,0,1]; # SO2
	
	Mw = PP.MolecularWeight();
	Ho = PP.IdealGasEnthalpyOfFormationAt25C();

	EQUATIONS	
	"Mass stream enthalpy"
	hmass*Mws = h;	

	"Mole fraction normalisation"
	sum(z) = 1;
	
	"Vapour stream"
	v = 1;
end



#*---------------------------------------------------------------------
*	Model of a specific stream
*--------------------------------------------------------------------*#	
Model stream_feed as stream_gasifier
	ATTRIBUTES
	Brief = "Model of a specific feed stream";
	
	
	PARAMETERS
	To	as temperature	(Brief="Ambient temperature", Default=298);
	Po	as pressure		(Brief="Atmospheric pressure", Default=1);
	
	
	VARIABLES
	Fuel	as raw_biomass	(Brief="Raw biomass", Symbol="_{fuel}");
	Air		as air_stream	(Brief="Air stream", Symbol="_{air}");
	
	phi		as Real			(Brief="Equivalence ratio", Lower=0, Symbol="\phi");
	Fw		as Real			(Brief="Mass air-fuel ratio", Lower=0);
	Fm		as Real			(Brief="Molar air-fuel ratio", Lower=0);
	
	
	EQUATIONS
	"Molar flow rate"
	F = Fuel.moisture.F + Air.F;

	"Mass stream feed flow rate"
	Fmass = Fuel.moisture.Fmass + Air.Fmass;
	
	
	"Molecular weight of stream"
	Mws*F = Fuel.moisture.F*Fuel.moisture.Mw + Air.F*Air.Mws;

	"Mass fraction"
	zmass*Mws = z*Mw;
	
	"Molar fraction by biomass"
	N*Fuel.dry.F = z*F;
	
	"Mass fraction by biomass"
	Nmass*Fuel.dry.Mw = N*Mw;

	
	"Molar stream enthalpy"
	h*F = Fuel.moisture.F*Fuel.moisture.H + Air.F*Air.h;
	
#*---------------------------------------------------------------------
*	C  + O2 -> CO2						Hc(C) < 0
*	H2 + O2 -> H2O						Hc(H) < 0
*	S  + O2 -> SO2						Hc(S) < 0
*	CO2 + H2O + N2 + SO2 -> CHNOS + O2	LHV > 0
*	C + H2 + N2 + O2 + S -> CHNOS		Hf = LHV + Hc(C+H+S)
*--------------------------------------------------------------------*#	
	"Biomass enthalpy"
	Fuel.dry.H = Fuel.dry.LHVmol 
		+ (Fuel.dry.na(1)*Ho(5) + 0.5*Fuel.dry.na(2)*Ho(1) 
		+ Fuel.dry.na(5)*Ho(8)); # CHONS CO2,H2O,SO2
# 	H2O, O2, N2, CO, CO2, CH4, H2, SO2
	
	"Input energy"
	E = Fuel.dry.LHVmol;
	
	"Temperature"
	T = To;
	"Pressure"
	P = Po;
	
	
	# [H2O, O2, N2]/[Initial dry biomass]
	
	"Mass water fraction"
	Nmass(1) = Air.sh*Air.Fmass/Fuel.dry.Fmass + Fuel.moisture.massfrac;
	
	"Mass oxygen fraction"
	Nmass(2) = Air.zmass(2)*Air.Fmass/Fuel.dry.Fmass;
	
	"Molar fraction" # only O2,N2,H2O in the air
	Nmass(4:NComp) = 0;
	

	"Equivalence ratio"
	phi = N(2)/(1 + 0.25*Fuel.dry.na(2) - 0.5*Fuel.dry.na(3) + Fuel.dry.na(5));
	
	"Equivalence rate of gasifying" # (Melgar,2007)
	Fm = Air.F/Fuel.dry.F; # molar base
	Fw = Air.Fmass/Fuel.dry.Fmass; # mass base
end



#*---------------------------------------------------------------------
*	Model of a product stream
*--------------------------------------------------------------------*#	
Model stream_products as stream_gasifier
	ATTRIBUTES
	Brief = "Model of a product stream";

	
	PARAMETERS
	LHV_syngas(3)	as enth_mass;
	
	
	VARIABLES
	Fvol	as flow_vol	(Brief="Volumetric flow rate (no solid)", Symbol="F_{vol}");
	vm		as volume_mol(Brief="Molar Volume");
	z_db(NComp-1)as fraction	(Brief="Molar fraction in dry base");
	
	LHV		as enth_mass(Brief="Mass low heat value", Upper=1e10, Symbol="LHV_{mass}");


	SET
	LHV_syngas = [12.622,35.814,10.788]*'MJ/kg'; # CO,CH4,H2
	
	
	EQUATIONS
	"Volumetric flow rate (no solid)"
	Fvol = F*vm;

	"Mass flow rate"
	Fmass = F*Mws;

	"Molar Volume"
	vm = PP.VapourVolume(T,P,z);
	
	"Molecular weight of stream"
	Mws = sum(z*Mw);
	
	"Mass fraction"
	zmass*Mws = z*Mw;

	"Mass fraction in dry base"
	z_db*(1 - z(1)) = z(2:NComp);
	
	
	"Low heat value of biomass syn-gas"
	LHV = (zmass(4)*LHV_syngas(1) + zmass(6)*LHV_syngas(2)
		+ zmass(7)*LHV_syngas(3)); # /(zmass(4) + zmass(6) + zmass(7));
end



#######################################################################
#							MAIN MODEL
#######################################################################

#*---------------------------------------------------------------------
*	Equilibrium model of a biomass gasifier
*--------------------------------------------------------------------*#
Model gasifier
	ATTRIBUTES
	Info = "
== Description ==
	Thermodynamic equilibrium modeling of a biomass gasifier.

== Assumptions ==
* thermodynamic equilibrium;
* steady-state;
* ideal gas relations;
* ambient and atmospheric input conditions (To,Po);
* global gasification reaction:
		CHONS + H2O + air -> CO + CO2 + CH4 + H2 + H2O + N2 + SO2
* no oxygen in the products.

== Specify ==
* ultimate biomass analysis (dry.massfrac);
* moisture biomass ratio (moisture.massfrac);
* equivalence ratio (phi);
* relative air humidity (air.rh);
* equilibrium temperature (Teq).
";	


	PARAMETERS
	PP	as Plugin(Brief="External physical properties",
		Type="PP", 
		Components = ["water","oxygen","nitrogen","carbon monoxide",
		"carbon dioxide","methane","hydrogen","sulfur dioxide"],
		LiquidModel = "IdealLiquid", 
		VapourModel = "Ideal"
	);
	
	NComp	as Integer 	(Brief="Number of components", Default=8);
	NElem	as Integer 	(Brief="Number of elements", Default=5);
	
	NReac 	as Integer 	(Brief="Number of reactions", Default=2);
	stoic(NComp,NReac) 	as Real (Brief="Stoichiometric matrix of considered reaction", Symbol="\nu");
	Rg		as Real		(Brief="Universal gas constant", Unit='J/mol/K', Default=8.314);	
	Po		as Integer	(Brief="Atmospheric pressure", Default=1);

	VARIABLES
	Inlet		as stream_feed	(Brief="Inlet stream", Symbol="^{in}", PosX=0.7165, PosY=0);
	Outlet		as stream_products(Brief="Outlet stream", Symbol="^{out}", PosX=0.7165, PosY=1);
	
	K(NReac)	as positive		(Brief="Equilibrium constant", Lower=0);
	G(NComp)	as enth_mol		(Brief="Component Gibbs free energy");	
	Teq			as temperature	(Brief="Equilibrium temperature", Symbol="T_{eq}");
	
	eff_cg		as Real	 		(Brief="Cold gas efficiency", Symbol="\eta_{cg}");
	
	
	SET
	NComp = PP.NumberOfComponents; # 8
	NElem = 5; # C,H,O,N,S
#*---------------------------------------------------------------------
*	Considered equilibrium reaction:
* 		C + 2H2 <-> CH4
*		CO + H2O <-> CO2 + H2
*--------------------------------------------------------------------*#	
	NReac = 2;
	stoic(:,1) = [ 0, 0, 0, 0, 0, 1,-2, 0]; # C + 2H2 <-> CH4
	stoic(:,2) = [-1, 0, 0,-1, 1, 0, 1, 0]; # CO + H2O <-> CO2 + H2
	
	
	EQUATIONS
	for k in [1:NElem]
	"Material balance"
		Inlet.F*sum(Inlet.z*Inlet.na(k,:))
			+ Inlet.Fuel.dry.Fmass*Inlet.Fuel.dry.massfrac(k)/Inlet.Fuel.dry.Mw_(k) 
			= Outlet.F*sum(Outlet.z*Outlet.na(k,:));
	end


	"Fraction normalisation"
	sum(Outlet.z) = 1;


	"Heat balance"
	Inlet.F*(Inlet.z(1)*Inlet.Fuel.moisture.H + sum(Inlet.z(2:NComp)*Inlet.Ho(2:NComp))) 
		+ Inlet.Fuel.dry.F*Inlet.Fuel.dry.H = Outlet.F*Outlet.h;
	
	
	for j in [1:NReac]
	"Equilibrium constant"
		K(j) = prod(Outlet.z^stoic(:,j));

	"Gibbs energy of reaction"
		sum(stoic(:,j)*G) = -Rg*Outlet.T*ln(K(j));
	end
	
	
	"Gibbs energy of component"
	G = PP.IdealGasGibbsOfFormation(Outlet.T);

	Outlet.h = sum(PP.IdealGasEnthalpyOfFormationAt25C()*Outlet.z)
		+ PP.VapourEnthalpy(Outlet.T,Outlet.P,Outlet.z);


	"Output energy"
	Outlet.E = Outlet.N(4)*(Inlet.Ho(4)-Inlet.Ho(5)) + Outlet.N(7)*(Inlet.Ho(7)
		-Inlet.Ho(1)) + Outlet.N(6)*(Inlet.Ho(6)-Inlet.Ho(5)-Inlet.Ho(1));
	
	"Cold gas efficiency"
	eff_cg = Outlet.E/Inlet.E;
	
	
	# Outlet stream
	
	"Equilibrium temperature"
	Outlet.T = Teq;

	"Mechanical equilibrium"
	Outlet.P = Inlet.P;
	
	
	"Molar fraction by initial biomass"
	Outlet.N*Inlet.Fuel.dry.F = Outlet.z*Outlet.F;
	
	"Mass fraction by initial biomass"
	Outlet.Nmass*Inlet.Fuel.dry.Mw = Outlet.N*Outlet.Mw;
end