skip to content
Aymen Hafeez

Given some medium, e.g. water, containing some material, say salt, the material will DIFFUSE from the concentrated area to the dilute area. Diffusion can be generalised by Fick’s law:

J=DdCdxJ = -D \frac{\text{d} C}{\text{d} x}

where JJ is molar flux and DD is diffusivity. Over a small control volume and in a short period of time δt\delta t we can perform a mass balance over the volume:

MOLES INMOLES OUT=Accumulated MOLES \text{MOLES IN} - \text{MOLES OUT} = \text{Accumulated MOLES} Control volume diagram DCxAδt[D(Cx+2Cx2δx)]Aδt=VΔC. -D\frac{\partial C}{\partial x} A \delta t - \left[ -D \left( \frac{\partial C}{\partial x} + \frac{\partial^2 C}{\partial x^2}\delta x \right) \right] A\delta t = V\Delta C.

Replacing V=AδxV = A\delta x, and after some manipulation, we get:

DCxAδt[D(Cx+2Cx2δx)]Aδt=AδxΔC. -D\frac{\partial C}{\partial x} A \delta t - \left[ -D \left( \frac{\partial C}{\partial x} + \frac{\partial^2 C}{\partial x^2}\delta x \right) \right] A\delta t = A\delta x \Delta C. DCxAδt+D(Cx+2Cx2δx)Aδt=AδxΔC. -D\frac{\partial C}{\partial x} A \delta t + D\left( \frac{\partial C}{\partial x} + \frac{\partial^2 C}{\partial x^2}\delta x \right) A\delta t = A\delta x \Delta C. DCx+DCx+D2Cx2δx=δxΔCδt. -D\frac{\partial C}{\partial x} + D\frac{\partial C}{\partial x} + D\frac{\partial^2 C}{\partial x^2}\delta x = \delta x \frac{\Delta C}{\delta t}. D2Cx2=ΔCδt. D\frac{\partial^2 C}{\partial x^2}= \frac{\Delta C}{\delta t}.

Taking the limit as the right hand side goes the zero we see that the above equation approaches

Ct=D2Cx2 \frac{\partial C}{\partial t} = D \frac{\partial^2 C}{\partial x^2}

This is the 1D diffusion equation over a small control volume VV.

Let us consider the case of transient diffusion through a semi-infinite geometry. The boundary conditions in this case are

  • at t<0t < 0, CA=CA,i,C_{A} = C_{A,i}, which is the initial molar concentration.
  • at t>0t > 0, CACA,iC_{A} \to C_{A,i} at x=0x = 0.
  • at t>0t > 0, CA=CA,sC_{A} = C_{A,s}.

Let us define a similarity variable, η\eta, that combines the variable xx and tt:

η=x4Dt \eta = \frac{x}{\sqrt{4Dt}}

We can now re-write the diffusion equation in terms of the similarity variable.

Cx=Cηηx \frac{\partial C}{\partial x} = \frac{\partial C}{\partial \eta}\frac{\partial \eta}{\partial x} ηx=14Dt \frac{\partial \eta}{\partial x} = \frac{1}{\sqrt{4Dt}}

And so, in terms of η\eta,

Cx=14DtCη \frac{\partial C}{\partial x} = \frac{1}{\sqrt{4Dt}}\frac{\partial C}{\partial \eta}

The second derivative term can be broken up as

2Cx2=x(Cx)=x(dCdηdηdx)=2Cη2dηdxdηdx=14Dt2Cη2\begin{equation*} \begin{aligned} \frac{\partial^2 C}{\partial x^2} &= \frac{\partial}{\partial x}\left(\frac{\partial C}{\partial x}\right) \\ &= \frac{\partial}{\partial x}\left(\frac{\text{d} C}{\text{d} \eta}\cdot\frac{\text{d} \eta}{\text{d} x}\right) \\ &= \frac{\partial^2 C}{\partial \eta^2} \frac{\text{d}\eta}{\text{d} x}\cdot\frac{\text{d}\eta}{\text{d} x} \\ &= \frac{1}{4Dt} \frac{\partial^2 C}{\partial \eta^2} \end{aligned} \end{equation*}

Taking the derivative of η\eta,

ηt=x4Dt322=x2t4Dt\begin{equation*} \begin{aligned} \frac{\partial \eta}{\partial t} &= \frac{x}{4D}\cdot\frac{-t^{-\frac{3}{2}}}{2} \\ &= -\frac{x}{2t\sqrt{4Dt}} \end{aligned} \end{equation*}

and so,

Ct=dCdηdηdt=dCdη(x2t4Dt)\begin{equation*} \begin{aligned} \frac{\partial C}{\partial t} &= \frac{\text{d} C}{\text{d} \eta} \frac{\text{d} \eta}{\text{d} t} \\ &= \frac{\text{d} C}{\text{d} \eta}\left(-\frac{x}{2t\sqrt{4Dt}}\right) \end{aligned} \end{equation*}

The diffusion equation thus becomes

14Dt2Cη2=1DdCdη(x2t4Dt)d2Cdη2=2x4DtdCdη\begin{equation*} \begin{aligned} \frac{1}{4Dt}\frac{\partial^2 C}{\partial \eta^2} &= \frac{1}{D}\frac{\text{d}C}{\text{d}\eta}\left(-\frac{x}{2t\sqrt{4Dt}}\right) \\ \frac{\text{d}^2 C}{\text{d} \eta^2} &= -\frac{2x}{\sqrt{4Dt}}\frac{\text{d}C}{\text{d}\eta} \end{aligned} \end{equation*}

Since x=η4Dtx = \eta\sqrt{4Dt}, the above becomes

d2Cdη=2ηdCdη\begin{equation*} \begin{aligned} \frac{\text{d}^2 C}{\text{d} \eta} = -2\eta \frac{\text{d} C}{\text{d} \eta} \end{aligned} \end{equation*}

Rewriting the boundary conditions in terms of η\eta gives

  • CA=CA,iC_{A} = C_{A, i} at t=0  ηt = 0 \ \therefore \ \eta \to \infty
  • CA=CA,sC_{A} = C_{A, s} at t>0  η=0t > 0 \ \therefore \ \eta = 0
  • CA=CA,iC_{A} = C_{A, i} at t>0t > 0 and as x  η=0x \to \infty \ \therefore \ \eta = 0

The similarity variable combines the variables xx and tt into a single variable η\eta. And so, the three boundary conditions collapse into two. In order to now solve for CC let U=dCdηU = \frac{\text{d}C}{\text{d}\eta}. The above equation then becomes,

dUdη=2ηU\begin{equation*} \begin{aligned} \frac{\text{d} U}{\text{d} \eta} = - 2\eta U \end{aligned} \end{equation*}

and now integrating both sides,

lnU=η2+C0U=C1eη2\begin{equation*} \begin{aligned} \ln U &= -\eta^2 + C_0 \\ U &= C_1e^{-\eta^2} \end{aligned} \end{equation*}

Substituting back in for UU we can solve for CC.

dCdη=C1eη2C=C1eη2dη+C2\begin{equation*} \begin{aligned} \frac{\text{d} C}{\text{d} \eta} &= C_1e^{-\eta^2} \\ C &= C_1\int e^{-\eta^2} \text{d}\eta + C_2 \end{aligned} \end{equation*}

Applying the boundary conditions:

  • C(η=0)=CsC2=CsC(\eta = 0) = C_s \Rightarrow C_2 = C_s
  • C(η)=CiCi=C10eη2dη+CsC(\eta \to \infty) = C_i \Rightarrow C_i = C_1\int_0^\infty e^{-\eta^2} \text{d} \eta + C_s

The above integral is half the Gaussian integral which can be evaluated to π2\frac{\sqrt{\pi}}{2} (see this post and the resources and the end to see why this is the case). We can, therefore, find the value of C1C_1:

Ci=π2C1+CsC1=2π(CiCs)\begin{equation*} \begin{aligned} C_i = \frac{\sqrt{\pi}}{2}C_1 + C_s \\ C_1 = \frac{2}{\sqrt{\pi}}(C_i - C_s) \end{aligned} \end{equation*}

And so, we get the expression for CC:

C=2π(CiCs)0ηeη2dη+CsCCsCiCs=2π0ηeη2dη\begin{equation*} \begin{aligned} C = \frac{2}{\sqrt{\pi}}(C_i - C_s)&\int_0^\eta e^{-\eta^2} \text{d} \eta + C_s \\ \frac{C - C_s}{C_i - C_s} = \frac{2}{\sqrt{\pi}}&\int_0^\eta e^{-\eta^2}\text{d}\eta \end{aligned} \end{equation*}

The integral 0ηeη2dη\int_0^\eta e^{-\eta^2} \text{d} \eta cannot be solved analytically (except for η\eta \to \infty), however, it can be solved numerically. The error function can be defined as

erf(η)2π0ηeη2dη\begin{equation*} \begin{aligned} \text{erf}(\eta) \equiv \frac{2}{\sqrt{\pi}} \int_0^\eta e^{-\eta^2} \text{d} \eta \end{aligned} \end{equation*}

with erf(1)=1\text{erf}(1) = 1 and erf(0)=0\text{erf}(0) = 0. And so, we get the expression for the concentration profile for semi-transient diffusion through a semi-infinite geometry,

CCsCiCs=erf(η)\begin{equation*} \begin{aligned} \boxed{ \frac{C - C_s}{C_i - C_s} = \text{erf}(\eta) } \end{aligned} \end{equation*}