1D Advection-Diffusion

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')

Answer 1

[Alan Stevens]

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')

[Mercury_80]

Thank you so much! That looks great!

[Hassan Abdullah Saleem]

Hi Alan Stevens I made a modification on the code your solution was using explicit... you will see that I resolve the problem using the Imlicit method which unconditionally stable eve if you increses V, dx, dt etc

So baiscally insted solving the problem using the prevous known solution C(x,t) I iterate on solving the problem using the unknown C(x,t+dt)

% modefied after Alan Stevens, 2020
% Concentration Distribution of contaminant (C)
T = 6000*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 = 50;            % apatial grid sections
M = 400;           % 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))
% Implicit
%C(x,t+dt) = {C(x,t) + (D*dt/dx^2)*(C(x+dx,t+dt)+C(x-dx,t+dt))
%                    - (V*dt/dx)*(-C(x-dx,t+dt))}/(2*(D*dt/dx^2)+(V*dt/dx)+1)
F = D*dt/dx^2;
G = V*dt/dx;
index = 1;          % puls location
C(index,:) = 100;   % inlet concentration set to 100 for all time
for iter = 1:400
for c = 1:M   % timesteps
   
    for r = 2:N  % spatial steps
        C(r,c+1) = (C(r,c)+F*(C(r+1,c+1)+C(r-1,c+1)) ...
                         -G*(-C(r-1,c+1)))/(G+2*F+1);            
    end
     
end
C(r,c) = C(r,c+1);
end
t = 0:dt:T;
x = 0:dx:L;
plot(x,C(:,10:50:M/2),x,C(:,250:50:M)),grid
xlabel('x'),ylabel('C')
legend(['C at time ' num2str(M*dt/(2*24)) ' days'],['C at time ' num2str(M*dt/24) ' days']) 
xlim([0 100])
figure
[t,x] = meshgrid(t,x);
surf(t,x,C)
xlabel('time [hrs]'),ylabel('distance [m]'),zlabel('Concentration')

[Uraz Sitil]

Hi. Which implicit method did you use in this question? If you answer this question, I will appreciated. Thank you.

[Kayemba Luwaga]

This appears to be a finite element method! Implicit and or analytical methods usually demand intense computational resources! Cheers