% 10.34 - Fall 2006 % Hw Set #7 - Problem #3 % Rob Ashcraft - Oct. 20, 2006 %equilibrium solving function p3_main clear; clc; moles_spec_in = [2; 3; .5; 0; 0; 1; 0]; options = optimset('Display','iter','TolFun',1e-8,'TolX',1e-14,'TolCon',1e-14); % Define inequality constraints (A*x <= b) A = eye(7)*(-1); b = ones(7,1)*1e-14; % define the C, H, and O conservation (equality constraints) Aeq = [1 0 1 1 0 0 0; 4 2 0 0 1 2 2; 0 1 1 2 2 0 2]; moles_atoms_in = Aeq*moles_spec_in; C_in = moles_atoms_in(1); %moles H_in = moles_atoms_in(2); O_in = moles_atoms_in(3); beq = [C_in; H_in; O_in]; lb = ones(7,1)*1e-14; ub = ones(7,1)*20; n_var0 = [ 0.0459 1.9934 1.2958 0.6583 0.0000 5.7176 0.1042]; [n_var, fval] = fmincon(@obj_fun,n_var0,A,b,Aeq,beq,lb,ub,[],options); n_total = sum(n_var); n_ch4 = n_var(1); n_h2o = n_var(2); n_co = n_var(3); n_co2 = n_var(4); n_ho2 = n_var(5); n_h2 = n_var(6); n_hooh = n_var(7); y_ch4 = n_ch4/n_total; y_h2o = n_h2o/n_total; y_co = n_co/n_total; y_co2 = n_co2/n_total; y_ho2 = n_ho2/n_total; y_h2 = n_h2/n_total; y_hooh = n_hooh/n_total; disp(' '); disp(' '); disp(['Total moles initially: ',num2str(sum(moles_spec_in))]); disp(['Initial Volume (L): ',num2str(sum(moles_spec_in)*0.08206*1000/1)]); disp(' '); disp(['Total moles at equilibrium: ',num2str(n_total)]); disp(['Final Volume (L): ',num2str(n_total*0.08206*1000/1)]); disp(' '); disp(['CH4: mole fraction: ',num2str(y_ch4,'%1.5f'),', number of moles: ',num2str(n_ch4,'%1.5f')]); disp(['H2O: mole fraction: ',num2str(y_h2o,'%1.5f'),', number of moles: ',num2str(n_h2o,'%1.5f')]); disp(['CO: mole fraction: ',num2str(y_co,'%1.5f'),', number of moles: ',num2str(n_co,'%1.5f')]); disp(['CO2: mole fraction: ',num2str(y_co2,'%1.5f'),', number of moles: ',num2str(n_co2,'%1.5f')]); disp(['H2: mole fraction: ',num2str(y_h2,'%1.5f'),', number of moles: ',num2str(n_h2,'%1.5f')]); disp(['HOOH: mole fraction: ',num2str(y_hooh,'%1.5f'),', number of moles: ',num2str(n_hooh,'%1.5f')]); disp(['HO2: mole fraction: ',num2str(y_ho2,'%1.2e'),', number of moles: ',num2str(n_ho2,'%1.2e')]); function cost = obj_fun(n_var) T = 1000; %kelvin R = 8.314; % J/mol-K n_total = sum(n_var); n_ch4 = n_var(1); n_h2o = n_var(2); n_co = n_var(3); n_co2 = n_var(4); n_ho2 = n_var(5); n_h2 = n_var(6); n_hooh = n_var(7); %Gibb free energies at 1000K in kJ/mole G_ch4 = 19720; G_h2o = -192420; G_co = -200240; G_co2 = -395790; G_ho2 = -227000; G_h2 = 0; G_hooh = -369060; % minimize the chemical potential of the system: cost = n_ch4*(G_ch4 + R*T*log(n_ch4/n_total)) +... n_h2o*(G_h2o + R*T*log(n_h2o/n_total)) +... n_co*(G_co + R*T*log(n_co/n_total)) +... n_co2*(G_co2 + R*T*log(n_co2/n_total)) +... n_ho2*(G_ho2 + R*T*log(n_ho2/n_total)) +... n_h2*(G_h2 + R*T*log(n_h2/n_total)) +... n_hooh*(G_hooh + R*T*log(n_hooh/n_total)); % cost = n_ch4*(G_ch4) +... % n_h2o*(G_h2o) +... % n_co*(G_co) +... % n_co2*(G_co2) +... % n_ho2*(G_ho2) +... % n_h2*(G_h2); return;