Problem with time plot from partial derivative equation solved by method of lines

1 次查看(过去 30 天)
Hello, I have a problem ,I tried to solve an pde equation using method of lines but appears a warning when time increases that dont let me raise the time and has a limit at time 7.5 h , dont know to solve it . I show the equations to better see the problem.Thanks a lot.
desorcion_final(0.4,0.5,2,0.02,20,1520,8,50.63,3.76,100,"si")
Warning: Failure at t=2.777457e+04. Unable to meet integration tolerances without reducing the step size below the smallest value allowed (5.820766e-11) at time t.
Warning: Imaginary parts of complex X and/or Y arguments ignored.
ans =
0.7441 - 0.0000i 7.5932 + 0.0000i
%tf means the time in hours
function [a]=desorcion_final(L,eps,u,k_f,c0,rhop,tf,Kf,n,q_inicial,sn)
Nt=60;
tf=tf*3600;
t=linspace(0,tf,Nt);
Nz=70;
z=linspace(0,L,Nz);
dz=z(2)-z(1);
q0star=Kf*c0^(1/n);
%Initial conditions
ICA=ones(1,Nz)*c0;
ICB=ones(1,Nz)*q_inicial;
IC=[ICA ICB];
[t, y]=ode15s(@fun_pde,t,IC,[],Nz,eps,n,Kf,k_f,u,rhop,dz);
%define value
cc=y(:,1:Nz);
qq=y(:,Nz+1:2*Nz);
%recalculate new limit conditions
cc(:,1)=0;
cc(:,end)=cc(:,end-1); %dC/dz=0
%plotting
cp=cc(:,end)./c0;
qp=qq(:,end-1)/q_inicial;
if sn=="si"
plot(t/3600,(1-qp))
title('Curve desorption conc')
xlabel('time(hora)')
ylabel('Convertion')
end
%%
%Return
a(1)=(1-qp(end));%conversion;
a(2)=t(end)/3600; %tiempo;
%%
% PDE function
function dydt=fun_pde(t,y,Nz,eps,n,KF,k_f,u,rhop,dz)
dcdt=zeros(Nz,1);
dqdt=zeros(Nz,1);
c=y(1:Nz);
q=y(Nz+1:2*Nz);
%BC
c(1)=0; %BC1
c(end)=c(end-1);%BC2
%interior
for i=2:Nz-1
qstar(i)=KF.*c(i).^(1/n);%Freundlich isotherm
dqdt(i)=k_f.*(qstar(i)-q(i));
dcdz(i)=(c(i+1)-c(i-1))./2./dz;
dcdt(i)=-u*dcdz(i)-rhop*((1-eps)./eps).*dqdt(i);
end
dydt=[dcdt;dqdt];
end
end

采纳的回答

Torsten
Torsten 2023-10-28
编辑:Torsten 2023-10-28
It's "Conversion", not "Convertion".
And setting a boundary condition at z = L is incorrect because you only have convection, no diffusion.
desorcion_final(0.4,0.5,2,0.02,20,1520,8,50.63,3.76,100,"si")
ans = 1×2
0.7497 8.0000
%tf means the time in hours
function [a]=desorcion_final(L,eps,u,k_f,c0,rhop,tf,Kf,n,q_inicial,sn)
Nt=60;
tf=tf*3600;
t=linspace(0,tf,Nt);
Nz=70;
z=linspace(0,L,Nz);
dz=z(2)-z(1);
q0star=Kf*c0^(1/n);
%Initial conditions
ICA=ones(1,Nz)*c0;
ICB=ones(1,Nz)*q_inicial;
IC=[ICA ICB];
[t, y]=ode15s(@fun_pde,t,IC,[],Nz,eps,n,Kf,k_f,u,rhop,dz);
%define value
cc=y(:,1:Nz);
qq=y(:,Nz+1:2*Nz);
%recalculate new limit conditions
cc(:,1)=0;
%cc(:,end)=cc(:,end-1); %dC/dz=0
%plotting
cp=cc(:,end)./c0;
qp=qq(:,end-1)/q_inicial;
if sn=="si"
plot(t/3600,(1-qp))
title('Curve desorption conc')
xlabel('time(hora)')
ylabel('Convertion')
end
%%
%Return
a(1)=(1-qp(end));%conversion;
a(2)=t(end)/3600; %tiempo;
%%
% PDE function
function dydt=fun_pde(t,y,Nz,eps,n,KF,k_f,u,rhop,dz)
dcdt=zeros(Nz,1);
dqdt=zeros(Nz,1);
c=y(1:Nz);
q=y(Nz+1:2*Nz);
%BC
c(1)=0; %BC1
%c(end)=c(end-1);%BC2
%interior
for i=2:Nz
qstar(i)=KF.*c(i).^(1/n);%Freundlich isotherm
dqdt(i)=k_f.*(qstar(i)-q(i));
if i < Nz
dcdz(i)=(c(i+1)-c(i-1))/2/dz;
else
dcdz(i)=(c(i)-c(i-1))/dz;
end
dcdt(i)=-u*dcdz(i)-rhop*((1-eps)./eps).*dqdt(i);
end
dydt=[dcdt;dqdt];
end
end

更多回答(0 个)

产品


版本

R2023a

Community Treasure Hunt

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

Start Hunting!

Translated by