Transient diffusion of mass within an infinitely long cylinder: The finite Hankel transform  
Author's Data: Housam BINOUS
Department of Chemical Engineering
National Institute of Applied Sciences and Technology
Tunis, TUNISIA
Email: binoushousam@yahoo.com

A Sheel balance gives the following PDE

The governing equations become when we introduce the following dimensionless quantities:

		u=(C/C0), s=(r/a)and τ=((D t)/(a^2))

We introduce the variable v=u-1 because we need homogeneous bounday conditions. The governing equations become:

Solution of the Eigenvalue problem:

eq=D[s D[ϕ[s],{s,1}],{s,1}]==-λ s ϕ[s]

sol1=DSolve[eq,ϕ[s],s]

ϕ[s] = 
      First[ϕ[s] /. sol1 /. C[1] -> C1 /. 
          C[2] -> C2]

C2 must be set equal to zero because ϕ is finite

ϕ[s] = ϕ[s] /. C2 -> 0

Computing the first nine eigenvalues (Subscript[λ, n])

For[i = 1, i < 10, 
      i += 1, {Print[(Subscript[λ, i]), "=", (Subscript[λ, i]) = FindRoot[
                  BesselJ[0, x] == 0, {x, i Pi}][[1, 
                  2]]^2]}]

Computing C1 in order to have an orthonormal set of eigenfunctions

Solve[{Integrate[s ϕ[s] ϕ[s], {s, 0, 1}] == 1, 
        BesselJ[0, (Sqrt[λ])] == 0}, C1]

ϕ[s] /. λ -> (Subscript[λ, n]) /. 
      C1 -> ((Sqrt[2])/BesselJ)[1, (Sqrt[(Subscript[λ, n])])]

(Subscript[ϕ, i]) for i = 
        1 to 9 are the first nine eigenfunctions

For[i = 0, i < 10, 
      i += 1, {n = 
          i, (Subscript[ϕ, i])[s] = ϕ[s] /. λ -> (Subscript[λ, n]) /. 
            C1 -> ((Sqrt[2])/BesselJ)[1, (Sqrt[(Subscript[λ, n])])]}]

Taking the finite Hankel transform of the governing equation leads to the following ordinary differential equation:

sol2 = 
      DSolve[D[v[τ], {τ, 1}] ==  -λ v[τ], 
        v[τ], τ]

First[v[τ] /. sol2 /. C[1] -> V1]

n=.

h[s]=Ci/C0-1

V1 = 
      Integrate[s h[s] ϕ[s] /. λ -> (Subscript[λ, n]) /. 
          C1 -> ((Sqrt[2])/BesselJ)[1, (Sqrt[(Subscript[λ, n])])] , {s, 0, 1}]

(Subscript[v, i]) for i = 1 to 9

For[i = 0, i < 10, 
      i += 1, {n = i, 
        (Subscript[v, i])[τ] = v[τ] /. sol2 /. 
              C[1] -> 
                V1 /. λ -> (Subscript[λ, n])}]

n=.

Taking the inverse transform gives the complete solution v(s,τ), u(s,τ) and c(s,τ) :

v[s_, τ_] := (∑_(i = 1)^9) (Subscript[v, i])[τ] (Subscript[ϕ, i])[
            s]

u[s_,τ_]:=1+v[s,τ]

cdim[r_,t_]:=First[C0 u[s,τ]/.τ->D t/a^2/.s->r/a]

Average concentration within the cylinder

CAverage=Integrate[cdim[r,t] 2 Pi r,{r,0,a}]/Integrate[2 Pi r,{r,0,a}]//Simplify

Plotting the dimensionless function of time X(t) defined as:

DC = Expand[CAverage - C0 // FullSimplify, C0 - Ci]

X[t_] := 
      Log[4 (Ci - C0)/((Subscript[λ, 1]) DC)] /. D/a^2 -> 0.1 // 
        Simplify

The plot of X(t) is the same as the one in figure 7.7 of Mathematical Methods in Chemical Engineering by Varma and Morbidelli, page 515, Oxford University Press, 1997.

Plot[X[t],{t,0,3},PlotStyle->RGBColor[1,0,0],PlotRange->{-1,2}]