Unsteady-State Diffusion and Reaction in a Semi-infinite Medium 
Author's Data: Housam BINOUS
Department of Chemical Engineering
National Institute of Applied Sciences and Technology
Tunis, TUNISIA
Email: binoushousam@yahoo.com

Reference:C. J. Geankoplis, Transport Processes and Unit Operations, 3rd edition, Prentice Hall, 1993.

Governing equation and boundary conditions are:

Analytical solution was derived by P. V. Danckwerts, Trans. Faraday Soc., 46, 300 (1950) :

Diffusion coefficient is (Subscript[D, AB]). Rate constant is k. c0 is surface concentration.

(Subscript[D, AB]) = 1.5 10^-9;

k=35;

c0=3 10^-2;

Case 1: Penetration depth is limited by reaction: ((Subscript[δ, d])/(Subscript[δ, r]))>>1

t=10^2;

(Subscript[δ, d]) = Sqrt[(Subscript[D, AB]) t]

(Subscript[δ, r]) = Sqrt[(Subscript[D, AB])/k]

((Subscript[δ, d])/(Subscript[δ, r]))

c=.;t=.;

Solving the governing equation with NDSolve using the following boundary condition c[10 (Subscript[δ, r])]=0 with (Subscript[δ, r])~(10^(-5)).

sol = 
      NDSolve[{D[c[t, x], {t, 1}] == 
            (Subscript[D, AB]) D[c[t, x], {x, 2}] - k c[t, x], 
          c[0, x] == c0 Exp[-10^20 x], c[t, 0] == c0, 
          c[t, 1 10^-4] == 0}, 
        c[t, x], {x, 0, 1 10^-4}, {t, 0, 100}]

plt1=Plot[Chop[c[t,x]/.sol/.t->10^2],{x,0,1 10^-4},PlotRange->{0,0.03},PlotStyle->RGBColor[1,0,0]]

Danckwerts' solution:

c = 
      c0 (0.5 Exp[-x (k/(Subscript[D, AB]))^0.5] Erfc[
                x/(2 (t (Subscript[D, AB]))^0.5) - (t k)^0.5] + 
            0.5 Exp[x (k/(Subscript[D, AB]))^0.5] Erfc[
                x/(2 (t (Subscript[D, AB]))^0.5) + (t k)^0.5])

plt2=Plot[Chop[c/.t->10^2],{x,0,10^-4},PlotRange->{0,0.03},PlotStyle->RGBColor[0,0,1]]

Both methods give the same concentration profile

Show[plt1,plt2]

Case 2: Penetration depth is limited by diffusion: ((Subscript[δ, d])/(Subscript[δ, r]))<<1

t=10^-7;

(Subscript[δ, d]) = Sqrt[(Subscript[D, AB]) t]

(Subscript[δ, r]) = Sqrt[(Subscript[D, AB])/k]

((Subscript[δ, d])/(Subscript[δ, r]))

c=.;t=.;

Solving the governing equation with NDSolve using the following boundary condition c[10 (Subscript[δ, d])]=0 with (Subscript[δ, d])~(10^(-8)).

sol = 
      NDSolve[{D[c[t, x], {t, 1}] == 
            (Subscript[D, AB]) D[c[t, x], {x, 2}] - k c[t, x], 
          c[0, x] == c0 Exp[-10^20 x], c[t, 0] == c0, 
          c[t, 1 10^-7] == 0}, 
        c[t, x], {x, 0, 1 10^-7}, {t, 0, 10^-5}]

plt1=Plot[Chop[c[t,x]/.sol/.t->10^-7],{x,0,1 10^-7},PlotRange->{0,0.03},PlotStyle->RGBColor[1,0,0]]

Danckwerts' solution:

c = 
      c0 (0.5 Exp[-x (k/(Subscript[D, AB]))^0.5] Erfc[
                x/(2 (t (Subscript[D, AB]))^0.5) - (t k)^0.5] + 
            0.5 Exp[x (k/(Subscript[D, AB]))^0.5] Erfc[
                x/(2 (t (Subscript[D, AB]))^0.5) + (t k)^0.5])

plt2=Plot[Chop[c/.t->10^-7],{x,0,1 10^-7},PlotRange->{0,0.03},PlotStyle->RGBColor[0,0,1]]

Both methods give the same concentration profile

Show[plt1,plt2]