Pdepe solver for diffusion equation in a packed column

[SimOpera]

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 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 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.
system_matrix=xlsread('mcl system physical constants.xlsx');

% 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. Thanks for your time.
-SimOpera

 

[mmwave]

Dear SimOpera,

Thank you for sharing your work with us. It is always exciting to see scientists working on complex problems and trying to find solutions. Your program seems to be very well-organized and comprehensive. However, I am not familiar with the specific equations and parameters you are using, so I am unable to provide specific feedback on the accuracy of your results.

However, here are some general suggestions that may help you improve your program:

1. Ensure that your input parameters are correct: Before running your program, it is important to double check that all the input parameters (such as velocity, length, diameter, porosity, etc.) are correct and accurately represent the system you are trying to model. Even small errors in these values can lead to significant discrepancies in your results.

2. Use appropriate initial and boundary conditions: The accuracy of your results also depends on the initial and boundary conditions that you have defined. Make sure that they are consistent with the system you are trying to model and that they are applied correctly in your program.

3. Check for errors in your equations: It is always a good idea to double check your equations to make sure they are correct. Even a small error in one of the equations can lead to significant errors in your results. You can try solving simpler cases (such as a single component system) and comparing your results with analytical solutions to verify the accuracy of your equations.

4. Use appropriate mesh sizes: The accuracy of your results also depends on the mesh sizes you have chosen for space and time. Make sure that they are small enough to capture the behavior of your system, but not too small that they make your program run very slowly.

5. Visualize your results: It is always helpful to visualize your results to get a better understanding of the behavior of your system. You can plot the concentration profiles for each component at different times to see how they change over time. You can also compare your results with published data or analytical solutions to verify their accuracy.

I hope these suggestions help you improve your program and get more accurate results. Keep up the good work and don't hesitate to reach out to other scientists for help and feedback. Good luck!