% 10.34 - Fall 2006 % HW Set#3, Problem #1 % Rob Ashcraft - Sept. 19, 2006 % Calculate the peak temperature for a water-gas shift reactor. This % function calculates the peak temperature for complete reaction and for a % reaction that proceeds to equilibrium. % CO + H2O <=> CO2 + H2 function pset3p1_main clear; clc; global CH4_coef O2_coef H2_coef H2O_coef CO_coef CO2_coef R Rj V % Get the Shomate coefficients for the species of interest [CH4_coef, O2_coef, H2_coef, H2O_coef, CO_coef, CO2_coef] = shomate_coef_wgs; R = 0.08206; % gas constant in atm-L/mole-K Rj = 8.314; % gas constant in J/mol-K V = 5; % reactor volume (not needed in this problem) T_in = [500:25:1500]'; % vector of Temps to evaluate things at % Inlet mole flow rates to the reactor (moles/sec) N_in = [0, 0, 2, 5, 1, 0.1]*1e-3; CH4_in = N_in(1); O2_in = N_in(2); H2_in = N_in(3); H2O_in = N_in(4); CO_in = N_in(5); CO2_in = N_in(6); % Calculate the thermo at each temperature of interest for i=1:length(T_in) [dH(i),TdS(i),dG(i)] = wgs_thermo(T_in(i)); end % This section will solve for the equilibrium conditions in the reactor at % each temperature of interest. for i=1:length(T_in) Tout_0 = T_in(i); % guess for Tout z1_0 = N_in(5)*0.9; % guess for extent of reaction var_0 = [Tout_0; z1_0]; options = optimset('Display','off','MaxFunEvals',10000,'MaxIter',1000,'TolX',1e-8,'TolFun',1e-8); % use fsolve to solve the problem var = fsolve(@wgs_solver,var_0,options,T_in(i),N_in); Tout(i,1) = var(1); z1(i,1) = var(2); conv_eq(i,1) = 100*(1 - (CO_in - z1(i))/CO_in); % calculate conversion from the extent end CH4 = CH4_in; O2 = O2_in; H2 = H2_in + z1; H2O = H2O_in - z1; CO = CO_in - z1; CO2 = CO2_in + z1; % This part solves for the outlet temperature for various combinations of % conversion and inlet temperature options2 = optimset('Display','off','MaxFunEvals',10000,'MaxIter',1000,'TolX',1e-8,'TolFun',1e-8); conv2 = linspace(0,100,21); for i=1:length(conv2) for j=1:length(T_in) T_mat(i,j) = T_in(j); %generate a matrix of temperatures used conv_mat(i,j) = conv2(i); %generate a matrix of conversions used z(i) = CO_in*conv2(i)/100; % convert the conversion to extent % solve for the outlet temperatures Tout_z(i,j) = fsolve(@wgs_solver_conv,T_in(j)*0.75,options2,T_in(j),N_in,z(i)); % Since we are asked to plot only the feasible region, we need to % compare the conversion at each point here to the equilibrium % conversion calculated above for each inlet temperature. If the % point is infeasible, then the outlet temp value is replaced with % NaN, which will not show up in a plot. if(conv_mat(i,j) > conv_eq(j)) Tout_z_constrained(i,j) = NaN; else Tout_z_constrained(i,j) = Tout_z(i,j); end end end % pull out the 100% conversion temperatures, as this was asked for in the % problem as well (the last row in Tout_z) Tout_100 = Tout_z(end,:)'; % generate the figures required by the problem figure; plot(T_in,dH,'-',T_in,TdS,'--',T_in,dG,'-.'); title('Thermochemistry for the Water-Gas Shift Reaction'); xlabel('Temperature [K]'); ylabel('Thermochemistry [kJ/mole]'); legend('dH(Rxn)','T*dS(Rxn)','dG(Rxn)'); set(gca,'xtick',[500:100:1500]); figure; set(get(0,'CurrentFigure'),'position',[275 594 1344 504]); subplot(1,2,1); plot(T_in,Tout_100 - T_in); title('Complete Conversion Temperatures'); xlabel('Inlet Temperature [K]'); ylabel('T_O_u_t - T_I_n [K]'); axis([500 1500 0 150]); set(gca,'xtick',[500:100:1500]); subplot(1,2,2); plot(T_in,Tout-T_in); title('Equilibrium Conversion Temperatures'); xlabel('Inlet Temperature [K]'); ylabel('T_O_u_t - T_I_n [K]'); axis([500 1500 0 150]); set(gca,'xtick',[500:100:1500]); figure; set(get(0,'CurrentFigure'),'position',[275 594 1344 504]); subplot(1,2,1); plot(T_in,Tout_100); title('Complete Conversion Temperatures'); xlabel('Inlet Temperature [K]'); ylabel('T_O_u_t [K]'); set(gca,'xtick',[500:100:1500]); subplot(1,2,2); plot(T_in,Tout); title('Equilibrium Conversion Temperatures'); xlabel('Inlet Temperature [K]'); ylabel('T_O_u_t [K]'); set(gca,'xtick',[500:100:1500]); figure; plot(T_in,conv_eq); title('Equilibrium Conversion'); xlabel('Inlet Temperature [K]'); ylabel('CO Conversion [%]'); set(gca,'xtick',[500:100:1500]); figure; surf(conv_mat,T_mat,Tout_z_constrained,'facecolor','interp'); %'EdgeColor','none', colormap 'jet'; colorbar; axis([0 1 0 1 500 1600]); axis 'auto xy'; title('Physically Realizable Reactor Conditions'); ylabel('T_I_n [K]'); xlabel('CO Conversion [%]'); zlabel('T_O_u_t [K]'); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Other functions used in this problem: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function [dH,TdS,dG] = wgs_thermo(T_in) % this function return the thermochemistry for the species as at the % requested temperature % input: T_in = temperature at which thermo will be calculated % output: dH: enthaply of reaction in kJ/mole % dG: Gibb energy of reaction in kJ/mole % TdS: Temp times the entropy of reaction in kJ/mole global CH4_coef O2_coef H2_coef H2O_coef CO_coef CO2_coef R Rj V [Cp_CH4, H_CH4, S_CH4] = Shomate(T_in,CH4_coef); [Cp_O2, H_O2, S_O2] = Shomate(T_in,O2_coef); [Cp_H2, H_H2, S_H2] = Shomate(T_in,H2_coef); [Cp_H2O, H_H2O, S_H2O] = Shomate(T_in,H2O_coef); [Cp_CO, H_CO, S_CO] = Shomate(T_in,CO_coef); [Cp_CO2, H_CO2, S_CO2] = Shomate(T_in,CO2_coef); G_CH4 = H_CH4 - T_in*S_CH4/1000; G_O2 = H_O2 - T_in*S_O2/1000; G_H2 = H_H2 - T_in*S_H2/1000; G_H2O = H_H2O - T_in*S_H2O/1000; G_CO = H_CO - T_in*S_CO/1000; G_CO2 = H_CO2 - T_in*S_CO2/1000; dH = H_CO2 + H_H2 - (H_CO + H_H2O); %kJ/mol TdS = T_in*(S_CO2 + S_H2 - (S_CO + S_H2O))/1000; %kJ/mol dG = G_CO2 + G_H2 - (G_CO + G_H2O); %kJ/mol return %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function resid = wgs_solver(var,T_in,N_in) % this is used to find the equilibrium conditions coming out of the reactor % for a given inlet temp and composition % inputs: % var: [Temp out (K), extent of reaction (moles)] % T_in: inlet temperature to reactor (K) % N_in: inlet mole flow rates (moles/sec) % outputs: % resid: residuals of the equations we are trying to solve global CH4_coef O2_coef H2_coef H2O_coef CO_coef CO2_coef R Rj V CH4_in = N_in(1); O2_in = N_in(2); H2_in = N_in(3); H2O_in = N_in(4); CO_in = N_in(5); CO2_in = N_in(6); Tout = var(1); z1 = var(2); % get the thermo parameters for inlet and outlet temps [Cp_CH4_in,H_CH4_in,S_CH4_in] = Shomate(T_in,CH4_coef); [Cp_O2_in,H_O2_in,S_O2_in] = Shomate(T_in,O2_coef); [Cp_H2_in,H_H2_in,S_H2_in] = Shomate(T_in,H2_coef); [Cp_H2O_in,H_H2O_in,S_H2O_in] = Shomate(T_in,H2O_coef); [Cp_CO_in,H_CO_in,S_CO_in] = Shomate(T_in,CO_coef); [Cp_CO2_in,H_CO2_in,S_CO2_in] = Shomate(T_in,CO2_coef); [Cp_CH4,H_CH4,S_CH4] = Shomate(Tout,CH4_coef); [Cp_O2,H_O2,S_O2] = Shomate(Tout,O2_coef); [Cp_H2,H_H2,S_H2] = Shomate(Tout,H2_coef); [Cp_H2O,H_H2O,S_H2O] = Shomate(Tout,H2O_coef); [Cp_CO,H_CO,S_CO] = Shomate(Tout,CO_coef); [Cp_CO2,H_CO2,S_CO2] = Shomate(Tout,CO2_coef); % calculate the Gibbs for each species in J/mole G_CH4 = 1000*H_CH4 - Tout*S_CH4; G_O2 = 1000*H_O2 - Tout*S_O2; G_H2 = 1000*H_H2 - Tout*S_H2; G_H2O = 1000*H_H2O - Tout*S_H2O; G_CO = 1000*H_CO - Tout*S_CO; G_CO2 = 1000*H_CO2 - Tout*S_CO2; % Define the molar flow rate out in terms of extents of reaction: CH4 = CH4_in; O2 = O2_in; H2 = H2_in + z1; H2O = H2O_in - z1; CO = CO_in - z1; CO2 = CO2_in + z1; % Enthalpy balance resid(1,1) = CH4_in*H_CH4_in + O2_in*H_O2_in + H2_in*H_H2_in +... H2O_in*H_H2O_in + CO_in*H_CO_in + CO2_in*H_CO2_in -... (CH4*H_CH4 + O2*H_O2 + H2*H_H2 + H2O*H_H2O + CO*H_CO + CO2*H_CO2); % equilibrium for Reaction dG = G_CO2 + G_H2 - (G_CO + G_H2O); dn = 0; resid(2,1) = (1/R/Tout)^dn*exp(-dG/Rj/Tout) - CO2*H2/(CO*H2O)*(1/V)^dn; return %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function resid = wgs_solver_conv(Tout,T_in,N_in,z1) % this is used to find the temperature coming out of the reactor % for a given conversion, inlet temp, and inlet composition % inputs: % Tout: Temp out (K) % T_in: inlet temperature to reactor (K) % N_in: inlet mole flow rates (moles/sec) % z1: given extent of reaction (moles) % outputs: % resid: residuals of the equations we are trying to solve global CH4_coef O2_coef H2_coef H2O_coef CO_coef CO2_coef R Rj V CH4_in = N_in(1); O2_in = N_in(2); H2_in = N_in(3); H2O_in = N_in(4); CO_in = N_in(5); CO2_in = N_in(6); [Cp_CH4_in_in,H_CH4_in,S_CH4_in] = Shomate(T_in,CH4_coef); [Cp_O2_in,H_O2_in,S_O2_in] = Shomate(T_in,O2_coef); [Cp_H2_in,H_H2_in,S_H2_in] = Shomate(T_in,H2_coef); [Cp_H2O_in,H_H2O_in,S_H2O_in] = Shomate(T_in,H2O_coef); [Cp_CO_in,H_CO_in,S_CO_in] = Shomate(T_in,CO_coef); [Cp_CO2_in,H_CO2_in,S_CO2_in] = Shomate(T_in,CO2_coef); [Cp_CH4,H_CH4,S_CH4] = Shomate(Tout,CH4_coef); [Cp_O2,H_O2,S_O2] = Shomate(Tout,O2_coef); [Cp_H2,H_H2,S_H2] = Shomate(Tout,H2_coef); [Cp_H2O,H_H2O,S_H2O] = Shomate(Tout,H2O_coef); [Cp_CO,H_CO,S_CO] = Shomate(Tout,CO_coef); [Cp_CO2,H_CO2,S_CO2] = Shomate(Tout,CO2_coef); CH4 = CH4_in; O2 = O2_in; H2 = H2_in + z1; H2O = H2O_in - z1; CO = CO_in - z1; CO2 = CO2_in + z1; % Enthalpy balance resid = CH4_in*H_CH4_in + O2_in*H_O2_in + H2_in*H_H2_in +... H2O_in*H_H2O_in + CO_in*H_CO_in + CO2_in*H_CO2_in -... (CH4*H_CH4 + O2*H_O2 + H2*H_H2 + H2O*H_H2O + CO*H_CO + CO2*H_CO2); return %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function [CH4_coef, O2_coef, H2_coef, H2O_coef, CO_coef, CO2_coef] = shomate_coef_wgs() % this just stores the shomate coefficients so I don't have to looks at % them, and returns the coef matrices that can be used in Shomate.m CH4_coef = [... -0.703029 85.81217 108.4773 11.26467 -42.52157 -2.114146 5.862788 0.138190 0.678565 -26.42221 -76.84376 -153.5327 158.7163 224.4143 300 1300 1300 6000]; O2_coef = [... 29.65900 6.137261 -1.186521 0.095780 -0.219663 -9.861391 237.9480 298 6000]; H2_coef = [... 33.066178 18.563083 43.413560 -11.363417 12.257357 -4.293079 11.432816 -2.859786 1.272428 -2.772874 0.268238 -0.096876 -0.158558 1.977990 -20.533862 -9.980797 -1.147438 -38.515158 172.707974 156.288133 162.081354 298 1000 2500 1000 2500 6000]; H2O_coef = [... 30.09200 41.96426 6.832514 8.622053 6.793435 -1.499780 -2.534480 0.098119 0.082139 -11.15764 -250.8810 -272.1797 223.3967 219.7809 300 1700 1700 6000]; CO_coef = [... 25.56759 35.15070 6.096130 1.300095 4.054656 -0.205921 -2.671301 0.013550 0.131021 -3.282780 -118.0089 -127.8375 227.3665 231.7120 298 1300 1300 6000]; CO2_coef = [... 24.99735 58.16639 55.18696 2.720074 -33.69137 -0.492289 7.948387 0.038844 -0.136638 -6.447293 -403.6075 -425.9186 228.2431 263.6125 298 1200 1200 6000]; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function [Cp,H,S]=Shomate(T,molecule_data) % Program to compute thermo stored as Shomate coefficients % W.H. Green, 9/27/02 % Outputs: % 1) constant-pressure heat capacity Cp in J/mol K % 2) enthalpy in kJ/mol % 3) entropy in J/mol K % Note, if T is a vector, these will all be vectors of the same shape. % Inputs: T in K, molecule_data = [A B C D E F G Tmin Tmax]' % A-G, Tmin, Tmax from NIST Webbook http://webbook.nist.gov % in Joule/mole Kelvin units % Cp0= A + B*t + C*t^2 + D*t^3 + E/t^2 % H0= A*t + B*t^2/2 + C*t^3/3 + D*t^4/4 - E/t + F % S0= A*ln(t) + B*t + C*t^2/2 + D*t^3/3 - E/(2*t^2) + G % Tmin and Tmax define the range over which % these A-G are accurate % Sometimes there are different parameter sets for different T ranges % if so, make molecule_data a matrix, with a column for each T range. % First: figure out how many parameter sets we have % and how many T's the user has input dim_molecule_data=size(molecule_data); number_of_parameter_sets=dim_molecule_data(2); dim_T=size(T); number_of_T=dim_T(1); Cp=-1e20*ones(dim_T); % if you see any -1e20's in output, something went wrong H=Cp; % most likely your T was outside the T range where the S=Cp; % Shomate fit is valid. for i=1:number_of_parameter_sets Tmin=molecule_data(8,i); Tmax=molecule_data(9,i); for j=1:number_of_T if((T(j)>=Tmin)&(T(j)<=Tmax)) A = molecule_data(1,i); B=molecule_data(2,i); C = molecule_data(3,i); D=molecule_data(4,i); E = molecule_data(5,i); F=molecule_data(6,i); G = molecule_data(7,i); t=0.001*T(j); % NIST t= T/1000 Kelvin Cp(j,1)=A+(B*t)+(C*t*t)+(D*(t^3))+(E/(t*t)); H(j,1)=(A*t)+(0.5*B*t*t)+(C*(t^3)/3)+(0.25*D*(t^4))-(E/t)+F; S(j,1)=(A*log(t))+(B*t)+(0.5*C*t*t)+(D*(t^3)/3)-(0.5*E/(t*t))+G; end end end % if any of the values are -1e20, it means that the Shomate polynomial % was not valid in the requested T range. % % Sample input % For H2 % H2_data=[33.10780 -11.50800 11.60930 -2.844400 -0.159665 -9.991971 172.7880 298 1500]'; % By running [Cp,H,S]=Shomate(300, H2_data); % % Sample output % » Cp % Cp = % 28.8494 % » H % H = % 0.0534 % » S % S = % 130.8586