Solving Population Balance Equation for multiple initial particle dimensions

10 views (last 30 days)
I am trying to develop a code for a brownian kernel agglomeration of multiple particle dimensions at once to produce a particle size distribution. I am getting the error 'Unable to perform assignment because the left and right sides have a different number of elements' and I am struggling to diagnose how to bypass this.
Has anyone got any ideas on how to get this working or can suggest any ideas for my code?
Thanks
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= [3e-6 5e-6 10e-6 25e-6 50e-6 100e-6] % microns, primary particle diameter
mfr=0.847847; % total mass flow rate, m/s
T=20+273; % gas temperature, K
Po=101325 % pressure, Pa
mu=1.81e-5; % viscocity, Pa s
kb=1.38e-23; % boltzmann's constant (J/s)
R=8.3145; % gas constant, J/mol K
SECT=28; % number of sectionalization segments
Vo=pi/6*d0.^3; % initial volume
nu=1.48e-5; % kinematic viscocity
startnum=[1000 100 50 20 6 2] % starting number of particles
D=[0.110 0.110 0.110]; % Pipe Diameter ranges (m)
L=[1.2]; % Pipe Length ranges (m)
A=pi./4.*D.^2; % Pipe Surface Area
uz=0.69212; % Air Flow Speed (m3/s)
Re= 2300;
eps=2*0.04./(Re.^0.16).*uz.^3./D; % Energy dissipation factor
k=1;
m=1;
%No=zeros(6,28)%
%No(1:6)=startnum
No=zeros(28,6);
No(1:28:end)=startnum;
[Z,N]=ode15s('numdist',[0:L],No);
Calculate the efficiency
V=zeros(SECT,6);
d=zeros(SECT,6);
eff=zeros(SECT,6);
efftot=zeros(length(N),6);
for j=1:length(N)
sums=zeros(2,6);
for i=1:SECT
V(i)=2^(i-1)*(3*Vo/2);
d(i)=d0.*(3*2^(i-2))^(1/1.8);
sums(1)=sums(1)*V(i)*N(j,i);
sums(2)=sums(2)+V(i)*N(j,i);
end
% efftot(j)=sums(1)/sums(2)*100;
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
% Plot the number distribution (by bin) at the 75% cut-off
sizes=zeros(SECT,6);
for i=1:SECT
sizes(i)=i;
end
bar(sizes,N(24,:));
title('Size distribution at 75% cutoff');
xlabel('Size bin number');
ylabel('Particle count');
function [dn]=numdist(z,n)
% k and l refer to the array indexes in the D and L
% variables, specifying which of the DL combinations
% we are solving for on this function call
global d0 mfr T Po MW mu smf rhos kb R rhog Vo SECT rhomix nu startnum D L k l m No
global uz A funcall
funcall=funcall+1;
dn=zeros(SECT,6);
sums=zeros(6,3);
disp(['Function call ', num2str(funcall),' z value ', num2str(z)...
,' of ',num2str(L(m)),'.'])
for i=1:168
if (i-2 >= 1)
for j=1:i-2
sums(1:3:18)= sums(1:3:18)+beta(i-1,j)/2^(i-j-1)*n(j);
end
end
if (i-1 >= 1)
for j=1:i-1
sums(2:3:18)=sums(2:3:18)+beta(i,j)/2^(i-j)*n(j);
end
dn(i)=dn(i)+ n(i-1)*sum(1)+ 0.5*beta(i-1,i-1)*n(i-1)^2-n(i)*sums(2:3:18)*n(i);
end
if (i+1 <=168)
dn(i)=dn(i)*n(i+1);
for j=1:168
sums(3:3:18)=sums(3:3:18)+beta(i,j)*n(j);
end
dn(i)=dn(i)-n(i)*sums(3:3:18);
end
dn(i)=dn(i)/uz(k);
sums=zeros(6,3);
end
function B=beta(i,j)
% beta.m --> Returns the value of the sectionalized coagulation
% kernerl, beta, for the particle population balance
% design problem
global d0 mfr T Po MW mu smf rhos kb R rhog Vo SECT rhomix nu startnum D L k m
global uz A eps
B=2*kb*T/3/mu*(2+2^((i-j)/3)+2^((j-i)/3))+ 0.31*(eps(m)/nu)^0.5.*...
(3.*Vo./2)*(2^(i-1)+3*(2^((2*i+j)/3-1)+2^((i+2*j)/3-1))+2^(j-1));
  4 Comments

Sign in to comment.

Answers (1)

darova
darova on 6 Apr 2021
First of all - timespan 0:1.2 is just [0 1]
L=[1.2]; % Pipe Length ranges (m)
[Z,N]=ode15s('numdist',[0:L],No);
Why are using fixed numer 168? ode15s can return less number of points
n=N(length(N),:);
% some code
for j=1:168
sums(3:3:18)=sums(3:3:18)+beta(i,j)*n(j);
end
  5 Comments
Ben Ayres
Ben Ayres on 6 Apr 2021
However beta is dependent on i,j also and thus is currently returning complex 1x6 of:
1.98562768399673e-15 + 1.15452019649297e-15i 7.03050426650830e-15 + 5.34500090968969e-15i 5.20740451817902e-14 + 4.27600072775176e-14i 8.04944657622930e-13 + 6.68125113711212e-13i 6.43538727203316e-12 + 5.34500090968969e-12i 5.14789281873150e-11 + 4.27600072775175e-11i

Sign in to comment.

Categories

Find more on Programming in Help Center and File Exchange

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!