I'm looking at the analytical solution of the convection-diffusion equation
$$\frac{\partial C}{\partial t} = D\frac{\partial ^2 C}{\partial x^2}-v\frac{\partial C}{\partial x}$$ with initial condition $$c(x,0) = 4000$$ and with Dirichlet boundary condition
$$C(x=0,t>0) = 4100$$
Neumann boundary condition $$\frac{\partial C}{\partial x}=0\text{ at } x=L ; t>0.$$
However, when the above analytical solution is coded I obtain negative values of C which is unrealistic.
function AnalyticalSoln2()
format long
R = 1;
D = 900;
v = 200;
L = 60;
co = 4100;
ci = 4000;
X = linspace(0,60,10);
t = 0:0.001:2;
sol=[];
for pos = 1:length(X)
x = X(pos);
A1 = 0.5*erfc((x.*R-t.*v)./(2*(t.*D*R).^0.5));
A2 = 0.5*exp(x.*v/D).*erfc((x.*R+t.*v)./(2*(D*t.*R).^0.5));
A31 = 0.5*(2+v*(2*L-x)/D + (t.*v^2)/(D*R))*exp(v*L/D);
A32 = erfc((R*(2*L-x)+t.*v)./(2*(D*t.*R).^0.5));
A41 = -((t.*v^2)./(pi*D*R)).^0.5;
A42 = exp((v*L/D)-(R./(4*t.*D)).*(2*L -x + t.*v/R).^2);
A = A1 + A2 + A31.*A32 + A41.*A42;
if t==0
C = ci + (co - ci)*A'
else
C = ci + (co - ci)*A' - co*A';
end
sol = horzcat(sol,C);
end
sol(1:100:end,:)
end
Whereas, the numerical solution obtained from the pdepe solver is non-negative.
Here's the solution obtained using pdepe
function DiffusionConvectionMATLAB
format short
global D
m = 0;
x = linspace(0,60,10);
t = 0.1:0.1:2;
D = 900;
sol = pdepe(m,@pdefun,@icfun,@bcfun,x,t)
function [g,f,s] = pdefun(x,t,c,DcDx)
v = 200;
g = 1;
f = D*DcDx;
s = -v*DcDx;
end
function c0 = icfun(x)
c0 = 4000;
end
function [pl,ql,pr,qr] = bcfun(xl,cl,xr,cr,t)
pl = cl-4100;
ql = 0;
pr = 0;
qr = 1;
end
end
Could someone suggest if there is any issue with the way I'm computing the solution form the analytical expression?
t
, but your if condition forC
is based onpos
, which represents a spatial location. $\endgroup$C
based ont
. Please have a look at the updated code. I'm still obtaining negative values. $\endgroup$if t==0
sincet
is a vector of all the time values. You need to change the for loop so that it loops over time instead of space, and at each time step you can evaluate the value of C at all spatial points through vectorization. See my answer below. However, setting $t_0 = 0$ doesn't seem to give the correct solutions. Are you sure it's zero? I get the right behavior if $t_0$ is large. $\endgroup$