# Computation of diffusion time

While simulating the diffusion of a substance in 1D, $$\frac{\partial C}{\partial t} = \nabla \cdot (D \nabla C).$$

I'd like to compute the diffusion time

In this link, the diffusion time is given by

$$t_D = \frac{l^2}{2D}$$

Say, l = 10 and the x-direction is discretized in steps of $$\Delta$$x = 1. i.e. xmesh = 0:$$\Delta$$x:10.

For the substance that is considered, $$t_D$$ = 5 seconds.

Does this mean the system that I am simulating will approach the steady-state concentration (at all nodes of xmesh) in a maximum of 5 seconds?

EDIT: Adding a check form numerical simulation. The following is the MATLAB code that simulates 1D diffusion system using pdepe solver.

function sol=so()
format short
global D nnode init_co find_index
m = 0;
xend = 5;
D = 500;
x = 0:1:xend;
find_index  = 0:1:xend;
t = 0:0.00001:0.5;
init_co = 1*ones(length(x),1);
nnode = length(x);
sol = pdepe(m,@pdefun,@icfun,@bcfun,x,t)
plot(t,sol)
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)
% Dirichlet at left (concentration at left boundary = 2, Neumann at right(dC/dx = 0)
pl = cl - 3;
ql = 0;
pr = 0;
qr = 1;
end
end


Result: Calculating $$t_D$$ as shown below, l = 5 nm and D = 500 $$nm^2/min$$, $$t_D$$ = 0.025 min. But from the plot of C vs time, all the curves reach a steady-state at 0.1 min. I would like to know if it is reasonable to compare $$t_D$$ and the time taken to reach a steady state. Could someone elaborate?

• “All curves reach a steady state at 0.1”. That’s not true. Jan 10, 2020 at 15:23

It's easy to derive that equation from Fick's law. You have this diffusion equation as:

$$\frac{\partial C}{\partial t} = D \nabla^{2} C$$

The mean square displacement weighted by the concentration profile is defined as:

$$\langle r^{2}(t) \rangle = \int_{\Omega} |\vec{r}|^{2} C(\vec{r},t) d^{3} \vec{r}$$

The time-evolution of this mean square displacement is calculated as :

$$\frac{d \langle r^{2} (t) \rangle}{d t} = \int_{\Omega} |\vec{r}|^{2} \frac{\partial C}{\partial t} d^{3} \vec{r}$$

Using your diffusion equation ($$\frac{\partial C}{\partial t} = D \nabla^{2} C$$), you have:

$$\frac{d \langle r^{2} (t) \rangle}{d t} = \int_{\Omega} D |\vec{r}|^{2} \nabla^{2} C d^{3} \vec{r}$$

Now, you know that: $$\nabla \cdot (|\vec{r}|^{2} \nabla C) = |\vec{r}|^{2} \nabla^{2} C + \nabla(|\vec{r}|^{2}) \cdot \nabla C$$, so:

$$\frac{d \langle r^{2} (t) \rangle}{d t} = \int_{\Omega} D \nabla \cdot (|\vec{r}|^{2} \nabla C) d^{3} \vec{r} - \int_{\Omega} D \nabla(|\vec{r}|^{2}) \cdot \nabla C d^{3} \vec{r}$$

or:

$$\frac{d \langle r^{2} (t) \rangle}{d t} = \int_{\partial \Omega} |\vec{r}|^{2} D\nabla C \cdot d \mathcal{A} - \int_{\Omega} D \nabla(|\vec{r}|^{2}) \cdot \nabla C d^{3} \vec{r}$$

Let's say you have this Neumann boundary condition as: $$-D\nabla C \cdot \mathbf{n} = 0$$ at $$\partial \Omega$$. So, the first would vanishes to zero. Finally:

$$\frac{d \langle r^{2} (t) \rangle}{d t} = -D \int_{\Omega} \nabla (|\vec{r}|^{2}) \cdot \nabla C d^{3} \vec{r} = -2D \int_{\Omega} \vec{r} \cdot \nabla C d^{3} \vec{r}$$

Furthermore, you have: $$\nabla \cdot (\vec{r} C) = \vec{r} \cdot \nabla C + (\nabla \cdot \vec{r}) C = \vec{r} \cdot \nabla C + 3 C$$, so:

$$\frac{d \langle r^{2} (t) \rangle}{d t} = -2D \Bigg(\int_{\Omega} \nabla \cdot (\vec{r} C) d^{3} \vec{r} - 3 \int_{\Omega} C(\vec{r},t) d^{3} \vec{r} \Bigg) = -2D \Bigg( \int_{\partial \Omega} C \vec{r} \cdot d \mathcal{A} - 3 \int_{\Omega} C d^{3} \vec{r} \Bigg)$$

The first term vanishes to zero when $$t \rightarrow \infty$$ due to the fact that mean displacement at the boundaries weighted by concentration must vanishes to zero. In fact, at $$t \rightarrow \infty$$, the net displacement of atoms or molecules that goes inside and outside the domain must be zero when the system reached an equilibrium. So:

$$\frac{d \langle r^{2} (t) \rangle}{d t} = 6 D \int_{\Omega} C(\vec{r},t) d^{3} \vec{r}$$

Note that $$\int_{\Omega} C(\vec{r},t) d^{3} \vec{r}$$, shows total mass, which is conserved in the domain and must be independent of time. If your concentration is defined as $$\frac{1}{\mathrm{m}^{3}}$$ or per volume, you can normalize the concentration based on total number of particles and safely take this integral as unity ($$\int_{\Omega} C(\vec{r},t) d^{3} \vec{r} = 1$$). Finally:

$$\frac{d \langle r^{2}(t) \rangle}{d t} = 6D$$

or:

$$\langle r^{2} (t) \rangle - \langle r^{2} (0) \rangle = 6Dt$$

By taking the initial condition as $$\langle r^{2} (0) \rangle = 0$$:

$$t = \frac{\langle r^{2} (t) \rangle }{6D}$$

This formula is derived for 3D. But for 1D, remember that $$\nabla \cdot \vec{r} = 1$$ instead of $$\nabla \cdot \vec{r} = 3$$ for 3D. So for 1D, you have:

$$\langle x^{2} (t) \rangle = 2Dt$$

or:

$$t = \frac{\langle x^{2} (t) \rangle}{2D}$$

So, yes, in fact the mean square displacement of atoms or particles, that eventually shows the diffusion, increases with time linearly. But, let's say your domain in 1D is a line with length $$L$$. So, that $$L$$ is your length scale and you want to know how long it will take for particles to diffuse in your whole domain, so your final diffusion time would be approximated as:

$$t_{f} = \frac{L^{2}}{2D}$$

Due to the fact that in a domain with length $$L$$, the maximum mean square displacement must be equal to $$\langle x^{2}(t_{f}) \rangle = L^{2}$$.

• Thanks a lot for the detailed illustration. Could you please have a look at the edit made in my original post? Dec 20, 2019 at 12:35