% 10.34 - Fall 2006 % HW Set#2, Problem #4a % Rob Ashcraft - Sept. 11, 2006 % This function solves for the reactor temperature and outlet stream % composition for a single reaction CSTR for the water-gas shift reaction. % The equations are written in terms of a single progress variable and the % reactor temperature. function [x_out, Treactor] = WaterGasShift(P,Tin,x_in) % inputs: % P = reactor pressure % Tin = inlet reactor temperature (K) % x_in = [CO H2O CO2 H2] (mole fractions) % outputs: % x_out = outlet mole fractions % Treactor = reactor temperature % define global variables for the constants in the problem global Vreactor res_time R Rj A Ea dHrxn dSrxn Cp_vec h x_CO_in = x_in(1); x_H2O_in = x_in(2); x_CO2_in = x_in(3); x_H2_in = x_in(4); PCOin = P*x_CO_in; PH2Oin = P*x_H2O_in; PCO2in = P*x_CO2_in; PH2in = P*x_H2_in; % Choose the initial guesses for Treactor and the progress variable (z). z % is the number of moles of consumed in the reactor (extent of reaction). % The initial reactor temperature is guessed to be 105% of the inlet. The % progress variable is guessed as 10% of the inlet CO composition. Treactor_0 = Tin*1.05; z_0 = P*Vreactor/R/Tin*x_in(1)*0.1; var_0 = [z_0; Treactor_0]; % Options for the solver... probably not needed for this case. options = optimset('MaxFunEvals',20000,'MaxIter',10000,'Display','final'); % call fsolve to solve the problem var = fsolve(@WGS_solver,var_0,options,P,Tin,x_in); z = var(1); Treactor = var(2); % determine the outlet mole fractions from the progress variable xCO = (PCOin - z*R*Treactor/Vreactor)/P; xH2O = (PH2Oin - z*R*Treactor/Vreactor)/P; xCO2 = (PCO2in + z*R*Treactor/Vreactor)/P; xH2 = (PH2in + z*R*Treactor/Vreactor)/P; x_out = [xCO; xH2O; xCO2; xH2]; %%%%%%%%%%%%%%%%%%%%%%%%% % This is the function that evaluates the residuals of the equations for % the system function resid = WGS_solver(var,P,Tin,x_in) % inputs: % var = [z, Treactor] = system variables to be solved % P, Tin, x_in = described above % outputs: % resid = [resid1, resid2] residuals for the two equations % define global variables for the constants in the problem global Vreactor res_time R Rj A Ea dHrxn dSrxn Cp_vec h z = var(1); Treactor = var(2); Cp_in = Cp_vec*x_in; dGrxn = dHrxn - Treactor*dSrxn; x_CO_in = x_in(1); x_H2O_in = x_in(2); x_CO2_in = x_in(3); x_H2_in = x_in(4); PCOin = P*x_CO_in; PH2Oin = P*x_H2O_in; PCO2in = P*x_CO2_in; PH2in = P*x_H2_in; % convert partial Pressures to progress variable PCO = PCOin - z*R*Treactor/Vreactor; PH2O = PH2Oin - z*R*Treactor/Vreactor; PCO2 = PCO2in + z*R*Treactor/Vreactor; PH2 = PH2in + z*R*Treactor/Vreactor; % Write out the rate equation for the system rate = A*exp(-Ea/Rj/Treactor)*(PCO*PH2O - PH2*PCO2*exp(dGrxn/Rj/Treactor))/(1 + (PCO/0.2)); % calculate the residuals for the conversion equation and energy balance resid_1(1,1) = z - rate*res_time; resid_1(2,1) = -dHrxn*rate - (P*Vreactor/R/Treactor)/res_time*(Cp_in*(Treactor-Tin)) - h*(Treactor-300); resid = (resid_1);