# Implicit solution using pdepe 1D advection diffusion reaction equation

22 views (last 30 days)

Show older comments

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:

##### 2 Comments

### Answers (1)

Alan Stevens
on 13 Jun 2023

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

##### 2 Comments

Torsten
on 13 Jun 2023

Alan Stevens
on 13 Jun 2023

### See Also

### Community Treasure Hunt

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

Start Hunting!