Solving partial differential equations using pdepe solver

Question

Ari Dillep Sai 2022-6-9

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1737425-solving-partial-differential-equations-using-pdepe-solver

评论： Rahmat Sunarya 2023-5-1

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.

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Torsten 2022-6-9

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1737425-solving-partial-differential-equations-using-pdepe-solver#answer_982565

编辑：Torsten 2022-6-9

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.

7 个评论
显示 5更早的评论隐藏 5更早的评论

Ari Dillep Sai 2022-6-10

在 MATLAB Online 中打开

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.

Ari Dillep Sai 2022-6-10

在 MATLAB Online 中打开

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 2022-6-10

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

Rahmat Sunarya 2023-5-1

在 MATLAB Online 中打开

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

请先登录，再进行评论。

Solving partial differential equations using pdepe solver

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

采纳的回答

7 个评论
显示 5更早的评论隐藏 5更早的评论

更多回答（0 个）

另请参阅

类别

标签

产品

Community Treasure Hunt

Solving partial differential equations using pdepe solver

0 个评论 显示 -2更早的评论隐藏 -2更早的评论

采纳的回答

7 个评论 显示 5更早的评论隐藏 5更早的评论

更多回答（0 个）

另请参阅

类别

标签

产品

Community Treasure Hunt

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

7 个评论
显示 5更早的评论隐藏 5更早的评论