% Population Balance Homework Problem || % Brendan Hoffman, Kelly Tipton, Yechun Wang, Matt McHale || % modifications by Sheryl Ehrman, Fall 2005 % ========================================================== clear close all global d0 mfr T Po MW mu smf rhos kb R rhog Vo SECT rhomix nu startnum D L k m global uz A funcall eps relerr=1e-5; funcall=0; % Set parameters up in a structure variable 'params' d0=0.25e-6; % microns, primary particle diameter mfr=3.5e4/2.02/3600; % total mass flow rate, kg/s T=200+273; % gas temperature, K Po=(10+14.7)*101325/14.7; % pressure, Pa MW=44.5e-3; % molecular weight, kg/mol mu=1.57e-5; % viscocity, Pa s smf=0.15; % mass fraction solids rhos=2200.0; % solids density, kg/m^3 kb=1.38e-23; % boltzmann's constant (J/s) R=8.3145; % gas constant, J/mol K rhog=Po*MW/R/T; % gas density, kg/m^3 Vo=pi/6*d0^3; % initial particle volume SECT=28; % number of sectionalization segments dpc=24e-6; % cut size, m rhomix=smf*rhos+(1-smf)*rhog; % mixture density nu=mu/rhomix; % kinematic viscocity startnum=mfr*smf/rhos/Vo; % starting number of particles %Pipe diameter ranges D=[0.254 0.381 0.508]; %Pipe Length ranges (m) L=609.6; %Velocities A=pi/4.*D.^2; uz=mfr./rhomix./A; Re=rhomix.*uz.*D/mu; eps=2*0.04./(Re.^0.16).*uz.^3./D; % Energy dissipation factor k=2; m=1; No=zeros(SECT,1); No(1)=startnum; [Z,N]=ode15s('numdist',[0 L(m)],No); % Calculate the efficiency V=zeros(SECT,1); d=zeros(SECT,1); eff=zeros(SECT,1); efftot=zeros(length(N),1); for j=1:length(N) sums=zeros(2,1); for i=1:SECT V(i)=2^(i-1)*(3*Vo/2); d(i)=d0*(3*2^(i-2))^(1/1.8); eff(i)=(d(i)/dpc)^2/(1+(d(i)/dpc)^2); sums(1)=sums(1)+eff(i)*V(i)*N(j,i); sums(2)=sums(2)+V(i)*N(j,i); end efftot(j)=sums(1)/sums(2)*100; end %Find where the efficiency is 75% (by index) for j=1:length(N)-1 if (efftot(j)-75 <=0 & efftot(j+1) - 75 >=0) J=j; break; end end % Check mass closure on last legnth segment m0=0; for i=1:SECT m0=m0+V(i)*N(length(N),i); end m0=m0/(3*Vo/2); n=N(length(N),:); check=0; for i=1:SECT check=check+2^(i-1)*n(i)/m0; end if (abs(check-1) < relerr) disp('Mass closure achieved'); end % Output the length for a chosen D Lcut=Z(J); disp(['For D = ', num2str(D(k)), ' 75% efficiency occurs at L = ', num2str(Lcut),'.']) % Plot the number distribution (by bin) at the 75% cut-off sizes=zeros(SECT,1); for i=1:SECT sizes(i)=i; end bar(sizes,N(J,:)); title('Size distribution at 75% cutoff'); xlabel('Size bin number'); ylabel('Particle count');