% 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');