Simulating 1D diffusion

Question

I'm trying to understand the influence of Neumann boundary condition while simulating 1D diffusion equation $$ \frac{\partial C}{\partial t} = \nabla \cdot (D \nabla C). $$

The initial value is set to 3 at the inlet node and the rest of the nodes are discretized to zero. Neumann boundary condition is used at both the left and right boundary to set diffusive flux to zero.

The following is a MATLAB implementation using pdepe solver.

function sol=check()
format short
m = 0;
delx = 0.25;
xend = 10; 
D = 500;
x = 0:delx:xend;
find_index  = x;
tspan = 0:0.00001:1;
init_co = [3 ; zeros(length(x)-1,1)];
nnode = length(x);

%% pdepe solver
sol = pdepe(m,@pdefun,@icfun,@bcfun,x,tspan);
figure(1)
plot(tspan,sol(:,end))
xlabel('time')
ylabel('c_{end}')
xlim([-0.01 0.5])
title('MATLAB - pdepe')
grid on


function [g,f,s] = pdefun(x,t,c,DcDx)
g = 1;
f = D*DcDx;
s = 0;
end

function c0 = icfun(x)
c0 = init_co(find(find_index==x));
end

function [pl,ql,pr,qr] = bcfun(xl,cl,xr,cr,t)
    pl = 0;
    ql = 1;
    pr = 0;
    qr = 1;
end
end

The following plot illustrates the change in C at the last node observed over time. Here, I'd like to understand why the value of C at the terminal node doesn't reach the value 3 (which is the initial value set at the left boundary at t=0).

EDIT: From the answer provided below I am trying to understand why the concentration boundary condition is infinity for the inconsistent inital and boundary condition that I've provided to my system.

At the inlet node:

$\frac{dC}{dt} = \frac{D}{\Delta x^2}(C_{i+1} - 2C{i} +C_{i-1})$

At the left boundary

$(C_{i+1} = C_{i-1})$, by equating diffusive flux equal to zero.

This implies, at the inlet, $\frac{dC}{dt} = \frac{D}{\Delta x^2}(2C_{i+1} - 2C{i})$

For the initial conditin that has been used, $\frac{dC}{dt} = \frac{D}{\Delta x^2}(2C_{i+1} - 2*3)$

i = 1 in the above.

I'm sorry for the stupid question. But I'd to understand how the concentration gradient is inferred to be infinity from the above.

EDIT2: From the suggestion received below, I tried the following initial condition.

$$C(x,0) = \begin{cases} C_{L} & 0 \leq x < \frac{L}{2} \\ C_{R} & \frac{L}{2} \leq x \leq L \end{cases}$$

$$C(x,0) = \begin{cases} 6 & 0 \leq x < \frac{L}{2} \\ 1 & \frac{L}{2} \leq x \leq L \end{cases}$$

is set by changing init_co = [3 ; zeros(length(x)-1,1)];

to init_co = [6*ones(20,1) ; ones(21,1)];

The following is the transient change in concentration that has been observed at the terminal node.

We can observe that steady state value of concentration is given by $C(x,t) \rightarrow \frac{C_{L} + C_{R}}{2}$ as t tends to $\infty$. The value observed in the plot exactly approaches $\frac{C_{L} + C_{R}}{2}$ when $\Delta x$ in the spatial direction is too small

But I am still trying to understand the infinite concentration gradient that has been mentioned in the answer provided below. The confusion here is even for the second initial condition that I've tried, the gradient in concentration at node positioned at L/2 and the subsequent node is quite high.

Answer 1

It seems still you don't specify your boundary conditions explicitly despite the suggestion given in one of your previous questions. As far as I understand from your MATLAB code, your boundary conditions don't make sense at all. You have two Neumann boundary conditions at the left and right sides of your 1D domain:

$$-D\frac{\partial C}{\partial x}|_{x=0} = 0$$

and

$$-D\frac{\partial C}{\partial x}|_{x=L} = 0$$

And this is a just pure diffusion equation:

$$\frac{\partial C}{\partial t} = \nabla \cdot (D \nabla C)$$

and your initial condition:

$$C(x,0) = C_{0}(x) = \begin{cases} C_{L} & x = 0 \\ 0 & \mathrm{otherwise} \end{cases}$$

Think about it physically. You blocked both sides of your computational domain to zero flux and it means mass can't go inside or outside of your domain. Also, you have an initial condition that obviously doesn't satisfy your boundary conditions and it makes the situation worse because of that infinity concentration gradient at $x = 0$ or left side of your domain which is in the obvious contrast to $-D\frac{\partial C}{\partial x}|_{x=0} = 0$ boundary condition. The other funny thing about your setup is that you put zero mass into your setup and then expect to have diffusion, which is nonsense. Why?:

$$M = \int_{0}^{L} C_{0}(x) dx = 0$$

Where $M$ is the mass at the $t = 0$ and this value should remain constant always for $t > 0$. It means there is no mass to diffuse here as well besides all other problem with your boundary conditions and initial condition.

Conclusion: You are forcing the system to go through a completely non-physical situation and that's the reason why it just show you something that is probably wrong. In my opinion, it should just give you that it can't find the solution and that's it.