# Pdepe solver for diffusion equation in a packed column

• MATLAB

## Main Question or Discussion Point

Consider I have a packed column of length L filled with known characteristic adsorbent. I am putting a mixture of N components in it and I am solving for concentration of each component in mobile phase at the outlet of the column. The equations which are to be generalised are as follows: An example of the 2 component systems is stated in the paper below which can help you to understand the equations better if following is unclear.http://www.sciencedirect.com/science/article/pii/S1369703X10000069?np=y

for component i:N:
Dcoe*∂2Ci/∂x2−v*∂Ci/∂x−∂qi/∂t=∂Ci/∂t
∂q1∂t=ki[Ciqmi(1−Σ(qi/qmi))−kd1q1]

[PLAIN]http://www.sciencedirect.com/sd/blank.gifwhere [Broken] the initial conditions are
t=0,x=0,Ci(0,0)=Cfi,qi(0,0)=0,
http://www.sciencedirect.com/sd/blank.gif
t=0,x≠0,Ci(x,0)=0,qi(x,0)=0.

[PLAIN]http://www.sciencedirect.com/sd/blank.gifand [Broken] the boundary conditions are
x=0,Ci(0,t)=Cfi,qi(0,t)=0
http://www.sciencedirect.com/sd/blank.gif
x=L, ∂Ci/∂x=0,∂qi/∂x=0.

I am solving these equations by using pdepe solver.

Following 4 files are used to solve pdepe: mcl_brkthrupdemain.m (pdepe main solution function), mcl_brkthrupde (stating equations), mcl_brkthruic.m (initial conditions), mcl_brkthrubc(boundary conditions)
%% file mcl_brkthrupdemain is used to define the input parameters and the pdepe main functions for solving the equations.
Code (Text):

function mcl_brkthrupdemain
mcl_globalinc;
m=0;
%% Importing system parameters from an excel sheets.

% input('enter the value for velocity V in cm/min:');
vel = system_matrix(1,1);

% input('enter the value for Length of the column in cm:');
L= system_matrix(1,2);

% input('enter the value for column diameter in cm:');
% Not used currently
dc= system_matrix(1,3);

% input('enter the value of porosity:');
epsilon= system_matrix(1,4);

%% importing parameters for component specific constants from excel sheets
matrix_in = xlsread('mcl Excelimport.xlsx'); % Main input matrix file from user
[nComp,nCin] = size(matrix_in); % defining number of components

k1 = zeros(nComp,1);
kd = zeros(nComp,1);
qm = zeros(nComp,1);
Cin = zeros(nComp,1);
Dcoe = zeros(nComp,1); % input('enter the value for diffusion coefficient Dcoe in cm^2/min:');

for w=1:nComp
k1(w) = matrix_in(w,1);
kd(w) = matrix_in(w,2);
qm(w) = matrix_in(w,3);
Cin(w) = matrix_in(w,4);
Dcoe(w) = matrix_in(w,5);
end
%% Defining mesh sizes for both space and time.
N=50;
Run_time=800;
% N =input('enter the value for number of meshing intervals:\n');
% Run_time= input('enter the value for maximum time of run:\n');

x=linspace(0,L,N); % spatial meshing
t=linspace(0,Run_time,N); % time meshing
tend=300;
%% solution and plotting
sol=pdepe(m,@mcl_brkthrupde,@mcl_brkthruic,@mcl_brkthrubc,x,t);
colorvec= hsv(5);
for g=1:nComp
C=sol(:,:,2*g-1);
solid_conc=sol(:,:,2*g);
Cfr=C/matrix_in(g,4);
disp(C);
disp(solid_conc);
% surface plots for conentrations
%surf(x,t,C);
%title(sprintf('surface plot for concentration of component %d',g))

% A solution profile for conc. vs length for all components.
figure, plot(x,C(end,:), 'color',colorvec(g, title(sprintf('solution at runtime end for Conc %d vs distance travelled',g))
xlabel('distance (cm)')
ylabel(sprintf ('component %d concentration',g))
hold on

%Plot conc. vs. time for component 1
figure, plot(t,Cfr)
title(sprintf('breakthrough for component %d at all x',g))
xlabel('Time (min)')
ylabel(sprintf('Component %d Concentration (mg/ml)',g))

% plot end concentration with time
figure, plot(t,Cfr(:,end))
title(sprintf('breakthrough for Component %d',g))
xlabel('Time (min)')
ylabel(sprintf('Concentration of component%d (mg/ml)',g))
end
end

%% file mcl_brkthrupde is used to define equations.
Code (Text):
function [c,f,s]= mcl_brkthrupde(x,t,u,DuDx) % inputs for equation coefficients
%% Redefining global variables
%% Defining system constants
mcl_globalinc;
neq=2*nComp;
ddf=x;
%% Usual variables
conc =zeros(nComp,1);
q =zeros(nComp,1);
dcdx =zeros(nComp,1);
dqdx =zeros(nComp,1);
dqdt =zeros(nComp,1);

for w=1:nComp
conc(w,1) = u(2*w-1);
q(w,1) = u(2*w);
dcdx(w,1) = DuDx(2*w-1);
dqdx(w,1) = DuDx(2*w);
end
%% writing c term and putting it into a column vector
c=ones(neq,1);
%disp(c);
%% writing f term and putting it into a column vector
f=zeros(neq,1);
count1=1;
while count1<(neq) % f term for all the equations
for p=1:nComp
f(count1,1)=Dcoe(p).*dcdx(p,1);
f(count1+1,1)=0.*dqdx(p,1);
count1=count1+2;
end
end
disp(f);
%% finding intermediate summation for q/qmax
s=zeros(neq,1);
SUM=0;
for w=1:nComp
SUM=SUM+q(w,1)/qm(w);
end
disp(SUM);
for r=1:nComp
dqdt(r,1)=k1(r)*conc(r,1)*qm(r)*(1-SUM)-k1(r)*kd(r).*q(r,1);
end
count2=1;
while count2<(neq) % f term for all the equations
for b=1:nComp
s(count2,1)=-vel*dcdx(b,1)-((1-epsilon)/epsilon).*dqdt(b,1);
s(count2+1,1)=dqdt(b,1);
count2=count2+2;
end
end
disp(s);
end
%% mcl_brkthruic describes initial conditions
Code (Text):
function u0= mcl_brkthruic(x) %function for inlet conditions
mcl_globalinc;
%% defining initial conditions for concentration for x=0 and 0<x<L
u0=zeros(2*nComp,1);
if x==0
for i=1:nComp
u0(2*i-1,1)= Cin(i);
u0(2*i,1)= 0;
end
end
disp(u0);
end

%% mcl_brkthrubc describes boundary conditions

Code (Text):
function [pl,ql,pr,qr]= mcl_brkthrubc(xl,ul,xr,ur,t) % function for boundary conditions
%% calling required variable set
mcl_globalinc;

neq=2*nComp;
%% defining usual variables
pl=zeros(neq,1);
ql=zeros(neq,1);
pr=zeros(neq,1);
qr=zeros(neq,1);

cleft=zeros(nComp);
qleft=zeros(nComp);
cright=zeros(nComp);
qright=zeros(nComp);

%% finding pl (left boundary) for terms without differential term and putting into a column vector
for i=1:nComp
cleft(i)=ul(2*i-1);
qleft(i)=ul(2*i);
cright(i)= ur(2*i-1);
qright(i)=ur(2*i);
end
for j=1:nComp;
pl(2*j-1,1)=cleft(j)-Cin(j); %left for conc in liquid phase
pl(2*j,1)=qleft(j); % left for conc in solid phase
ql(2*j-1,1)=0; % left for differential term for conc in liquid phase
ql(2*j,1)=0; % left for differential term for conc in solid phase

pr(2*j-1,1)=cright(j); %right for conc in liquid phase
pr(2*j,1)=qright(j); % right for conc in solid phase
qr(2*j-1,1)=1; % right for differential term for conc in liquid phase
qr(2*j,1)=1; % right for differential term for conc in solid phase
end
disp(pl);
disp(ql);
disp(pr);
disp(qr);
end

As you can see that the paper mentions a bi component system and I am trying to make it a multicomponent system.
This program is taking values from input xlsx files "mcl_system physical constants.xlsx" and "mcl Excelimport.xlsx". You can take the values you like to start with.
It is giving me the profile but not the correct profile.
Can anyone please have a look?
Let me know if anything is unclear.

-SimOpera

Last edited by a moderator: