Problem: 3D Diffusion Equation with Sinkterm

9 次查看（过去 30 天）
Yorick de
Yorick de 2019-12-18
回答： Yorick de 2019-12-19
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

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

采纳的回答

Yorick de
Yorick de 2019-12-19
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;

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Eigenvalue Problems 的更多信息

标签

产品


版本

R2018a

Community Treasure Hunt

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

Start Hunting!

Translated by