Solving partial differential equations using pdepe solver

2022 6 9
1 个回答
回答已采纳
更新时间：2023 5 1
4 次查看（30 天）

0 个投票

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.

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

 采纳的回答

Torsten
Torsten 2022-6-9
编辑：Torsten 2022-6-9

0 个投票

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
Ari Dillep Sai 2022-6-10
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
Ari Dillep Sai 2022-6-10
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
Torsten 2022-6-10
编辑：Torsten 2022-6-10
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
Ari Dillep Sai 2022-6-10
Ok. I will look into the parameters.
Torsten
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
Rahmat Sunarya 2023-5-1
Ran in:
Ran in:
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

请先登录，再进行评论。

更多回答（0 个）

类别

在 帮助中心 和 File Exchange 中查找有关 Ordinary Differential Equations 的更多信息

产品

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by