Clear Filters
Clear Filters

Implicit solution using pdepe 1D advection diffusion reaction equation

22 views (last 30 days)
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:

Answers (1)

Alan Stevens
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
Torsten on 13 Jun 2023
Nice, only the inflow condition at x=0 seems to be set a bit different according to the output graph.
Alan Stevens
Alan Stevens on 13 Jun 2023
@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).

Sign in to comment.

Categories

Find more on MATLAB in Help Center and File Exchange

Products


Release

R2022a

Community Treasure Hunt

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

Start Hunting!