MATLAB: 1D Advection-Diffusion

1dadvectioncontaminant transportdiffusiongroundwaterpde

Consider the 1-dimensional advection-diffusion equation for a chemical constituent, C, with a constant concentration (which can represent contamination) of 100 at x = 0 m andconcentration of 0 at x = 100.
Using finite difference methods, this equation can be applied to a variety of environmental problems. You should begin this assignment by writing out the governing equation, the finite difference formulation, and the appropriate boundary and initial conditions on paper. This will guide you as you write your computer code. Assume a diffusion constant, D, of 10-6m2/s and velocity, V, of 10-7m/s.
Problem 1.Write a computer code of your finite difference formulation using time steps and grid spacings that are appropriate for the problem. Plot the concentration profiles at time intervals that allow you to see evolution of the concentration profile. (You can do this either with a 2-dimensional plot with various lines or as a 3-dimensional plot i.e., time-space-concentration).
Problem 2. Evaluate the finite difference stability criteria. Run the program from (a) using and submit the code. What happens and why?
I've been working on this problem for a while now, but can't figure out how to set up the equaiton. Could someone help me figure out how to set it up? This is the code I have so far:
% Concentration Distribution of contaminant (C)
dt=86400; % time spacing per step (seconds)
dx=1;   % distance spacing per step
D=10E-6;     % diffusion constant 
V= 10E-7; % velocity 
C=zeros(101,101); % set up matrix C
        % Boundary Conditions 
C(1,:)=0;       % concentration  
C(101,:)=0;   % end of boundary
C(1,1)=10;    % concentrtation at source
for n=1:99  % time step 
    for i=1:99 % space step
        C(i+1,n+1)= V*dt/dx + D*dt/dx.^2;
        % 
    end 
end
figure(1);
subplot(2,1,1)
t=[0:100];
d=[0:100];
contourf(t,d,C)
c=colorbar;
title('Concentration Distribution')
c.Label.String = 'Concentration';
xlabel('Time (years)')
ylabel('Distance (m)')
colormap('cool')

Best Answer

Perhaps this will help
% Concentration Distribution of contaminant (C)
T = 4000*24;        % [hrs] total time 
D = (1E-6)*3600;   % [m^2/hr] diffusion constant 
V = (1E-7)*3600;   % [m/hr] velocity 
L = 100;           % [m] Total length
N = 20;      % apatial grid sections
M = 40;      % temporal grid sections
dx = L/N;     % spatial spacing
dt = T/M;     % time spacing
C=zeros(N+1,M+1); % allocate space for concentrations
% (C(x,t+dt) - C(x,t))/dt = -V*(C(x,t)-C(x-dx,t))/dx 
%                           +D*(C(x+dx,t)-2*C(x,t)+C(x-dx,t))/dx^2
%
% C(x,t+dt) = C(x,t) + (D*dt/dx^2)*(C(x+dx,t)-2C(x,t)+C(x-dx,t))
%                    - (V*dt/dx)*(C(x,t)-C(x-dx,t))
F = D*dt/dx^2;
G = V*dt/dx;
C(1,:) = 100;   % inlet concentration set to 100 for all time
for c = 1:M   % timesteps
   
    for r = 2:N  % spatial steps
        C(r,c+1) = C(r,c)+F*(C(r+1,c)-2*C(r,c)+C(r-1,c)) ...
                         -G*(C(r,c)-C(r-1,c));            
    end
     
end
t = 0:dt:T;
x = 0:dx:L;
plot(x,C(:,M/2),x,C(:,M)),grid
xlabel('x'),ylabel('C')
legend(['C at time ' num2str(M*dt/(2*24)) ' days'],['C at time ' num2str(M*dt/24) ' days']) 
figure
[t,x] = meshgrid(t,x);
surf(t,x,C)
xlabel('time [hrs]'),ylabel('distance [m]'),zlabel('Concentration')
Related Question