Produced using MATLAB® software. % CSTR_ODE\CSTR_ODE Page 1 of 52 % CSTR_ODE\CSTR_ODE.m % % function iflag_main = CSTR_ODE(); % % This MATLAB program calculates the dynanamic concentration and % temperature profiles within a continuous stirred-tank reactor. % This program models a CSTR in overflow mode (no volume change), % or, if the flow rate is set to zero, the model simulates a batch % reactor. At the moment, semi-batch reactors with a volume % that changes with time is not supported. % % This program uses the MATLAB ODE-IVP solver ode15s, a multistep % method that works well for stiff systems. % % K. Beers % MIT ChE % 10/28/2001 function iflag_main = CSTR_ODE(); iflag_main = 0; % First, read in the problem definition data from the keyboard or % from an input file. imode_input = input('Select mode of input (0 = keyboard, 1 = file) : '); if(~imode_input) [ProbDim,Reactor,Inlet,Physical,Rxn,StateInit,iflag] = ... read_program_input; if(iflag <= 0) iflag_main = -1; error(['CSTR_ODE: Error returned from read_program_input = ', ... int2str(iflag)]); end else file_name = input('Enter input file name (without .m) : ','s'); [ProbDim,Reactor,Inlet,Physical,Rxn,StateInit,iflag] = ... feval(file_name); end % We now stack the initial state data into the state vector. 7/16/2002% CSTR_ODE\CSTR_ODE Page 2 of 52 x_init = stack_state(StateInit,ProbDim.num_species, ... Reactor.nonisothermal); % Then, the options for the MATLAB ODE-IVP solver are set. options = odeset('RelTol',1e-4,'AbsTol',1e-6, ... 'Refine',10,'Stats','on'); % The user is now asked to the length of the time period to % be simulated. % time_final prompt = 'Enter final simulation time : '; check_real=1; check_sign=1; check_int=0; time_final = get_input_scalar(prompt, ... check_real,check_sign,check_int); % Now, the ODE-IVP solver is called. [t,x_traj] = ode15s(@CSTR_ODE_calc_f, ... [0 time_final],x_init, ... options,ProbDim,Reactor,Inlet,Physical,Rxn); % Plot the simulation results. for ispecies=1:ProbDim.num_species figure; plot(t,x_traj(:,ispecies)); title(['Concentration of species # ', ... int2str(ispecies)]); xlabel('time'); ylabel(['c (',int2str(ispecies),')']); end if(Reactor.nonisothermal) figure; plot(t,x_traj(:,ProbDim.num_species+1)); title('Reactor Temperature'); xlabel('time'); ylabel('T'); end % Save the results to a binary output file. save CSTR_ODE_results.mat; iflag_main = 1; return; 7/16/2002% CSTR_ODE\CSTR_ODE Page 3 of 52 % CSTR_ODE\read_program_input.m % % function [ProbDim,Reactor,Inlet,Physical,Rxn,StateInit,iflag] = ... % read_program_input(); % % This procedure reads in the simulation parameters that are % required to define a single CSTR dynamic simulation. % % Kenneth Beers % Massachusetts Institute of Technology % Department of Chemical Engineering % 10/28/2001 % % Version as of 10/28/2001 function [ProbDim,Reactor,Inlet,Physical,Rxn,StateInit,iflag] = ... read_program_input(); func_name = 'read_program_input'; iflag = 0; disp(' '); disp(' '); disp('Input the system parameters : '); disp('The parameters are input in real units, where '); disp(' L = unit of length'); disp(' M = unit of mass'); disp(' t = unit of time'); disp(' E = unit of energy'); disp(' T = unit of temperature'); % REACTOR DATA -----------------------------------------% PDL> Input first the reactor volume and heat transfer data : disp(' '); disp(' '); disp('Input the reactor data : '); disp(' '); % Perform assertion that a real scalar positive number has % been entered. This is performed by a function assert_scalar 7/16/2002% CSTR_ODE\CSTR_ODE Page 4 of 52 % that gives first the value and name of the variable, the name % of the function that is making the assertion, and values of % 1 for the three flags that tell the assertion routine to make % sure the value is real, positive, and not to check that it is % an integer. disp(' '); disp('Reactor size information : '); % Reactor.volume check_real=1; check_sign=1; check_int=0; prompt = 'Input the volume of the reactor (L^3) : '; Reactor.volume = get_input_scalar(prompt, ... check_real,check_sign,check_int); % Reactor.F_vol check_real=1; check_sign=1; check_int=0; prompt = 'Input volumetric flowrate through reactor (L^3/t) : '; Reactor.F_vol = get_input_scalar(prompt, ... check_real,check_sign,check_int); % Now, the user selects whether to perform a non-isothermal % calculation or to keep the temperature constant at the % value of Temp_cool. % Reactor.nonisothermal check_real=1; check_sign=1; check_int=2; prompt = 'Simulate isothermal (0) or non-isothermal (1) reactor? : '; Reactor.nonisothermal = get_input_scalar(prompt, ... check_real,check_sign,check_int); % If operate in isothermal mode, only input set temperature. if(~Reactor.nonisothermal) % Reactor.Temp_cool check_real=1; check_sign=1; check_int=0; prompt = 'Enter the constant reactor temperature : '; Reactor.Temp_cool = get_input_scalar(prompt, ... check_real,check_sign,check_int); % We now set the other non-needed coolant information % to dummy values. Reactor.F_cool = 1; Reactor.U_HT = 1; Reactor.A_HT = 1; Reactor.Cp_cool = 1; % otherwise, if simulating non-isothermal reactor, input 7/16/2002% CSTR_ODE\CSTR_ODE Page 5 of 52 % the data for the coolant jacket. else disp(' '); disp('Reactor coolant information : '); % Reactor.F_cool check_real=1; check_sign=1; check_int=0; prompt = 'Input the coolant flowrate (L^3/t) : '; Reactor.F_cool = get_input_scalar(prompt, ... check_real,check_sign,check_int); % Reactor.Temp_cool check_real=1; check_sign=1; check_int=0; prompt = 'Input the coolant inlet temperature (T) : '; Reactor.T_cool = get_input_scalar(prompt, ... check_real,check_sign,check_int); % Reactor.U_HT check_real=1; check_sign=1; check_int=0; prompt = 'Input the jacket heat transfer coefficient (E/t/(L^2)/T) : '; Reactor.U_HT = get_input_scalar(prompt, ... check_real,check_sign,check_int); % Reactor.A_HT check_real=1; check_sign=1; check_int=0; prompt = 'Input the jacket heat transfer area (L^2) : '; Reactor.A_HT = get_input_scalar(prompt, ... check_real,check_sign,check_int); % Reactor.Cp_cool check_real=1; check_sign=1; check_int=0; prompt = 'Input the reactor coolant heat capacity (E/mol/T) : '; Reactor.Cp_cool = get_input_scalar(prompt, ... check_real,check_sign,check_int); end % PDL> Input number of species, ProbDim.num_species disp(' '); disp(' '); % ProbDim.num_species check_real=1; check_sign=1; check_int=1; prompt = 'Input the number of species : '; ProbDim.num_species = get_input_scalar(prompt, ...
View Full Document
Unlocking...