~  Heat conduction partial differential equation
 Explicit method from Carnahan, Luther & Wilkes,
 "Applied Numerical Methods",  Wiley, 1969, p. 432
 Programmed by Victor Bloomfield (victor@biosci.cbs.umn.edu)
 11/11/96
~

--Heat conduction in a rod of length L with ends held at
--T_1 and T_2 and initial internal temperature T_in
--Reduced time tau=a*t/L^2, a=thermal diffusivity k/(rho*cp)
--Reduced position = x/L
--To calculate T(x,t) from the heat conduction equation in
--reduced coordinates: dT/dt = d^2 T/ dx^2.

T_in=70; T_1=70; T_2=0

nx=13 -- number of spatial grid points, including ends
dx=1/(nx-1) -- length step in reduced length units x/L
tmax=.5 -- in reduced time units
nt=201 -- number of time steps
dt=tmax/(nt-1) -- time step in reduced units
r=dt/(dx)^2;
r:0.360;  -- lambda in Carnahan et al (should be 2.5 for convergence)

--Set up T matrix, including initial condition T=T_in everywhere
zeros2[i,j] = T_in dim[nx,nt]
T:=zeros2:

--Boundary conditions
T1:=T_1[:nt]:; T2:=T_2[:nt]:
k:=1,(T[1,k]:=T1[k],T[nx,k]:=T2[k],k:=k+1) while k<=nt:

--Step through grid
 k:=2,
 ((i:=2,
  (T[i,k]:=r*T[i-1,k-1]+(1-2*r)*T[i,k-1]+r*T[i+1,k-1],
   i:=i+1) while i<nx),
  k:=k+1) while k<=nt:

a=.02; L=1 -- Brick wall L ft thick
maxtime=tmax*L^2/a; maxtime:25.000--hr

--Plot temperature as f(x) at various times
Xlabel="x"; Ylabel="T"
Xmin=0; Xmax=L; Xsteps=nx; Ymin=0
tau={.1,.2,.5}*tmax
tau*L^2/a:{2.500,5.000,12.500}-- time in hrs
plot T[X/L*(nx-1)+1,nt*tau[1]]
plot T[X/L*(nx-1)+1,nt*tau[2]]
plot T[X/L*(nx-1)+1,nt*tau[3]]

--Plot temperature as f(t) at various positions
newaxis
Xmin:=0:; Xmax:=maxtime:; Xsteps:=nt:
Xlabel:="time":; Ylabel:="T":
hrs=X*L^2/a
plot T[nx/4,(X/maxtime*(nt-1)+1)]
plot T[nx/2,(X/maxtime*(nt-1)+1)]
plot T[3*nx/4,(X/maxtime*(nt-1)+1)]

~ Thermal diffusivities a of some common materials in ft^2/hr
To convert to SI units (m^2/s) multiply by 2.58e-5.
brick:   0.02 (glass and limestone are similar)
copper   4.42
wood     0.004
water    0.00554 (at 20 C)
air      0.8587 (at 27 C)
~