please help plot and visualize the results

Question

Ofentse Maila 2023-7-16

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1996543-please-help-plot-and-visualize-the-results

编辑： dpb 2023-7-16

clear 
clc
close all
%input data
c0=2;
L=0.5;
eps=0.62;
u=0.06;
De=9.53e-4;
rhop=1940;
kl=0.090;
qm=200;
ks=0.00425;
%Length
Nz=201;
z=linspace(0,L,Nz);
dz=z(2)-z(1);
%time
t=[0 1800];
%initial conditions
ICA=zeros(1,Nz); %for PDE dc/dt...
ICB=zeros(1,Nz); %for ODE dq/dt...
IC=[ICA ICB];
%odesolver
[t, y]=ode15s(@f,t,IC,[],eps,u,kl,c0,rhop,Nz,dz,De,qm,ks);
%define value or recalculation using coding function
%define value
%dont forget :
c=y(:,1:Nz);
qstar=y(:,Nz+1:2*Nz);
%boundary conditions
c(:,1)=c0;
c(:,end)=(4*c(:,end-1)-c(:,end-2))./3; % using backward second order
%function
function dydt=f(t,y,a,eps,u,kl,c0,rhop,Nz,dz,De,qm,ks)
dydt=zeros(length(y),1);
dcdt=zeros(Nz,1);
dqdt=zeros(Nz,1);
dcdz=zeros(Nz,1);
dc2dz2=zeros(Nz,1);
%define value
c=y(1:Nz);
qstar=y(Nz+1:2*Nz);
%boundary conditions
c(1)=c0;
c(end)=(4*c(end-1)-c(end-2))./3; % using backward second order
%interior
for i=2:Nz-1
    qstar(i)=qm.*ks.*c(i)./((1+(ks.*c(i))^1/a)); %Toth model
    dqdt(i)=kl.*(qstar(i)-q(i));
    dcdz(i)=(c(i+1)-c(i-1))./2./dz; %centered finite diff
    dc2dz2(i)=(c(i+1)-2*c(i)+c(i-1))./dz.^2; %centred second derivative
    dcdt(i)=De.*dc2dz2(i)-u*dcdz(i)-rhop.*((1-eps)./eps).*dqdt(i); %trabsport equation
end
dydt=[dcdt;dqdt];
end