% 10.34 - Fall 2006
% HW Set#2, Problem #4
% Rob Ashcraft - Sept. 11, 2006

% This script uses two other functions: WaterGasShift and WGS_contour
% This takes a vector of inlet reactor temperatures and attempts to solve
% for the steady-state reactor temperature and outlet mole fractions for
% the simple reaction: CO + H2O <=> CO2 + H2
% It then extracts the results at Tin = 700K for the given conditions:
%   x_in = [CO H2O CO2 H2] = [0.2, 0.8, 0, 0] (mole fractions)
%   P = 10 atm, Vreactor = 5 L 

clear; clc;
start_time = cputime;

% define global variables for the constants in the problem
global Vreactor res_time R Rj A Ea dHrxn dSrxn Cp_vec h

P = 10;  % atm
Vreactor = 5;   % Liters
res_time = 1.5;  % seconds
h = 1.2;  % heat tranfer coefficient in Watts/K

R = 0.08206;  % atm-L/mole-K
Rj = 8.314;  % J/mole-K
A = 32;
Ea = 53000;
dHrxn = -41000;
dSrxn = -42;
Cp_vec = [3.5*Rj, 4*Rj, 4*Rj, 3.5*Rj];

% This is the vector of inlet that we want to investigate in order to find
% the best operating inlet Temperature and the resulting conversion
Tin = [500:25:1500]';

% inlet mole fraction values for this simulation
x_in = [0.2;0.8;0;0];

% loop over the temperature vector in order to get the results at each
% temp, also keep track of the number of steps and time on screen
for i=1:length(Tin)
    timein = cputime;
    [x_out(:,i), Treactor(i,1)] = WaterGasShift(P,Tin(i),x_in);
    conversion(i) = (x_in(1) - x_out(1,i))/x_in(1);
    
    step_time = cputime - timein;
    total_time = cputime - start_time;
    step_num = num2str(i);
    disp(['Completed Step Number ',step_num,' of ',num2str(length(Tin))]);
    disp(['Step Time (s): ',num2str(step_time)]);
    disp(['Total Time (s): ',num2str(total_time)]);
    disp('  ');
    
end

% extract the data for Tin = 700 K to report as asked for in the HW
loc_700 = find(Tin == 700);
x_out_700 = x_out(:,loc_700);
Treactor_700 = Treactor(loc_700);

disp('  ');
disp('Results for Inlet Temperature = 700 K');
disp(['Reactor Temperature = ',num2str(Treactor_700),' K']);
disp(['Outlet Mole Fraction of CO = ',num2str(x_out_700(1))]);
disp(['Outlet Mole Fraction of H2O = ',num2str(x_out_700(2))]);
disp(['Outlet Mole Fraction of CO2 = ',num2str(x_out_700(3))]);
disp(['Outlet Mole Fraction of H2 = ',num2str(x_out_700(4))]);
disp('  ');


%determine the maximum conversion of CO and the operating conditions
[max_conv,loc] = max(conversion);
Tin_conv_max = Tin(loc);
Treactor_conv_max = Treactor(loc);

disp('Operating Condition to Achieve Maximum CO conversion for the current conditions');
disp(['Inlet Temperature = ',num2str(Tin_conv_max),' K']);
disp(['Reactor Temperature = ',num2str(Treactor_conv_max),' K']);
disp(['Maximum Conversion of CO = ',num2str(max_conv*100),' %']);
disp(['Outlet Mole Fraction of H2O = ',num2str(x_out_700(2))]);
disp('  ');

% Plot various figures of interest
figure;
plot(Tin,conversion);
xlabel('Inlet Temperature (K)'); ylabel('CO Conversion Fraction');

figure;
plot(Treactor,conversion);
xlabel('Reactor Temperature (K)'); ylabel('CO Conversion Fraction');

figure;
plot(Tin,Treactor-Tin,'-',[500 1500],[0 0],':');
xlabel('Inlet Temperature  (K)'); ylabel('T_R_e_a_c_t_o_r - T_I_n_l_e_t  (K)');

% This part was not required by the problem, but is an interesting way to
% look at a 2-D problem.
%
% Look at the problem graphicly, by plotting contours of the residual
% equations.  The solution(s) will be where the zero contours intersect

T_in = 700;  % arbitrary inlet temperatue for this case

% specify a grid size.  creates an N x N grid of Treactor and z points at
% which to evaluate the residual equations.  The results will be a matrix
% that can be plotted as a surface of contour plot.  
gridsize = 200; 
n_CO_in = x_in(1)*P*Vreactor/0.08206/T_in;
z_vec = linspace(0,n_CO_in*1.5,gridsize);
Tr_vec = linspace(500,3*T_in,gridsize);

WGS_contour(z_vec,Tr_vec,P,T_in,x_in);