% 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);