Analytical Solution of Transport Equation

Question

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?

Answer 1

Here's one way to implement the analytic solution. However, I don't seem to be getting the correct plots if I use t0 = 0 (as you mentioned in the comment). It works if you set t0 to something greater than the final time you have (like t0 = 10 for example).

I also used function handles so that A(x,t) can be evaluated at any x and t value.

format long
R = 1;
D = 900;
v = 200;
L = 60; 
co = 4100;
ci = 4000;
x = linspace(0,60,100);
tvec = 0.0:0.01:2;
sol=[];

%Use function handles so that A(x,t) can be evaluated for any x and t
A1 = @(x,t) 0.5*erfc((x.*R-t.*v)./(2*(t.*D*R).^0.5));
A2 = @(x,t) 0.5*exp(x.*v/D).*erfc((x.*R+t.*v)./(2*(D*t.*R).^0.5));
A31 = @(x,t) 0.5*(2+v*(2*L-x)/D + (t.*v^2)/(D*R))*exp(v*L/D);
A32 = @(x,t) erfc((R*(2*L-x)+t.*v)./(2*(D*t.*R).^0.5));
A41 = @(x,t) -((t.*v^2)./(pi*D*R)).^0.5;
A42 = @(x,t) exp((v*L/D)-(R./(4*t.*D)).*(2*L -x + t.*v/R).^2);
A = @(x,t) A1(x,t) + A2(x,t) + A31(x,t).*A32(x,t) + A41(x,t).*A42(x,t);

t0 = 0; %this doesn't seem to give the right solution
%t0 = 10; %this seems to work

for t = tvec %loop over time
    if (t <= t0)
        C = ci + (co - ci)*A(x,t);
    else
        C = ci + (co - ci)*A(x,t) - co*A(x,t-t0);
    end

    sol = [sol; C];
end