% tubular_reactor_2rxn_SS.m % This MATLAB program computes the concentration profiles within a % tubular rector for the system of two elementary reactions % A + B --> C % C + B --> D % The model includes both convection and dispersion, and uses % upwind finite differences to discretize the system at a grid % of N uniformly-spaced points within the reactor. % K. Beers. MIT ChE. 10/6/03 function iflag_main = tubular_reactor_2rxn_SS(); iflag_main = 0; % set problem specifications Param.L = 10; % reactor length Param.k = [1; 1]; % rate constants for rxn. 1 and 2 Param.c_in = [1; 1; 0; 0]; % inlet conc. of A,B,C,D Param.v = 1; % axial velocity in reactor Param.D = 1e-4; % dispersion coefficient Param.N_f = 4; % total number of fields % Next, give name of user-supplied routine that computes the reaction rates % and the derivatives of the reaction rates with respect to each % concentration. Param.scheme_rxn = 'scheme_2rxn'; % set problem dimensions Grid.N = 250; % 100 grid points Grid.dz = Param.L/(Grid.N + 1); % uniform grid spacing Grid.z = linspace(Grid.dz,Param.L-Grid.dz,Grid.N)'; % set uniform grid Grid.N_var = Grid.N*Param.N_f; % total # of unknowns Grid.Pe_loc = Param.v*Grid.dz/Param.D; % local Peclet number of grid % as initial guess, use uniform profiles at inlet values x0 = zeros(Grid.N_var,1); % allocate memory for ifield = 1:Param.N_f % for each field for k=1:Grid.N % for each grid point n = tubular_reactor_get_label(ifield,k,Param.N_f); x0(n) = Param.c_in(ifield); end end % call fsolve to find the steady state concentration profiles Options_NAE = optimset('LargeScale','off','Jacobian','on'); [x,f_val,exitflag] = fsolve(@tubular_reactor_2rxn_SS_calc_f,x0,Options_NAE,Param,Grid); if(exitflag <= 0) error(['fsolve returned exitflag = ', int2str(exitflag)]); end % plot results figure; % extract and plot concentration profile of A (use BC to generate % values at the boundaries) z_plot = [0; Grid.z; Param.L]; % extended to include inlet, outlet positions field_val = extract_field(1,x,Grid,Param); plot(z_plot,field_val); hold on; % extract and plot concentration profile of B field_val = extract_field(2,x,Grid,Param); plot(z_plot,field_val,'--'); % extract and plot concentration profile of C field_val = extract_field(3,x,Grid,Param); plot(z_plot,field_val,'-.'); % extract and plot concentration profile of D field_val = extract_field(4,x,Grid,Param); plot(z_plot,field_val,':'); legend('[A]','[B]','[C]','[D]'); xlabel('z'); ylabel('species concentrations'); title('Tubular reactor with A + B --> C, B + C --> D'); iflag_main = 1; return; % --------------------------------------------------------------------- % tubular_reactor_get_label() % This routine returns the index that identifies the location in the % state vector of the concentation at a specific interior grid point function label = tubular_reactor_get_label(ifield,k,N_f); label = N_f * (k - 1) + ifield; return; % --------------------------------------------------------------------- % extract_field() % This routine extracts from the state vector the concentration values % for a particular species from the master state vector. It also uses % the Danckwert's BC to compute the inlet and outlet concentrations. function field_val = extract_field(ifield,x,Grid,Param); field_val = zeros(Grid.N + 2, 1); % allocate space in memory for k=1:Grid.N n = tubular_reactor_get_label(ifield,k,Param.N_f); field_val(1+k) = x(n); end % add outlet value from BC field_val(Grid.N + 2) = field_val(Grid.N + 1); % add inlet value from BC field_val(1) = ( Grid.Pe_loc * Param.c_in(ifield) + field_val(2) ) / (1 + Grid.Pe_loc); return; %---------------------------------------------------------------------- % This routine computes the function vector and the Jacobian of the % tubular reactor model. It calls a user-supplied routine that % computes the rates of each reaction, the derivatives of the % reaction rates with respect to each component, and returns the % stoichiometric coefficients. From this information, the routine % computes the function vector and the Jacobian. function [f,Jac] = tubular_reactor_2rxn_SS_calc_f(x,Param,Grid); % allocate space in memory f = zeros(Grid.N_var,1); nnz_Jac_max = Grid.N_var * (Param.N_f + 2); Jac = spalloc(Grid.N_var,Grid.N_var,nnz_Jac_max); % set common coefficients for lower and upper grid points for discretizing % the transport terms A_lo = Param.v/Grid.dz + Param.D/(Grid.dz^2); A_mid = -Param.v/Grid.dz - 2*Param.D/(Grid.dz^2); A_hi = Param.D/(Grid.dz^2); % for each point for k=1:Grid.N % extract concentrations of each field variable at grid point % conc_loc = ( [A], [B], [C], [D] ) conc_loc = zeros(Param.N_f,1); for ifield=1:Param.N_f n = tubular_reactor_get_label(ifield,k,Param.N_f); conc_loc(ifield) = x(n); end % compute the reaction rates and the partial derivatives of each % rate with respect to each local concentration. This is done in % the routine supplied by the user whose name is stored in the % character string Param.rxn_scheme. [num_rxn,rxn_rate,stoich_coeff,rxn_rate_deriv] = ... feval(Param.scheme_rxn,conc_loc,Param); % compute local source term values of each species source = zeros(Param.N_f,1); for ifield=1:Param.N_f for irxn=1:num_rxn source(ifield) = source(ifield) + rxn_rate(irxn)*stoich_coeff(ifield,irxn); end end % compute source term Jacobianm matrix for which element (i,j) % is the partial derivative of the source term for field # i % with respect to the local value of field # j source_Jac = zeros(Param.N_f,Param.N_f); for ifield=1:Param.N_f for irxn=1:num_rxn for jfield=1:Param.N_f source_Jac(ifield,jfield) = source_Jac(ifield,jfield) + ... rxn_rate_deriv(jfield,irxn)*stoich_coeff(ifield,irxn); end end end % Now, initialize function vector and Jacobian elements for this grid % point to contain the source term and the source term Jacobian for ifield=1:Param.N_f n = tubular_reactor_get_label(ifield,k,Param.N_f); f(n) = source(ifield); for jfield=1:Param.N_f % for all fields m = tubular_reactor_get_label(jfield,k,Param.N_f); Jac(n,m) = source_Jac(ifield,jfield); end end % Now, get transport contribution to function and Jacobian values for ifield=1:Param.N_f n = tubular_reactor_get_label(ifield,k,Param.N_f); % dependence of transport on field value itself f(n) = f(n) + A_mid*x(n); Jac(n,n) = Jac(n,n) + A_mid; % dependence on value to the left if(k > 1) % if not the left-most point m = tubular_reactor_get_label(ifield,k-1,Param.N_f); % value of same field at left point f(n) = f(n) + A_lo*x(m); Jac(n,m) = Jac(n,m) + A_lo; else % use Danckwerts' BC to estimate conc. at inlet c_lo = (Grid.Pe_loc*Param.c_in(ifield) + x(n) ) / (1 + Grid.Pe_loc); f(n) = f(n) + A_lo*c_lo; Jac(n,n) = Jac(n,n) + A_lo/(1+Grid.Pe_loc); end % dependence on value to the right if(k < Grid.N) % if not the right-most point m = tubular_reactor_get_label(ifield,k+1,Param.N_f); % value of same field at right point f(n) = f(n) + A_hi*x(m); Jac(n,m) + Jac(n,m) + A_hi; else % use outlet BC to estimate conc. at outlet c_hi = conc_loc(ifield); % same as current value f(n) = f(n) + A_hi*c_hi; Jac(n,n) = Jac(n,n) + A_hi; end end end return; % ---------------------------------------------------------------------- % This routine computes the rates of each reaction, sets the stoichiometric % coefficients, and computes the derivative of each reaction rate with % respect to the concentrations of each species. function [num_rxn,rxn_rate,stoich_coeff,rxn_rate_deriv] = ... scheme_2rxn(conc_loc,Param); % compute rates of each reaction and set stiochoimetric coefficients num_rxn = length(Param.k); rxn_rate = zeros(1,num_rxn); stoich_coeff = zeros(Param.N_f,num_rxn); % reaction 1, A + B --> C, (1) + (2) --> (3) rxn_rate(1) = Param.k(1)*conc_loc(1)*conc_loc(2); % local rate of rxn. 1 stoich_coeff(:,1) = [-1; -1; 1; 0]; % reaction 2, B + C --> D, (2) + (3) --> (4) rxn_rate(2) = Param.k(2)*conc_loc(2)*conc_loc(3); % local rate of rxn. 2 stoich_coeff(:,2) = [0; -1; -1; 1]; % compute partials of each rxn. rate with respect to concentrations. rxn_rate_deriv = zeros(Param.N_f,num_rxn); % reaction 1, A + B --> C, (1) + (2) --> (3) rxn_rate_deriv(1,1) = Param.k(1)*conc_loc(2); % deriv. of rate # 1 w.r.t. [A] rxn_rate_deriv(2,1) = Param.k(1)*conc_loc(1); % deriv. of rate # 1 w.r.t. [B] % reaction 2, B + C --> D, (2) + (3) --> (4) rxn_rate_deriv(2,2) = Param.k(2)*conc_loc(3); % deriv. or rate # 2 w.r.t. [B] rxn_rate_deriv(3,2) = Param.k(2)*conc_loc(2); % deriv. of rate # 2 w.r.t. [C] return;