Implicit solution using pdepe 1D advection diffusion reaction equation

Hello, I'm trying to solve this question wits using pdepe. However I couldn't succeed it.

I tried to discreatize to solve it in python however my output graphs does not coincide with the result of the problem. If you could help me with this problem I would be appreciate it. General equation like the given: If you need the discreatize version of the equation it is giving like that:

# Constants

L = 10 # Length of the reactor (m)

W = 2 # Width of the reactor (m)

D = 1 # Depth of the reactor (m)

E = 2 # Hydraulic characteristics (m^2/hr)

U = 1 # Velocity (m/hr)

Cin = 100 # Initial concentration (mg/L)

k = 0.2 # Decay rate (hr^-1)

A = W * D # Area (m^2)

Q = A * U # Flow (m^3/hr)

V = dx*A # Volume of the each step (m^3)

Ep = (E * A) / dx

E' is defined as the Ep in the constants part

-----------------------------------------------------------------------

Output graph should be as follows:

[Torsten]

Show us your attempt with pdepe.

[Sait Mutlu Karahan]

I tried to apply on Python not in Matlab thats why Im going to try it on Matlab now.

Answer 1

[Alan Stevens]

I think you've overcomplicated the situation! Should be more like the following (which you should check very carefully, as I did it rather quickly!!):

% dC/dt = -U*dC/dx + E*d2C/dx2 -k*C
% Central difference equations
% dC/dx = (C(x+dx,t) - C(x-dx,t))/(2dx)
% d2C/dx2 = (C(x+dx,t) - 2*C(x,t) + C(x-dx,t))/dx^2
% Explicit in time
% dC/dt = (C(x,t)-C(x,t-dt))/dt

% Data 
L = 10;         % m
E = 2;          % m^2/hr
U = 1;          % m/hr
k = 0.2;        % hr^-1
cm = 100;       % mg/L
dt = 0.002/60;  % hr
dx = 0.25;      % m
tend = 10;      % hr

nx = L/dx + 1;  % number of gateposts
nt = tend/dt;   % number of timesteps

C = zeros(nx,nt);   % Initialize
C(1,:) = cm;

for i = 2:nt  % step through time
    
   for j = 2:nx-1   % step through space
       
       dCdx = (C(j+1,i-1) - C(j-1,i-1))/(2*dx);
       d2Cdx2 = (C(j+1,i-1) - 2*C(j,i-1) + C(j-1,i-1))/dx^2;
       C(j,i) = C(j,i-1) + dt*(-U*dCdx + E*d2Cdx2 - k*C(j,i-1));  
       
   end
   
   dCdx = (C(nx,i-1) - C(nx-1,i-1))/dx;
   d2Cdx2 = (C(nx,i-1) - 2*C(nx-1,i-1) + C(nx-2,i-1))/dx^2;
   C(nx,i) =  C(nx,i-1) + dt*(-U*dCdx + E*d2Cdx2 - k*C(nx,i-1)); 
   
end

x = 0:dx:L;
t = [1 2 4 6 10]; % hr
iplot = t/dt;

plot(x, C(:,iplot)),grid
xlabel('x [m]'), ylabel('C [mg/L]')
axis([0 10 0 100])
legend('1hr','2hrs','4hrs','6hrs','10hrs')

[Torsten]

Nice, only the inflow condition at x=0 seems to be set a bit different according to the output graph.

[Alan Stevens]

@Torsten: Yes, the graphs don't really match if you look closely! As I said, I did it quickly - I'll leave it to Sait to do a careful check. (It looks like I didn't let cm at the inlet decay with time).