please help plot and visualize the results

3 次查看（过去 30 天）
Ofentse Maila
Ofentse Maila 2023-7-16
编辑： 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

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

回答（1 个）

Sahaj
Sahaj 2023-7-16
编辑：Sahaj 2023-7-16
Hi Ofentse.
You can add the following code after the ODE solver section to plot and visualiize your results.
% Plot concentration profile
subplot(2,1,1);
plot(z, c(end,:), 'b-', 'LineWidth', 2);
xlabel('Distance (z)');
ylabel('Concentration (c)');
title('Concentration Profile');
% Plot deposition rate profile
subplot(2,1,2);
plot(z, qstar(end,:), 'r-', 'LineWidth', 2);
xlabel('Distance (z)');
ylabel('Deposition Rate (q*)');
title('Deposition Rate Profile');

类别

在 Help Center 和 File Exchange 中查找有关 Geometry and Mesh 的更多信息

标签

产品


版本

R2023a

Community Treasure Hunt

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

Start Hunting!

Translated by