Episode 05.
Radial Diffusion in Cylindrical Solids - Unsteady State
Summary
In this episode, we describe the diffusion of a cylindrical solid in an unsteady state. First, we explain the Bessel equation to find the solution for Fick’s 2nd law in cylindrical coordinates. Then, we discuss how to calculate the concentration for various diffusion scenarios when the cylindrical solid diffuses only in the radial direction.
1. Introduction
In this episode, we discuss the diffusion of cylindrical solids in unsteady states. To describe the diffusion in the cylindrical shape, we transform Fick’s 2nd Law into cylindrical coordinates, as shown in (Eq 2.3).
\( {\small\displaystyle\begin{align} \frac{\partial C}{\partial t} \quad=\quad \frac{D}{r}\left[\frac{\partial}{\partial r}\left(r\frac{\partial C}{\partial r}\right)\,+\,\frac{\partial}{\partial\theta}\left(\frac{1}{r}\frac{\partial C}{\partial\theta}\right)\,+\,\frac{\partial}{\partial z}\left(r\frac{\partial C}{\partial z}\right)\right] \tag{Eq 2.3} \end{align} } \)
When dealing with cylindrical coordinates, we often encounter the form of a Bessel equation. So, what is the Bessel equation? Let’s briefly explain this before delving into the discussions on each diffusion scenario.
Bessel Equation
\( {\small\displaystyle\begin{align} x^2y^{\prime\prime}\,+\,xy^{\prime}\,+\,\left(x^2\,-\,n^2\right)y \quad=\quad 0 \tag*{} \end{align} } \)
This type of equation is referred to as the Bessel equation of order n. The general solution for this equation is as follows.
\( {\small\displaystyle\begin{align} y(x) \quad=\quad \alpha J_n(x)\,+\,\beta Y_n(x) \tag*{} \end{align} } \)
Applying the Frobenius method to the Bessel equation, we can obtain a particular solution, \({\small\displaystyle J_n(x)}\). Here, we refer to \({\small\displaystyle J_n(x)}\) as the Bessel function of the first kind.
\( {\small\displaystyle\begin{align} J_n(x) \quad=\quad \sum_{m=0}^{\infty}\,\frac{(-1)^m}{m!(n+m)!}\left(\frac{x}{2}\right)^{2m+n} \tag*{} \end{align} } \)
\({\small\displaystyle Y_n(x)}\) is the Bessel function of the second kind. Since it diverges as x approaches 0, we set the coefficient \({\small\displaystyle \beta}\) to 0 when dealing with various diffusion cases below.
2. Various cases of radial diffusion in cylindrical solids
Let’s assume the component within the cylindrical solid diffuses only in the radial direction. Then, we can describe various diffusion cases based on the following equation.(The diffusion coefficient is constant.)
\( {\small\displaystyle\begin{align} \frac{\partial C}{\partial t} \quad=\quad D\left(\frac{1}{r}\frac{\partial C}{\partial r}\,+\,\frac{\partial^2 C}{\partial r^2}\right) \tag*{} \end{align} } \)
(Figure 5.1)
Case 1: when the boundary condition is 0
\( {\small\displaystyle\begin{align} \frac{\partial C}{\partial t} \quad&=\quad D\left(\frac{1}{r}\frac{\partial C}{\partial r}\,+\,\frac{\partial^2 C}{\partial r^2}\right),\qquad(0\lt r \lt a) \\[15pt] C(a,t) \quad&=\quad 0,\qquad\text{for all }\,t \\[15pt] C(r,0) \quad&=\quad f(r) \tag*{} \end{align} } \)
First, let’s consider the case where the concentration at the boundary is 0. Similar to Episode 04, we can use the separation of variables method to derive the concentration for this case. Now, let’s use this method to find the solution. Here, R is a function of the radius, and T is a function of time.
\( {\small\displaystyle\begin{align} C(r,t) \quad=\quad R(r)\cdot T(t) \tag*{} \end{align} } \)
| Step 1: Separation of variables | ||||
|---|---|---|---|---|
| Suppose C | \( {\small\displaystyle\begin{align} C(r,t) \quad=\quad R(r)T(t) \end{align} } \) | \( {\small\displaystyle\begin{align} \longrightarrow \end{align} } \) | \( {\small\displaystyle\begin{align} RT^{\prime} \quad=\quad \left(\frac{R^{\prime}}{r}\,+\,R^{\prime\prime}\right)DT \end{align} } \) | |
| \( {\small\displaystyle\begin{align} \frac{1}{D}\frac{T^{\prime}}{T} \,=\, \frac{R^{\prime}}{rR}\,+\,\frac{R^{\prime\prime}}{R} \,=\, -\sigma^2 \end{align} } \) | \( {\small\displaystyle\begin{align} \longrightarrow \end{align} } \) | \( {\small\displaystyle\begin{align} \mathrm{i})&\quad T^{\prime}\,+\,D\sigma^2T \quad=\quad 0 \\[10pt] \mathrm{ii})&\quad R^{\prime\prime}\,+\,\frac{R^{\prime}}{r}\,+\,\sigma^2R \quad=\quad 0 \end{align} } \) | ||
| \( {\small\displaystyle\begin{align} \mathrm{i})&\quad T\quad=\quad pe^{-D\sigma^2t} \\[10pt] \mathrm{ii})&\quad R\quad=\quad \alpha J_0(\sigma r) \end{align} } \) | \( {\small\displaystyle\begin{align} \longrightarrow \end{align} } \) | \( {\small\displaystyle\begin{align} C\quad=\quad AJ_0(\sigma r)\exp\left(-D\sigma^2t\right) \end{align} } \) | \( {\small\displaystyle\begin{align} &\mathrm{ii})\, \text{is the} \\[8pt] &\text{Bessel equation} \\[8pt] &\text{of order zero.} \end{align} } \) | |
| \( {\small\displaystyle\begin{align} C \quad=\quad \sum_{n=1}^{\infty}\,A_nJ_0(\sigma_nr)\exp\left(-D\sigma_n^2t\right) \end{align} } \) | ||||
| Step 2: Using BC | ||||
| \( {\small\displaystyle\begin{align} \sum_{n=1}^{\infty}\,A_nJ_0(\sigma_na)\exp\left(-D\sigma_n^2t\right) \quad=\quad 0 \end{align} } \) | \( {\small\displaystyle\begin{align} \longrightarrow \end{align} } \) | \( {\small\displaystyle\begin{align} J_0(\sigma_na) \quad=\quad 0 \end{align} } \) | ||
| \( {\small\displaystyle\begin{align} J_0(\sigma_na) \quad=\quad 0 \end{align} } \) | \( {\small\displaystyle\begin{align} \longrightarrow \end{align} } \) | \( {\small\displaystyle\begin{align} \int_0^{a}\,r\left[J_0(\sigma_nr)\right]^2\,dr \quad=\quad \frac{a^2}{2}\left[J_0^{\prime}(\sigma_na)\right]^2 \end{align} } \) | ||
| Step 3: Using IC | ||||
| \( {\small\displaystyle\begin{align} C(r,0) \quad=\quad \sum_{n=1}^{\infty}A_nJ_0(\sigma_nr) \end{align} } \) | \( {\small\displaystyle\begin{align} \longrightarrow \end{align} } \) | \( {\small\displaystyle\begin{align} f(r)\quad=\quad A_1J_0(\sigma_1r)\,+\,A_2J_0(\sigma_2r)\,+\,A_3J_0(\sigma_3r)\,+\cdots \end{align} } \) | ||
| \( {\small\displaystyle\begin{align} \int_0^{a}rf(r)J_0(\sigma_nr)\,dr \quad=\quad A_n\int_0^{a}\,r\left[J_0(\sigma_nr)\right]^2\,dr \end{align} } \) | \( {\small\displaystyle\begin{align} \longrightarrow \end{align} } \) | \( {\small\displaystyle\begin{align} A_n \quad=\quad \frac{2\int_0^{a}\,rf(r)J_0(\sigma_nr)\,dr}{a^2\left[J_0^{\prime}(\sigma_nr)\right]^2} \end{align} } \) | ||
\( {\small\displaystyle\begin{align} C(r,t) \quad=\quad \frac{2}{a^2} \sum_{n=1}^{\infty}\frac{\int_{0}^{a}\,rf(r)J_0\left(\sigma_nr\right)dr}{\left[J_0^{\prime}\left(\sigma_na\right)\right]^2}J_0\left(\sigma_nr\right)\exp\left(-D\sigma^2_nt\right) \tag{Eq 5.1} \end{align} } \)
Case 2: when the boundary condition is constant
\( {\small\displaystyle\begin{align} \frac{\partial C}{\partial t} \quad&=\quad D\left(\frac{1}{r}\frac{\partial C}{\partial r}\,+\,\frac{\partial^2 C}{\partial r^2}\right),\qquad(0\lt r \lt a) \\[15pt] C(a,t) \quad&=\quad C_0,\qquad\text{for all }\,t \\[15pt] C(r,0) \quad&=\quad f(r) \tag*{} \end{align} } \)
In this case, as in the previous episode, assuming that the concentration C is composed of the sum of u and v, we can solve the PDE. (Table 5.1) summarizes the conditions for each variable.
| Variables | C | u | v |
|---|---|---|---|
| Equations | \( {\small\displaystyle\begin{align} \frac{\partial C}{\partial t} \quad=\quad D\left(\frac{1}{r}\frac{\partial C}{\partial r}\,+\,\frac{\partial^2 C}{\partial r^2}\right) \end{align} } \) | \( {\small\displaystyle\begin{align} \frac{\partial u}{\partial t} \quad=\quad D\left(\frac{1}{r}\frac{\partial u}{\partial r}\,+\,\frac{\partial^2 u}{\partial r^2}\right) \end{align} } \) | \( {\small\displaystyle\begin{align} \frac{1}{r}\frac{\partial v}{\partial r}\,+\,\frac{\partial^2 v}{\partial r^2} \quad=\quad 0 \end{align} } \) |
| Boundary condition | \( {\small\displaystyle\begin{align} C(a,t) \quad=\quad C_0 \end{align} } \) | \( {\small\displaystyle\begin{align} u(a,t) \quad=\quad 0 \end{align} } \) | \( {\small\displaystyle\begin{align} v(a,t) \quad=\quad C_0 \end{align} } \) |
| Initial condition | \( {\small\displaystyle\begin{align} C(r,0) \quad=\quad f(r) \end{align} } \) | \( {\small\displaystyle\begin{align} u(r,0) \quad=\quad f(r)\,-\,v \end{align} } \) | - |
Using (Eq 5.1), we can express the concentration C as follows.
\( {\small\displaystyle\begin{align} C(r,t) \quad&=\quad u\,+\,v \\[15pt] C(r,t) \quad&=\quad \frac{2}{a^2}\sum_{n=1}^{\infty}\frac{\int_0^{a}\,r(f(r)\,-\,v)J_0(\sigma_nr)dr}{\left[J_0^{\prime}\left(\sigma_n a\right)\right]^2}J_0(\sigma_nr)\exp\left(-D\sigma_n^2t\right) \,+\, v \tag*{} \end{align} } \)
Based on the given boundary condition, v becomes:
\( {\small\displaystyle\begin{align} v \quad&=\quad C_0 \\[15pt] \text{Thus}\qquad& \\ C(r,t) \quad&=\quad C_0\,+\,\frac{2}{a^2}\sum_{n=1}^{\infty}\frac{\int_0^{a}\left[f(r)\,-\,C_0\right]rJ_0(\sigma_nr)\,dr}{\left[J_0^{\prime}(\sigma_na)\right]^2}J_0(\sigma_nr)\exp\left(-D\sigma_n^2t\right) \tag{Eq 5.2} \end{align} } \)
Let’s take a closer look at the following equation. Due to the properties of Bessel functions, we can derive the following relation for this equation.
\( {\small\displaystyle\begin{align} \int_0^{a}\,rJ_0(\sigma_nr)\,dr \quad&=\quad \left[\frac{r}{\sigma_n}\sum_{m=0}^{\infty}\frac{(-1)^m}{m!(m\,+\,1)!}\left(\frac{\sigma_nr}{2}\right)^{2m+1}\right]_0^{a} \\[10pt] &=\quad\frac{a}{\sigma_n}J_1(\sigma_na) \\[10pt] &=\quad \frac{-a}{\sigma_n}J_0^{\prime}(\sigma_na) \tag*{} \end{align} } \)
Based on the above relation, we can rearrange (Eq 5.2) as (Eq 5.3).
\( {\small\displaystyle\begin{align} C(r,t) \quad=\quad C_0\left[1\,+\,\frac{2}{a}\sum_{n=1}^{\infty}\frac{1}{\sigma_n}\frac{J_0(\sigma_nr)}{J_0^{\prime}(\sigma_na)}\exp\left(-D\sigma^2_nt\right)\right]& \\ \tag{Eq 5.3} \\ \,+\, \frac{2}{a^2}\sum_{n=1}^{\infty}& \frac{\int_0^a\,rf(r)J_0(\sigma_nr)dr}{\left[J_0^{\prime}(\sigma_na)\right]^2}J_0(\sigma_nr)\exp(-D\sigma^2_nt) \end{align} } \)
Case 3: when the boundary condition is vary with time
\( {\small\displaystyle\begin{align} \frac{\partial C}{\partial t} \quad&=\quad D\left(\frac{1}{r}\frac{\partial C}{\partial r}\,+\,\frac{\partial^2 C}{\partial r^2}\right),\qquad(0\lt r \lt a) \\[15pt] C(a,t) \quad&=\quad \phi(t) \\[15pt] C(r,0) \quad&=\quad f(r) \tag*{} \end{align} } \)
Finally, let’s consider the diffusion scenario where the boundary conditions vary with time. Similar to Case 2, assuming C is composed of the sum of u and v, we can organize each variable as shown in (Table 5.2).
| Variables | C | u | v |
|---|---|---|---|
| Equations | \( {\small\displaystyle\begin{align} \frac{\partial C}{\partial t} \quad=\quad D\left(\frac{1}{r}\frac{\partial C}{\partial r}\,+\,\frac{\partial^2 C}{\partial r^2}\right) \end{align} } \) | \( {\small\displaystyle\begin{align} \frac{\partial u}{\partial t} \quad=\quad D\left(\frac{1}{r}\frac{\partial u}{\partial r}\,+\,\frac{\partial^2 u}{\partial r^2}\right) \end{align} } \) | \( {\small\displaystyle\begin{align} \frac{\partial v}{\partial t} \quad=\quad D\left(\frac{1}{r}\frac{\partial v}{\partial r}\,+\,\frac{\partial^2 v}{\partial r^2}\right) \end{align} } \) |
| Boundary condition | \( {\small\displaystyle\begin{align} C(a,t) \quad=\quad \phi(t) \end{align} } \) | \( {\small\displaystyle\begin{align} u(a,t) \quad=\quad 0 \end{align} } \) | \( {\small\displaystyle\begin{align} v(a,t) \quad=\quad \phi(t) \end{align} } \) |
| Initial condition | \( {\small\displaystyle\begin{align} C(r,0) \quad=\quad f(r) \end{align} } \) | \( {\small\displaystyle\begin{align} u(r,0) \quad=\quad f(r) \end{align} } \) | \( {\small\displaystyle\begin{align} v(r,0) \quad=\quad 0 \end{align} } \) |
For u, we can apply (Eq 5.1) since the conditions are the same as in Case 1. For v, we can obtain the solution using the Laplace transform method.
\( {\small\displaystyle\begin{align} \mathscr{L}\left[\frac{\partial v}{\partial t}\right] \quad=\quad \mathscr{L}\left[D\left(\frac{1}{r}\frac{\partial v}{\partial r}\,+\,\frac{\partial^2 v}{\partial r^2}\right)\right] \tag{Eq 5.4} \end{align} } \)
| Step 1: Laplace transform | ||||
|---|---|---|---|---|
| From (Eq 5.4) | \( {\small\displaystyle\begin{align} sV\,-\,v(r,0) \quad=\quad D\left(\frac{\partial^2 V}{\partial r^2}\,+\,\frac{1}{r}\frac{\partial V}{\partial r}\right) \end{align} } \) | \( {\small\displaystyle\begin{align} \longrightarrow \end{align} } \) | \( {\small\displaystyle\begin{align} \frac{\partial^2 V}{\partial r^2}\,+\,\frac{1}{r}\frac{\partial V}{\partial r}\,-\,\frac{s}{D}V \quad=\quad 0 \end{align} } \) | \( {\small\displaystyle\begin{align} \mathscr{L}[v]\,&=\,V \\[10pt] v(r,0)\,&=\,0 \end{align} } \) |
| Step 2: Solve ODE | ||||
| \( {\small\displaystyle\begin{align} V(r,s)\,=\,AJ_0\left(ir\sqrt{\frac{s}{D}}\right) \,=\, CI_0\left(r\sqrt{\frac{s}{D}}\right) \end{align} } \) | \( {\small\displaystyle\begin{align} \longrightarrow \end{align} } \) | \( {\small\displaystyle\begin{align} C\quad=\quad \frac{\Phi(s)}{I_0\left(r\sqrt{\frac{s}{D}}\right)},\quad\text{where }V(a,s)\,=\,\Phi(s) \end{align} } \) | \( {\small\displaystyle\begin{align} \mathscr{L}[\phi(t)] \,=\, \Phi(s) \end{align} } \) | |
| \( {\small\displaystyle\begin{align} \text{Thus},\qquad V \,=\, \frac{\Phi(s)I_0\left(r\sqrt{\dfrac{s}{D}}\right)}{I_0\left(a\sqrt{\dfrac{s}{D}}\right)} \,=\, \Phi(s)\cdot s\cdot V_0,\qquad \text{where}\quad V_0\,=\, \frac{I_0\left(r\sqrt{\dfrac{s}{D}}\right)}{sI_0\left(a\sqrt{\dfrac{s}{D}}\right)} \end{align} } \) | ||||
| Step 3: Inverse laplace transform | ||||
| \( {\small\displaystyle\begin{align} \mathscr{L}^{-1}[V] \quad=\quad \mathscr{L}^{-1}[\Phi(s)\cdot sV_0] \end{align} } \) | \( {\small\displaystyle\begin{align} \longrightarrow \end{align} } \) | \( {\small\displaystyle\begin{align} v\quad=\quad \mathscr{L}^{-1}\left[\Phi(s)\cdot\mathscr{L}\left[\frac{\partial v_0}{\partial t}\right]\right] \end{align} } \) | ||
| Convolution | \( {\small\displaystyle\begin{align} v \quad=\quad \phi(t)*\left[\frac{\partial v_0}{\partial t}\right] \end{align} } \) | \( {\small\displaystyle\begin{align} \longrightarrow \end{align} } \) | \( {\small\displaystyle\begin{align} v \quad=\quad \int_0^{t}\,v_0^{\prime}(r,t\,-\,\lambda)\phi(\lambda)\,d\lambda \end{align} } \) | |
| *For V0 | ||||
| Condition | \( {\small\displaystyle\begin{align} \phi(t) \quad=\quad 1 \end{align} } \) | \( {\small\displaystyle\begin{align} \longrightarrow \end{align} } \) | \( {\small\displaystyle\begin{align} V_0\,=\, \frac{I_0\left(r\sqrt{\dfrac{s}{D}}\right)}{sI_0\left(a\sqrt{\dfrac{s}{D}}\right)} \end{align} } \) | \( {\small\displaystyle\begin{align} \Phi(s)\,=\,\frac{1}{s} \end{align} } \) |
| By (Eq 5.3) | \( {\small\displaystyle\begin{align} v_0 \quad&=\quad 1\,+\,\frac{2}{a}\sum_{n=1}^{\infty}\,\frac{J_0(\sigma_nr)}{\sigma_n J_0^{\prime}(\sigma_na)}\exp\left(-D\sigma_n^2t\right) \\[10pt] \frac{\partial v_0}{\partial t} \quad&=\quad -\frac{2D}{a}\sum_{n=1}^{\infty}\frac{\sigma_nJ_0(\sigma_nr)}{J_0^{\prime}(\sigma_na)}\exp\left(-D\sigma_n^2t\right) \end{align} } \) | |||
\( {\small\displaystyle\begin{align} v \quad=\quad -\,\frac{2D}{a}\sum_{n=1}^{\infty}\,\frac{\sigma_nJ_0(\sigma_nr)}{J_0^{\prime}(\sigma_na)}\exp\left(-D\sigma_n^2t\right)\int_0^{t}\,\exp\left(D\sigma_n^2\lambda\right)\phi(\lambda)\,d\lambda \tag{Eq 5.5} \end{align} } \)
As previously mentioned, since C is the sum of u and v, C becomes:
\( {\small\displaystyle\begin{align} C \quad=\quad \frac{2}{a}\sum_{n=1}^{\infty}J_0(\sigma_nr)\exp\left(-D\sigma_n^2t\right)\left[\frac{\int_0^{a}\,rf(r)J_0(\sigma_nr)\,dr}{a\left[J_0^{\prime}(\sigma_na)\right]^2}\,-\,\frac{D\sigma_n}{J_0^{\prime}(\sigma_na)}\int_0^{t}\exp\left(D\sigma_n^2\lambda\right)\phi(\lambda)d\lambda\right] \tag{Eq 5.6} \end{align} } \)
Nomenclature
\({\small\displaystyle C}\) : concentration
\({\small\displaystyle D}\) : diffuision coefficient
\({\small\displaystyle I_n}\): modified Bessel function
\({\small\displaystyle J_n}\): Bessel function of the first kind
\({\small\displaystyle r}\) : radius
\({\small\displaystyle t}\) : time
\({\small\displaystyle u,\,v}\) : variables
\({\small\displaystyle Y_n}\): Bessel function of the second kind
\({\small\displaystyle \mathscr{L} }\): Laplace transform