Solving Population Balance Equation for multiple initial particle dimensions

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

[darova]

Can you please show the line the erorr occurs?

[Ben Ayres]

The error occurs in the for loop:

Line - sums (3:3:18) +beta(i,j)*n(j)

However it seems like it is all of the for loop equations under ‘sums’

Thanks

[Ben Ayres]

Thanks Darova, these probes helped me solve the issue

[KAJAL KUMARI]

can you provide me the whole code? I am also working in a similar project.

Answer 1

[darova]

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

[Ben Ayres]

I have since changed my span to 0:0.01:1.2

I used 168 instead of writing SECT*6 so that all of the values were covered in the 28x6 matrix

[darova]

What is size of beta?

[Ben Ayres]

currently a 1x6 for each particle diameter

[darova]

But you are trying to access beta(168,168). It doesn't exist

[Ben Ayres]

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