% Computes the equilibrium composition of a reacting mixture global ng ns nr tempin prin ctot xg nu global keq298 keq delhr global n ntot global rgas %number of species and number of reactions ng = 4 ; nr = 2 ; % initialize the vectors. keq298 = zeros(nr, 1 ); keq = zeros(nr, 1 ); delhr = zeros(nr, 1 ); xg = zeros(ng,1); nu = zeros ( nr, ng); zeta = zeros(nr,1) ; % provide values for stoichiometry and eq constants. nu(1,1) = -1.0 ; nu(1,2) = -1.0; nu(1,3) = 2.0 ; nu(1,4)= 0.0; nu(2,1) = 0.0 ; nu(2,2) = -0.5; nu(2,3) = -1.0 ; nu(2,4)= 1.0; keq298(1) = 4.069E-31; delhr(1) = 180500.0; keq298(2) = 1.22E+06 ; delhr(2) = -57070.0; %assumed constant (not a function of T) rgas = 8.314472 ; % gas constant; use consistent units % provide the feed conditions tempin = 1800.0; prin = 1.0e+05; % initial moles xg(1) = 96.7; xg(2) = 3.3; xg(3) = 0.0; xg(4) =0.0; % provide initial values for extent of reactions zeta(1) = 0.4; zeta(2) = 0.1; zeta = fsolve ('fun_eq', zeta) ---------------------------------------------------- function fvec = fun_eq( zeta) % Defines the (nr) functions to be solved global ng ns nr tempin prin ctot xg nu global keq298 keq delhr global n ntot global rgas fvec = zeros (nr,1) ; ptotatm = prin /1.0e+05 ; % find K at the desired conditions for i = 1: nr keq(i)=keq298(i)*exp( delhr(i) /rgas *(1./298-1./tempin)) ; end %number of moles of each species for the given extents % nj0 = 1 for these calculations for j = 1: ng n(j) = xg(j); for i=1:nr n(j) = n(j) + nu(i,j) * zeta(i); end end % find ntot ntot = 0.0; for j = 1: ng ntot = ntot + n(j); end % find the mole fractions for each species pp = n /ntot * ptotatm; % set up discrepancy function for each reaction for i = 1: nr prod = 1.0; for j = 1: ng prod = prod * pp(j)^nu(i,j); end fvec(i) = keq(i) - prod end