Solving partial differential equations using pdepe solver

I have to solve two coupled partial differential equations which I have attached below. I used pdepe solver for those equations but I am getting error. I would like to know whether I can solve such equations with pdepe solver.

I have written the following code for the above equations.

% paramters
global P R T L epsf pel tc epstar beta Keq cref yin sh gamma p

% paramters
L = 6e-2;                           % in m
epsf = 0.43;    
beta = 46.65;
Keq = 75;
epstar = epsf + Keq*beta;
pel = 21.73;
tc = 140e-3;                        % in sec
P = 1e5;                            % in pascals
R = 8.314;                          % in j/mol-k
T = 323;                            % in K
cref = P/(R*T);                     % in mol/m3
yin = 2420;
sh = 4.51;
gamma = 0.111;
p = 0.0326;

% defining mesh
xmesh = linspace(0,1,300);
tmesh = linspace(0,1500,500);

% using pdepe solver
m = 0;
sol = pdepe(m,@adsorp,@adsorpic,@adsorpbc,xmesh,tmesh);

% % plotting results
% plot(tmesh,sol(:,end)*1e6*R*T/P)
% xlabel('time (sec)')
% ylabel('concentration (ppm)')

function [c,f,s] = adsorp(x,t,u,dudx)
global epsf beta pel sh p gamma Keq

c = [epsf; beta];
f = [((1/pel)*dudx(1));0 ];
s = [-dudx(1)-(sh/p)*(u(1)-(gamma*u(1)+beta*u(2)/Keq)/(gamma+beta*(1-u(2)))); (sh/p)*(u(1)-(gamma*u(1)+beta*u(2)/Keq)/(gamma+beta*(1-u(2))))];
end

function u0 = adsorpic(x)
u0 = [0;0];
end

function [pl,ql,pr,qr] = adsorpbc(xl,ul,xr,ur,t)
global yin

pl = [ul(1)-yin;0];
ql = [1;0];
pr = [0;0];
qr = [1;0];
end

when I am running above code I am getting error as follows:

Error using pdepe (line 293)
Spatial discretization has failed. Discretization supports only parabolic and elliptic equations, with flux term
involving spatial derivative.

Error in systemofpdesusingpdepe (line 27)
sol = pdepe(m,@adsorp,@adsorpic,@adsorpbc,xmesh,tmesh);

Can someone please tell whether the pdepe solver can be used for these type of problems? Also suggestions on the method to solve above equations is very much appreciated.

Answer 1

[Torsten]

Some people say, one can take pdepe for this kind of problem (e.g. Bill Greene), I'd say one shouldn't because it's problematic.

The second equation is a pure ordinary differential equation without spatial derivatives. Thus boundary conditions for it cannot be imposed. And you did this: you didn't impose boundary conditions. But pdepe expects them at both ends. So one way to artificially deal with this problem is to set pl(2) = pr(2) = 0, ql(2) = qr(2) = 1. Try if the results come out to be reasonable.

And I think your boundary condition for y_f in the code has the wrong sign (at least the sign compared to your attached image is reversed).

Why do you try pdepe ? Your code using the method of lines was ok in my opinion. It's not a problem with the code, but with the adsorption parameters that caused that the breakthrough was too early.

[Ari Dillep Sai]

Hi Torsten,

Thanks for your reply.

Why do you try pdepe ? Your code using the method of lines was ok in my opinion. It's not a problem with the code, but with the adsorption parameters that caused that the breakthrough was too early.

I have used the method of lines code for this. I am attaching the code below. But when I am running the code, I am getting concentration same as initial concentration in the complete bed.

clc
clear
close all

% parameters
L=6e-2;
epsf=0.43;
Keq=74;
pel = 21.73;
Def=11.8e-4;
P = 1e5;
R = 8.314;
T = 50+273;
yin = 2420;
Cin=(yin*1e-6*P)/(R*T);
uf=0.43;
sh =4.51;
p = 0.0326;
gamma = 0.111;
beta = 46.65;


% creating step size
Nz=301;
z=linspace(0,L,Nz);
dz=z(2)-z(1);


% creating time steps 
t=linspace(0,15000,500);

% initial Conditions
IC_A=zeros(1,Nz); %for C
IC_B=zeros(1,Nz); %for q
IC=[IC_A, IC_B];

% using ode solver
[t,y]=ode15s(@f12,t,IC);

% plotting results
plot(t,y(:,end)*1e6*R*T/P)

function dydt=f12(~,y)
L=6e-2;
epsf=0.43;
Keq=74;
pel = 21.73;
Def=11.8e-4;
P = 1e5;
R = 8.314;
T = 50+273;
yin = 2420;
Cin=(yin*1e-6*P)/(R*T);
uf=0.43;
sh =4.51;
p = 0.0326;
gamma = 0.111;
beta = 46.65;

Nz=301;
z=linspace(0,L,Nz);
dz=z(2)-z(1);

dCdt=zeros(Nz,1);
dqdt=zeros(Nz,1);

C=y(1:Nz); %for C
q=y(Nz+1:2*Nz); %for q

% B.C at x = 0
C(1) = Cin;

% creating interior points
for i=2:Nz-1
    dqdt(i)=(sh/(p*beta))*(C(i)-(gamma*C(i)+beta*q(i)/Keq)/(gamma+beta*(1-q(i))));
    dCdt(i)=(-1/(2*dz*epsf))*(C(i+1)-C(i-1)) + (1/(dz^2*epsf*pel))*(C(i+1)-2*C(i)+C(i-1)) - (beta/epsf)*dqdt(i);
end
dCdt(Nz) = dCdt(Nz-1);

dydt=[dCdt;dqdt];

end

I am not able to figure out what is going wrong in the above code. So I tried solving with the pdepe solver.

If you are able to figure out any mistakes in the above code, please let me know.

Also in the above code for simplicity I have taken boundary condition at x = 0 as C(x=0) = Cin.

[Torsten]

But Davides code worked:

https://de.mathworks.com/matlabcentral/answers/1730075-how-to-solve-coupled-partial-differential-equations-with-method-of-lines?s_tid=prof_contriblnk

Anything wrong with it ?

[Ari Dillep Sai]

I tried with that code as well.

But I am not getting the correct results with that as well. I am attaching that code below.

% Parameters
p.epsf         = 0.43;                           % bed porosity  
p.beta         = 46.65;                           % superficial velocity in m/s  
p.L            = 1;                          % bed length
p.Keq          = 74;
p.epstar      = p.epsf + p.Keq*p.beta;
p.pel = 21.73;
p.sh = 4.51;
p.pe = 0.0326;
p.gamma = 0.111;
p.tc = 140e-3;                        % in sec
p.P = 1e5;                            % in pascals
p.R = 8.314;                          % in j/mol-k
p.T = 323;                            % in K
p.cref = p.P/(p.R*p.T);                     % in mol/m3
p.yin = 2420;

% using method of lines
tf      = 15000;                               % End time
p.nz    = 400;                                  % Number of nodes
z       = linspace(0,p.L,p.nz);                 % Space domain
dz      = diff(z); p.dz = dz(1);                % Space increment


% Initial conditions
c       = zeros(p.nz,1);
c(1)    = p.yin;
q       = zeros(p.nz,1);
y0      = [c;q];   
tspan   = linspace(0,tf,500);

% Mass matrix (for DAE system)
M               = eye(2*p.nz);
M(1,1)          = 0;
M(p.nz,p.nz)    = 0;
options         = odeset('Mass',M,'MassSingular','yes');

[t,y]           = ode15s(@(t,y)rates(t,y,p),tspan,y0,options);

% Plot results
plot(t,y(:,p.nz))
grid

function out = rates(~,y,p)
% Variables
c       = y(1:p.nz,1);
q       = y(p.nz+1:2*p.nz,1);

% Space derivatives
dcdz            = zeros(p.nz,1);
d2cdz2          = zeros(p.nz,1);
dcdz(2:end-1)   = (c(3:end)-c(1:end-2))/(2*p.dz);
d2cdz2(2:end-1) = (c(3:end)-2*c(2:end-1)+c(1:end-2))/(p.dz^2);

% equilibrium concentration
cs = (p.gamma*c+p.beta*q/p.Keq)./(p.gamma+p.beta*(1-q));

% Governing equations
dqdt = (p.sh/(p.pe*p.beta))*(c-cs); 
dcdt = (1/(p.pel*p.epsf))*d2cdz2-(1/p.epsf)*dcdz-(p.beta/p.epsf)*dqdt;
       
% Boundary conditions
dcdt(1)     = (-1/p.pel)*(c(2)-c(1))/p.dz+(c(1)-p.yin);
dcdt(end)   = c(end)-c(end-1);

% Output
out = [dcdt;dqdt];

end

[Torsten]

As I already said: the reason that you don't get the "correct" results is not your code, but your parameters in the adsorption model or the equilibrium model itself.

[Ari Dillep Sai]

Ok. I will look into the parameters.

[Torsten]

E.g. if you multiply p.K0 by a factor of 8, you get a breakthrough time at about 300 minutes.

[Rahmat Sunarya]

I see my style of coding here similar to inscript below :)

please visit my channel for more: https://www.youtube.com/channel/UCirsOp2CQzbHOseEvj9qcsw

Thank You

clc
clear
close all

% parameters
L=6e-2;
epsf=0.43;
Keq=74;
pel = 21.73;
Def=11.8e-4;
P = 1e5;
R = 8.314;
T = 50+273;
yin = 2420;
Cin=(yin*1e-6*P)/(R*T);
uf=0.43;
sh =4.51;
p = 0.0326;
gamma = 0.111;
beta = 46.65;


% creating step size
Nz=301;
z=linspace(0,L,Nz);
dz=z(2)-z(1);


% creating time steps 
t=linspace(0,15000,500);

% initial Conditions
IC_A=zeros(1,Nz); %for C
IC_B=zeros(1,Nz); %for q
IC=[IC_A, IC_B];

% using ode solver
[t,y]=ode15s(@f12,t,IC);

% plotting results
plot(t,y(:,end)*1e6*R*T/P)

function dydt=f12(~,y)
L=6e-2;
epsf=0.43;
Keq=74;
pel = 21.73;
Def=11.8e-4;
P = 1e5;
R = 8.314;
T = 50+273;
yin = 2420;
Cin=(yin*1e-6*P)/(R*T);
uf=0.43;
sh =4.51;
p = 0.0326;
gamma = 0.111;
beta = 46.65;

Nz=301;
z=linspace(0,L,Nz);
dz=z(2)-z(1);

dCdt=zeros(Nz,1);
dqdt=zeros(Nz,1);

C=y(1:Nz); %for C
q=y(Nz+1:2*Nz); %for q

% B.C at x = 0
C(1) = Cin;

% creating interior points
for i=2:Nz-1
    dqdt(i)=(sh/(p*beta))*(C(i)-(gamma*C(i)+beta*q(i)/Keq)/(gamma+beta*(1-q(i))));
    dCdt(i)=(-1/(2*dz*epsf))*(C(i+1)-C(i-1)) + (1/(dz^2*epsf*pel))*(C(i+1)-2*C(i)+C(i-1)) - (beta/epsf)*dqdt(i);
end
dCdt(Nz) = dCdt(Nz-1);

dydt=[dCdt;dqdt];

end