2
$\begingroup$

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.$$

enter image description here

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?

$\endgroup$
  • $\begingroup$ It seems like you are not implementing the analytic solution correctly. What is the value of $t_0$ in your problem? The analytic solution has two branches based on the value of t, but your if condition for C is based on pos, which represents a spatial location. $\endgroup$ – Savithru Jan 1 at 17:39
  • $\begingroup$ @Savithru Sorry, that was my mistake. The value of $t_0$ in my problem is 0. I'm now modifying the condition for C based on t. Please have a look at the updated code. I'm still obtaining negative values. $\endgroup$ – Natasha Jan 2 at 3:38
  • $\begingroup$ You cannot do if t==0 since t 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$ – Savithru Jan 2 at 4:55
1
$\begingroup$

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
$\endgroup$
  • $\begingroup$ Thanks for correcting another silly mistake,t0=10 works. I'm confused on how to compare the solutions obtained from MATLAB against the analytical solution obtained using t0=10. Also, could you please clarify how the xmesh option in MATLAB works. Say, x = linspace(0,60,100); /x = linspace(0,60,10); .Does MATLAB discretize the space according to how x is defined by the user? Or the user-defined x is just the point in space at which the solution is output? $\endgroup$ – Natasha Jan 2 at 5:34
  • $\begingroup$ t=10 doesn't work too. Since , tvec = 0.0:0.01:2; is till 2 , t0 cannot be greater than 2. I also tried changing tvec = 0.0:0.01:20; and setting t0=10 .This gives negatives values again. $\endgroup$ – Natasha Jan 2 at 5:45
  • 1
    $\begingroup$ I found the source from where you got the analytical solution. If you read the boundary condition definition in that link, it states that $t_0$ is the time until which the left Dirichlet BC enforces a value of $C_0$. When $t > t_0$, the Dirichlet condition becomes $c(0,t) = 0$. But in your problem, the Dirichlet condition is $c(0,t) = C_0$ for all time. This means that $t_0$ is effectively infinity in your case. $\endgroup$ – Savithru Jan 2 at 19:34
  • 1
    $\begingroup$ So you can either set t0 = Inf in the code, or better, remove the if condition inside the time loop and always use C = ci + (co - ci)*A(x,t). $\endgroup$ – Savithru Jan 2 at 19:36
  • $\begingroup$ A silly doubt, The numerical solution is in accordance with the analytical solution for the convection-diffusion equation. But, when I try to compare the results for $$\frac{\partial C}{\partial t} = D\frac{\partial ^2 C}{\partial x^2}$$ alone by setting v=0 in the analytical solutions, the results are not comparable. Wouldn't the analytical solution of the diffusion equation be obtained when v= 0 in the above analytical solution? $\endgroup$ – Natasha Jan 14 at 6:10

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.