cFeed = 10; % Feed concentration
L = 0.01; % Column length
t0 = 0; % Initial Time
tf = 90; % Final time %%%%%%%%%%% Shorter time
dt = 1; % Time step %%%%%%%%%%% Smaller timestep
z = 0:0.001:L; % Mesh generation %%%%%%%%%%% Finer mesh
t = t0:dt:tf;% Time vector
n = numel(z); % Size of mesh grid
%Initial Conditions / Vector Creation
c0 = zeros(n,1);% t = 0, c = 0
c0(1)=cFeed;q0 = zeros(n,1); % t = 0, q = 0 for all z,
q0(1) = zeros; y0 = [c0 ; q0]; % Appends conditions together
%ODE15S Solver
[T, Y] = ode15s(@(t,y) Myfun(t,y,z,n),[t0 tf],y0);plot(T, Y(:,n)/cFeed,'b')xlabel('time')ylabel('exit conc.')Q = trapz(T,Y(:,n));Efficiency_ads = Q/(tf*cFeed);function DyDt=Myfun(~,y,z,n)% Defining Constants
D = 8.76e-4; % Axial Dispersion coefficient
v = 0.0574; % Superficial velocity
epsilon = 0.383; % Voidage fraction
Kldf = 1.38e-3; % Mass Transfer Coefficient
n = 0.5396; % Langmuir parameter
qm = 28.75;% Saturation capacity
b = 2.052e-3;rho = 480;%Tc = zeros(n,1);
%q = zeros(n,1);
DcDt = zeros(n,1);DqDt = zeros(n,1);%DyDt = zeros(2*n,1);
zhalf = zeros(n-1,1);DcDz = zeros(n,1);D2cDz2 = zeros(n,1);c = y(1:n);q = y(n+1:2*n);% Interior mesh points
zhalf(1:n-1)=(z(1:n-1)+z(2:n))/2;for i=2:n-1 DcDz(i) = ((z(i)-z(i-1))/(z(i+1)-z(i))*(c(i+1)-c(i))+(z(i+1)-z(i))/(z(i)-z(i-1))*(c(i)-c(i-1)))/(z(i+1)-z(i-1)); D2cDz2(i) = (zhalf(i)*(c(i+1)-c(i))/(z(i+1)-z(i))-zhalf(i-1)*(c(i)-c(i-1))/(z(i)-z(i-1)))/(zhalf(i)-zhalf(i-1));end% Calculate DcDz and D2cDz2 at z=L for boundary condition dc/dz = 0
DcDz(n) = 0; D2cDz2(n) = -1.0/(z(n)-zhalf(n-1))*(c(n)-c(n-1))/(z(n)-z(n-1));% Set time derivatives at z=0
% DcDt = 0 since c=cFeed remains constant over time and is already set as initial condition
% Standard setting for q
DqDt(1) = Kldf*((qm*b*c^1/n(1))/(1+(b*c^1/n(1))) - q(1));DcDt(1) = 0.0;% Set time derivatives in remaining points
for i=2:n %Equation: dq/dt = k(q*-q) where q* = qs * (b*c)/(1+b*c)
DqDt(i) = Kldf*((qm*b*c^1/n(i))/(1+(b*c^1/n(i))) - q(i)); %Equation: dc/dt = D * d2c/dz2 - v*dc/dz - ((1-e)/(e))*dq/dt
DcDt(i) = D*D2cDz2(i) - v*DcDz(i) - ((1-epsilon)/(epsilon))*rho*DqDt(i); end% Concatenate vector of time derivatives
DyDt = [DcDt;DqDt];end
MATLAB: When I attempt to run this it produces an ‘error using zeros, size inputs must be integers’ for the line DcDt = zeros(n,1). The code is supposed to produce a breakthrough curve. Any help is appreciated.
adsorption columnerror with zeros
Best Answer