Problem with time plot from partial derivative equation solved by method of lines
1 visualizzazione (ultimi 30 giorni)
Mostra commenti meno recenti
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.
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1523476/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1523481/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1523486/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1523491/image.png)
desorcion_final(0.4,0.5,2,0.02,20,1520,8,50.63,3.76,100,"si")
%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
0 Commenti
Risposta accettata
Torsten
il 28 Ott 2023
Modificato: Torsten
il 28 Ott 2023
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")
%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 Commenti
Più risposte (0)
Vedere anche
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!