HGS. Chemical Equilibrium Calculation

by

 

19 Mar 2013 (Updated )

Computation of chemical equilibrium problems related to combustion and isentropic expansion.

Ex07_isentropic_comparison_rpa
%***********************************************************************************************************
%* HGS 1.1 
%* By Arnau Miro, Pau Manent, Eva Astrid Jara
%* Supervised by Manel Soria and Ramon Carreras
%* LLOP, ETSEIAT UPC          
%***********************************************************************************************************
%
% Example 07: isentropic comparison with RPA software

function Ex07_isentropic_comparison_rpa

clear; clc;

shifting = true % conditions of shifting flow are assumed

species={'H','H2','H2O','H2O2','HO2','O','O2','OH'}

% Nozzle inlet mole fractions (RPA)
ni_rpa=[ 0.0644043;...  % H
    0.1402066;...       % H2
    0.6176198;...       % H2O
    0.0000024;...       % H2O2
    0.0000367;...       % HO2
    0.0268578;...       % O
    0.0467048;...       % O2
    0.1041674];         % OH
    
Tc=3027.58  % K (RPA)
Pc=1        % bar (RPA)

% 1-Verification of properties at given Tc and Pc
disp('1-Verification of properties at given Tc and Pc');

[Cp,Cv,MM,Rg,gamma,a,H,~,S]=hgsprop(species,ni_rpa,Tc,Pc);
n=sum(ni_rpa);  % total number of mols in the mixture (1)
m=n*MM*1e-3;    % mixture total mass in kg

Pc
Tc
h=H/m   % kJ/kg
s=S/m   % kJ/kgK
cp=Cp/m % kJ/kgK
cv=Cv/m
Rg
MM
gamma
rho=Pc/(Rg*Tc)
a

% Results are more or less ok except cp and cv which are quite different.
% However if a delta t of 1 K, the difference of h is similar to cp.

% 2-Verification of the equilibrium composition given T and P
disp('Verification of the equilibrium composition given T and P');

ni_calc=hgseq(species,ni_rpa,Tc,Pc)

error_tantpercent=100*(ni_calc-ni_rpa) / norm(ni_rpa)

% 3-Verification of the insentropic expansion at a certain pressure Pt
disp('3-Verification of the insentropic expansion at a certain pressure Pt');

Pt=0.1 
    function DeltaS=DeltaS(T)
        if shifting
            nt=hgseq(species,ni_rpa,T,Pt); % Shifting, equilibrium conditions at T,Pe
        else
            nt=ni_rpa; % frozen
        end
        [~,~,MM2,~,~,~,~,~,S2]=hgsprop(species,nt,T,Pt); % enthropy at T,Pt
        n2=sum(nt); % mixture total number of mols (1)
        m2=n2*MM2*1e-3; % mixture mass in kg
        s2=S2/m2; % kJ/kgK
        DeltaS=s2-s;
    end

Tt=fzero(@DeltaS,3000,optimset('Display','iter'))
if shifting
    nt=hgseq(species,ni_rpa,Tt,Pt) 
else
    nt=ni_rpa
end

[~,~,MM2,~,~,a2,H2,G2,S2]=hgsprop(species,nt,Tt,Pt);
n2=sum(nt); % mitxture total number of mols (1)
m2=n2*MM2*1e-3; % mixtur mass in kg
s
s2=S2/m2 % kJ/kgK
S
S2

h2=H2/m2

vt=sqrt(2*1000*(h-h2)) % Enthalpy en J !

mt=vt/a2 % Mach 

Is=vt/9.81 % Is (optimal expansion, Pe=Pambient)

end

Contact us