Computation of diffusion time

Question

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?

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

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