Solving Population Balance Equation for multiple initial particle dimensions

9 visualizaciones (últimos 30 días)
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 comentarios
Ben Ayres
Ben Ayres el 7 de Abr. de 2021
Thanks Darova, these probes helped me solve the issue
KAJAL KUMARI
KAJAL KUMARI el 25 de Jun. de 2021
can you provide me the whole code? I am also working in a similar project.

Iniciar sesión para comentar.

Respuestas (1)

darova
darova el 6 de Abr. de 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 comentarios
darova
darova el 6 de Abr. de 2021
But you are trying to access beta(168,168). It doesn't exist
Ben Ayres
Ben Ayres el 6 de Abr. de 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

Iniciar sesión para comentar.

Categorías

Más información sobre Programming en Help Center y File Exchange.

Community Treasure Hunt

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

Start Hunting!

Translated by