# Simulating 1D diffusion

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.

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

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

• Thanks a lot for the clarification. I'd expected the mass to remain constant by imposing zero diffusive flux at both the boundaries. I thought the mass at the inlet node will diffuse to the other region of the domain (In my other post, I had used dirichlet boundary condition at the inlet node, this defines a contant flux at the inlet at all time.But what I want to try here is provide inlet flux only at t=0). I understand I was completely wrong. I've a few more doubts. Could you please have an edit made in my original post? – Natasha Feb 17 '20 at 2:45
• @Natasha I think you have some misunderstanding here. The infinite gradient is in your initial condition when you jump from zero everywhere to 3 at $x=0$. You don't have any transition region and that's the reason why an infinite concentration gradient is in your initial condition. Also, there is no mass to diffuse here at all. In a 1D domain, a point concentration is meaningless, because of that integral that is equal to mass is zero. So, you need to change your boundary or initial conditions or both. By the way, your analysis in your edited question by finite difference is wrong. – Alone Programmer Feb 17 '20 at 3:42
• Sorry for being dumb :'(. But I still don't understand what you mean by  infinite gradient is in your initial condition when you jump from zero everywhere to 3 at x=0. Could you please elaborate? I'm not able to figure out what's wrong in my edit too. I've used central difference for the second derivative at the inlet node and central difference for the left boundary (first derivative) – Natasha Feb 17 '20 at 3:51
• @Natasha Your initial condition is just equal to 3 at $x=0$ and suddenly without any transition region becomes zero everywhere except at $x=0$. This means your initial condition is not a continuous function. In fact, the gradient of this function at $x=0$ is undefined but you could think of it as an infinitely large gradient because without any change in $x$, concentration jumps from 3 to 0 and that gives you an infinite gradient and that's the problem of your formulation. – Alone Programmer Feb 17 '20 at 13:54
• @Natasha Just for sake of training, try this 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}$$ and verify your solution by seeing that for $t \rightarrow \infty$, concentration should go to $C(x,t) \rightarrow \frac{C_{L} + C_{R}}{2}$. – Alone Programmer Feb 17 '20 at 16:30