Unsteady-State Conduction 
and the Schmidt Numerical Method 
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.

We use the Schmidt numerical method with M=2. 

Thermal diffusivity is α=2 (10^(-5)) (m^2)/s. 

The slab is 1m thick. We choose Δx=0.05m. Thus, 1m/0.05m=20 slices will be used. 

We want to compute the temperature profile at (Subscript[t, final]) = 9600  s. Since M=2=\(TraditionalForm\`Δx\^2\/\(α\ Δt\)\), we get Δt=62.5s. Hence, (Subscript[t, final])/Δt=9600/62.5=96 time increments will be used. 

Since M=2, temperature at the next time step is computed using the following relation: T[t,x]==1/2 (T[t-1,x-1]+T[t-1,x+1])

Temperature at x=0 is fixed equal to 0°C.

The slab is insulated at x=1m.

Initial temperature is 100°C expect at x=0 where we take it equal to (100+0)/2=50°C.

Eq[t_,x_]:=T[t,x]==1/2 (T[t-1,x-1]+T[t-1,x+1])

Eq[t_,0]:=T[t,0]==0.

Eq[t_,20]:=T[t,20]==T[t-1,19]

T[0,x_]:=If[x>0,100.,50.]

t=1;

While[t<97,x=0;While[x<21,sol=Solve[Eq[t,x],T[t,x]];T[t,x]=T[t,x]/.sol[[1]];x+=1];t+=1]

Results are the same as those presented in Table 5.4-1 page 356 C. J. Geankoplis, Transport Processes and Unit Operations, 3rd edition, Prentice Hall, 1993.

Print[T[96,0]," ",T[96,4]," ",T[96,8]," ",T[96,12]," ",T[96,16]," ",T[96,20]]

Ploting the temperature profile at different time increments.

tbl[i_]:=Table[{x,T[10*i+6,x]},{x,0,20}];

plt[i_]:=ListPlot[tbl[i],PlotJoined->True,PlotStyle->Hue[i/10],DisplayFunction->Identity]

Show[Table[plt[i],{i,0,9}],DisplayFunction->$DisplayFunction]

Solving the same problem using NDSolve

t=.;x=.;

solution=u[x, t]/.First[NDSolve[{\(∂\_t u[x, t] == 
                      2\ 10^\(-5\)\ \ ∂\_{x, 2}u[x, t]\),u[x, 0] == 100,u[0, t] == 100 Exp[-10000 t],Derivative[1, 0][u][1, t]==0},u[x, t],{x, 0, 1},{t, 0, 6000}]];

Off[General::spell]

plot=Plot[solution/.t->6000/.x->N 0.05,{N,0,20},PlotStyle->RGBColor[0,1,0]]

Comparing the result obtained using the Schmidt numerical method and NDSolve.

Show[plot,plt[9]]