Implicit solution using pdepe 1D advection diffusion reaction equation
4 views (last 30 days)
Show older comments
Hello, I'm trying to solve this question wits using pdepe. However I couldn't succeed it. ![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1409049/image.jpeg)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1409049/image.jpeg)
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:
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1409089/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1409059/image.jpeg)
# 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:![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1409064/image.jpeg)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1409064/image.jpeg)
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
Nice, only the inflow condition at x=0 seems to be set a bit different according to the output graph.
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).
See Also
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!