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 = − D d C d x J = -D \frac{\text{d} C}{\text{d} x} J = − D d x d C
where J J J is molar flux and D D D is diffusivity. Over a small control volume and
in a short period of time δ t \delta t δ t we can perform a mass balance over the
volume:
MOLES IN − MOLES OUT = Accumulated MOLES \text{MOLES IN} - \text{MOLES OUT} = \text{Accumulated MOLES} MOLES IN − MOLES OUT = Accumulated MOLES
− D ∂ C ∂ x A δ t − [ − D ( ∂ C ∂ x + ∂ 2 C ∂ x 2 δ 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. − D ∂ x ∂ C A δ t − [ − D ( ∂ x ∂ C + ∂ x 2 ∂ 2 C δ x ) ] A δ t = V Δ C .
Replacing V = A δ x V = A\delta x V = A δ x , and after some manipulation, we get:
− D ∂ C ∂ x A δ t − [ − D ( ∂ C ∂ x + ∂ 2 C ∂ x 2 δ 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. − D ∂ x ∂ C A δ t − [ − D ( ∂ x ∂ C + ∂ x 2 ∂ 2 C δ x ) ] A δ t = A δ x Δ C .
− D ∂ C ∂ x A δ t + D ( ∂ C ∂ x + ∂ 2 C ∂ x 2 δ 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. − D ∂ x ∂ C A δ t + D ( ∂ x ∂ C + ∂ x 2 ∂ 2 C δ x ) A δ t = A δ x Δ C .
− D ∂ C ∂ x + D ∂ C ∂ x + D ∂ 2 C ∂ x 2 δ 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}. − D ∂ x ∂ C + D ∂ x ∂ C + D ∂ x 2 ∂ 2 C δ x = δ x δ t Δ C .
D ∂ 2 C ∂ x 2 = Δ C δ t . D\frac{\partial^2 C}{\partial x^2}= \frac{\Delta C}{\delta
t}. D ∂ x 2 ∂ 2 C = δ t Δ C .
Taking the limit as the right hand side goes the zero we see that the above
equation approaches
∂ C ∂ t = D ∂ 2 C ∂ x 2 \frac{\partial C}{\partial t} = D \frac{\partial^2 C}{\partial x^2} ∂ t ∂ C = D ∂ x 2 ∂ 2 C
This is the 1D diffusion equation over a small control volume V V V .
Let us consider the case of transient diffusion through a semi-infinite
geometry. The boundary conditions in this case are
at t < 0 t < 0 t < 0 , C A = C A , i , C_{A} = C_{A,i}, C A = C A , i , which is the initial molar concentration.
at t > 0 t > 0 t > 0 , C A → C A , i C_{A} \to C_{A,i} C A → C A , i at x = 0 x = 0 x = 0 .
at t > 0 t > 0 t > 0 , C A = C A , s C_{A} = C_{A,s} C A = C A , s .
Let us define a similarity variable, η \eta η , that combines the variable x x x and
t t t :
η = x 4 D t \eta = \frac{x}{\sqrt{4Dt}} η = 4 D t x
We can now re-write the diffusion equation in terms of the similarity variable.
∂ C ∂ x = ∂ C ∂ η ∂ η ∂ x \frac{\partial C}{\partial x} = \frac{\partial C}{\partial \eta}\frac{\partial \eta}{\partial x} ∂ x ∂ C = ∂ η ∂ C ∂ x ∂ η
∂ η ∂ x = 1 4 D t \frac{\partial \eta}{\partial x} = \frac{1}{\sqrt{4Dt}} ∂ x ∂ η = 4 D t 1
And so, in terms of η \eta η ,
∂ C ∂ x = 1 4 D t ∂ C ∂ η \frac{\partial C}{\partial x} = \frac{1}{\sqrt{4Dt}}\frac{\partial C}{\partial \eta} ∂ x ∂ C = 4 D t 1 ∂ η ∂ C
The second derivative term can be broken up as
∂ 2 C ∂ x 2 = ∂ ∂ x ( ∂ C ∂ x ) = ∂ ∂ x ( d C d η ⋅ d η d x ) = ∂ 2 C ∂ η 2 d η d x ⋅ d η d x = 1 4 D t ∂ 2 C ∂ η 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*} ∂ x 2 ∂ 2 C = ∂ x ∂ ( ∂ x ∂ C ) = ∂ x ∂ ( d η d C ⋅ d x d η ) = ∂ η 2 ∂ 2 C d x d η ⋅ d x d η = 4 D t 1 ∂ η 2 ∂ 2 C
Taking the derivative of η \eta η ,
∂ η ∂ t = x 4 D ⋅ − t − 3 2 2 = − x 2 t 4 D t \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*} ∂ t ∂ η = 4 D x ⋅ 2 − t − 2 3 = − 2 t 4 D t x
and so,
∂ C ∂ t = d C d η d η d t = d C d η ( − x 2 t 4 D t ) \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*} ∂ t ∂ C = d η d C d t d η = d η d C ( − 2 t 4 D t x )
The diffusion equation thus becomes
1 4 D t ∂ 2 C ∂ η 2 = 1 D d C d η ( − x 2 t 4 D t ) d 2 C d η 2 = − 2 x 4 D t d C d η \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*} 4 D t 1 ∂ η 2 ∂ 2 C d η 2 d 2 C = D 1 d η d C ( − 2 t 4 D t x ) = − 4 D t 2 x d η d C
Since x = η 4 D t x = \eta\sqrt{4Dt} x = η 4 D t , the above becomes
d 2 C d η = − 2 η d C d η \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*} d η d 2 C = − 2 η d η d C
Rewriting the boundary conditions in terms of η \eta η gives
C A = C A , i C_{A} = C_{A, i} C A = C A , i at t = 0 ∴ η → ∞ t = 0 \ \therefore \ \eta \to \infty t = 0 ∴ η → ∞
C A = C A , s C_{A} = C_{A, s} C A = C A , s at t > 0 ∴ η = 0 t > 0 \ \therefore \ \eta = 0 t > 0 ∴ η = 0
C A = C A , i C_{A} = C_{A, i} C A = C A , i at t > 0 t > 0 t > 0 and as x → ∞ ∴ η = 0 x \to \infty \ \therefore \ \eta = 0 x → ∞ ∴ η = 0
The similarity variable combines the variables x x x and t t t into a single
variable η \eta η . And so, the three boundary conditions collapse into two. In
order to now solve for C C C let U = d C d η U = \frac{\text{d}C}{\text{d}\eta} U = d η d C . The above
equation then becomes,
d U d η = − 2 η U \begin{equation*}
\begin{aligned}
\frac{\text{d} U}{\text{d} \eta} = - 2\eta U
\end{aligned}
\end{equation*} d η d U = − 2 η U
and now integrating both sides,
ln U = − η 2 + C 0 U = C 1 e − η 2 \begin{equation*}
\begin{aligned}
\ln U &= -\eta^2 + C_0 \\
U &= C_1e^{-\eta^2}
\end{aligned}
\end{equation*} ln U U = − η 2 + C 0 = C 1 e − η 2
Substituting back in for U U U we can solve for C C C .
d C d η = C 1 e − η 2 C = C 1 ∫ e − η 2 d η + C 2 \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*} d η d C C = C 1 e − η 2 = C 1 ∫ e − η 2 d η + C 2
Applying the boundary conditions:
C ( η = 0 ) = C s ⇒ C 2 = C s C(\eta = 0) = C_s \Rightarrow C_2 = C_s C ( η = 0 ) = C s ⇒ C 2 = C s
C ( η → ∞ ) = C i ⇒ C i = C 1 ∫ 0 ∞ e − η 2 d η + C s C(\eta \to \infty) = C_i \Rightarrow C_i = C_1\int_0^\infty e^{-\eta^2}
\text{d} \eta + C_s C ( η → ∞ ) = C i ⇒ C i = C 1 ∫ 0 ∞ e − η 2 d η + C s
The above integral is half the Gaussian integral which can be evaluated to
π 2 \frac{\sqrt{\pi}}{2} 2 π (see this post and the resources and the end to see why
this is the case). We can, therefore, find the value of C 1 C_1 C 1 :
C i = π 2 C 1 + C s C 1 = 2 π ( C i − C s ) \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*} C i = 2 π C 1 + C s C 1 = π 2 ( C i − C s )
And so, we get the expression for C C C :
C = 2 π ( C i − C s ) ∫ 0 η e − η 2 d η + C s C − C s C i − C s = 2 π ∫ 0 η e − η 2 d η \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*} C = π 2 ( C i − C s ) C i − C s C − C s = π 2 ∫ 0 η e − η 2 d η + C s ∫ 0 η e − η 2 d η
The integral ∫ 0 η e − η 2 d η \int_0^\eta e^{-\eta^2} \text{d} \eta ∫ 0 η e − η 2 d η 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 − η 2 d η \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*} erf ( η ) ≡ π 2 ∫ 0 η e − η 2 d η
with erf ( 1 ) = 1 \text{erf}(1) = 1 erf ( 1 ) = 1 and erf ( 0 ) = 0 \text{erf}(0) = 0 erf ( 0 ) = 0 . And so, we get the expression
for the concentration profile for semi-transient diffusion through a
semi-infinite geometry,
C − C s C i − C s = erf ( η ) \begin{equation*}
\begin{aligned}
\boxed{
\frac{C - C_s}{C_i - C_s} = \text{erf}(\eta)
}
\end{aligned}
\end{equation*} C i − C s C − C s = erf ( η )