% filename rxnFD.m % This code solves the equations for reaction-diffusion % in a porous catalyst using the finite difference method. % a = 1, 2, 3 for planar, cylindrical, and spherical geometry % If the problem is linear, set linear = 1 to avoid a second iteration. geometry = 3 linear = 0 ngrid = 6 %set ngrid before calling % ngrid must be a number*2+1 (for Simpson's rule) % set phi before calling %phi = 10 phi2 = phi*phi; dr = 1/(ngrid-1); dr2=dr*dr; % iterate tol = 1.e-12; change = 1.; iter=0; % initial guess for i=1:ngrid c(i) = 0.0; rad(i) = (i-1)/(ngrid-1); end while change>tol iter=iter+1; % Set up the matrices % set the matrices to zero for i=1:ngrid f(i) = 0.; for j=1:ngrid aa(i,j) = 0.; end end % set the matrices for the differential equation for j=2:ngrid-1 ans = rate(c(j)); raterxn = ans(1); drate = ans(2); aa(j,j-1) = 1 - 2.*dr/rad(j); aa(j,j) = -2 - phi2*dr2*drate; aa(j,j+1) = 1 + 2.*dr/rad(j); f(j) = phi2*dr2*(raterxn - drate*c(j)); end f(ncol+1) = 1; ans = rate(c(1)); aa(1,1) = -1 - dr2*phi2*ans(2)/6.; aa(1,2) = 1; f(1) = phi2*dr2*(ans(1) - ans(2)*c(1))/6.; aa(ngrid,ngrid) = 1.; f(ngrid)=1.; % Solve one iteration. c11 = aa\f'; % Calculate the criterion to stop the iterations. sum=0.; for i=1:ngrid sum = sum + abs(c11(i) - c(i)); c(i) = c11(i); end if linear==1,break,end change = sum end % Plot the solution. plot(rad,c,'*-b') xlabel('r') ylabel('c') c % calculate the effectiveness factor for i=2:2:ngrid-1 w(i) = 4.; w(i+1) = 2.; end w(1) = 1; w(ngrid) = 1; sum = 0.; for i=1:ngrid ans = rate(c(i)); raterxn = ans(1); num1 = rad(i)^(geometry-1); sum = sum + w(i)*raterxn*num1; end eta = geometry*(dr/3)*sum % rate.m % calculate the rate and derivative for a given concentration function ans = rate(c) % use for linear reaction %ans(1) = c; %ans(2) = 1.; % use for second order reaction ans(1) = c*c; ans(2) = 2.*c; % use for non-isothermal reaction %ggg = 18.; bbb = 0.3; %ggg = 30; bbb = 0.4; %t = 1 + bbb - bbb*c; %eee = exp(ggg*(1 - 1/t)); %ans(1) = c*eee; %ans(2) = eee*(1 - bbb*ggg*c/(t*t));