MATLAB: Problem: 3D Diffusion Equation with Sinkterm

3d3d plotsconcentrationdiffusiondiffusion equationdiffusion visualizationMATLABmeshgridpdeproblemprofile

Hi guys, I am working on a 3d simulation which shows the concentration profile in a 1m^3 box. In this box I placed a filter which filters out a concentration of substance X. The initial conditions for this problem are: At r (radius) = R (radius filter), c = cs (concentration = surface concentration of the filter), at r = inf, c = c0 (concentration = initial concentration), and for t=0, c =c0. The problem I encounter is that the concentration profile in the box changes too quickly to be believable. Also changing the diffusion coefficient has no result on the outcome of the simulation.
clear all
close all
%paramters
diffusion  = 4.058*10^-5; %calculated for substance X
%dimensions
Lx = 10;
Ly = 10;
Lz = 10;
Nx = 21; Nt  = 400000; %amount of iterations
Ny = 21;
Nz = 21;
dx = Lx/(Nx-1);
dy = Ly/(Ny-1);
dz = Lz/(Nz-1);
%CFL conditions, determines how big dt should be for the equation to
%converge
c = 1;
C=0.075; %C<1
dt = dx*C/c;
%field variables
cn=zeros(Nx,Ny,Nz); %concentration
x=linspace(0,Lx,Nx); %xdistance
y=linspace(0,Ly,Ny); %ydistance
z=linspace(0,Lz,Nz); %zdistance
[X,Y,Z]=meshgrid(x,y,z); %defining 3D grid
%Value of Diffusion
D = ones(Nx,Ny,Nz)+diff;
D([1 end],:,:) = 10^-15; %insulated problem which means that the diffusion is almost zero at the boundaries 
D(:,[1 end],:) = 10^-15;
D(:,:,[1 end]) = 10^-15;
%initial condition
cn(:,:,:) = 0.15*10^-9; %initial value of c
t=0;
for n=1:Nt
   cc = cn; 
    t=t+dt; %new time
     
%New temperature
for i=2:Nx-1
    for j=2:Ny-1
        for k=2:Nz-1
            cn(i,j,k) = cc(i,j,k)+ dt*D(i,j,k)*...
            ((cc(i+1,j,k) - 2*cc(i,j,k) + cc(i-1,j,k))/dx/dx +...
            (cc(i,j+1,k) - 2*cc(i,j,k) + cc(i,j-1,k))/dy/dy +...
            (cc(i,j,k+1) - 2*cc(i,j,k) + cc(i,j,k-1))/dz/dz);
        end
    end
end
%Insulation conditions
cn(1,:,:)=cn(2,:,:); %walls


cn(end,:,:) = cn(end-1,:,:); %walls
cn(:,1,:)=cn(:,2,:); %walls
cn(:,end,:) = cn(:,end-1,:); %walls 
cn(:,:,1)=cn(:,:,2); %floor
cn(:,:,end) = cn(:,:,end-1); %roof
%sink term Filter
cn(10,10,5) = 4*10^-11; 
%Visualization 
    if(mod(n,600) == 0) %updates image every 0.25 minutes, this has been done to speed up computation 
    slice(X,Y,Z,cn,5,5,0); colormap(flipud(hsv(256))); colorbar; caxis([3*10^-11 1.5*10^-10]);
    title(sprintf('time = %f minutes', t/60))
    pause(0.01);
    end
end

Best Answer

For the people who are interested. This line is wrong: D = ones(Nx,Ny,Nz)+diff; Instead, first define
D = ones(Nx,Ny,Nz) and next line say D(:,:,:) = diff; I think this gives an accurate representation of the diffusion. Here is the code:
%field variables
cn=zeros(Nx,Ny,Nz); %concentration matrix
D = ones(Nx,Ny,Nz); %Diffusion matrix
x=linspace(0,Lx,Nx); %xdistance
y=linspace(0,Ly,Ny); %ydistance
z=linspace(0,Lz,Nz); %zdistance
[X,Y,Z]=meshgrid(x,y,z); %defining 3D grid
%Value of Diffusion
D(:,:,:) = diff;
D([1 end],:,:) = c_border; %insulated problem which means that Diffusion is almost zero at the boundaries end does both sides
D(:,[1 end],:) = c_border;
D(:,:,[1 end]) = c_border;
Related Question