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)
function [g,f,s] = pdefun(x,t,c,DcDx)
g = 1;
f = D*DcDx;
s = 0;

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

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;

Result: enter image description here 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?

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

1 Answer 1


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


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


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


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

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

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.