%I can't find where I've made a mistake in the code modelling advection diffusion in 1d. I'm trying to get 2 radial plots that show diffusion of the concentrations at t=0 and t=end but i get the exceed array bounds error. I have pasted the exact code below:
clf
clear
clc
% Inputs
D=1.04e-6; % Thermal diffusivity, m2/s, eg 0.01
r_in=0;
r_out=0.05;
m=35; %number of divided sections between r_out and r_in
delta_r=(r_out-r_in)/m; %section length, m
nr=m+1; %total no. of radial points
nphi=36;%no. of angle steps
delta_phi=2*pi/nphi; %angle step
t=40000; % total time, s
u=delta_r/t; % veloctiy in x direction, m/s
nt=450000; % total no. of time steps
delta_t=t/nt; % timestep, s
Crin=0; % Concentration at r_in, ppm
Crout=1000; % Concentration at r_out, ppm
Cin=0; %initial concentration at r_in<r<r_out, ppm
Cmax=max([Crin,Crout]);
% Solution
for i=1:nr
r(i)=r_in+(i-1)*delta_r;
end
r;
disp(u)
phideg=[0:360/nphi:360+360/nphi];
phi=phideg.*pi/180;
c=u*delta_t/delta_r; % Courant Number
d=D*delta_t/delta_r^2; % diffusion number
d1=D*delta_t/(2*delta_r);
Pe=u*r/D; % Peclet Number
%fprintf('\n')
if d<=0.5
fprintf('solution stable\nd=%10.7f\n', d)
else
fprintf('solution unstable\nd=%10.7f\n', d)
end
fprintf('\n')
fprintf('\n')
if c^2<=2*d
fprintf('solution stable\n c^2 (%4.2f)<=2*d(%4.2f)\n', c^2, 2*d)
else
fprintf('solution unstable\n c^2 (%4.2f)>(%4.2f)\n', c^2, 2*d)
end
% Solution
for i=1:nr
r(i)=r_in+(i-1)*delta_r;
end
% Creating initial boundary conditions
C=zeros(nr,nphi+2,nt);
for k = 1:nt+1
for j = 1:nphi+2
for i=1:nr
if(i==1)
C(i,j,k)=Crin;
elseif(i==nr)
C(i,j,k)=Crout;
else
C(i,j,k)=Cin;
end
end
end
end
C;
% Creating C matrix
for k=1:nt
for j=1:nphi+2
for i= 2:nr-1
C(i,j,k+1)= C(i,j,k) +(c/2)*(C(i+1,j,k)-C(i-1,j,k)) + d*(C(i+1,j,k) - 2*C(i,j,k) + C(i-1,j,k))+(d1/r(i))*(C(i+1,j,k)-C(i-1,j,k));
end
end
end
C(:,:,1);
C(:,:,nt+1);
C;
C(:,1,1);
C(:,1,nt+1);
C(:,1,:);
%Initial Temperature profile
C1=C(:,:,1);
subplot(2,2,1)
r=r;
[Phi,R]=meshgrid(phi,r);
Z=C1; %C1.*Phi./Phi;
[X,Y,Z]=pol2cart(Phi,R,Z);
%figure(1)
[C,h]=contourf(X,Y,Z,300);
set(h,'LineColor','none')
grid on
axis equal
colorbar
title({'Initial Concentration Profile - 1D - Circular coordinates';'C in ppm'})
xlabel(['r in m'])
%final C
C2=C(:,:,nt+1); <----------------------------------------------------------------------- Where the error appears
subplot(2,2,3)
r=r;
[Phi,R]=meshgrid(phi,r);
Z=C2.*Phi./Phi;
[X,Y,Z]=pol2cart(Phi,R,Z);
%figure(1)
[C,h]=contourf(X,Y,Z,300);
set(h,'LineColor','none')
grid on
axis equal
colorbar
title({'Initial Concentration Profile - 1D - Circular coordinates';'C in ppm'})
xlabel(['r in m'])
