Solving Population Balance Equation for multiple initial particle dimensions

5 次查看（过去 30 天）
Ben Ayres
Ben Ayres 2021-4-6
评论： KAJAL KUMARI 2021-6-25
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 个评论
显示 2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

回答（1 个）

darova
darova 2021-4-6
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 个评论
显示 3更早的评论
Ben Ayres
Ben Ayres 2021-4-6
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

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Particle & Nuclear Physics 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by