Heat diffusion in a layered medium

Heat diffusion in a uniform one-dimensional medium is a basic problem that is treated in undergraduate texts on analysis. Its solutions are decaying sines, i.e. functions of the form sin(k x) e ^{-k2 t}. The negative exponentional reflects the fact that temperature distributions gradually smooth out; the sines come in handy, because they are orthogonal wrt. to the L₂ inner product. We can find a solution to a problem, by expanding the initial conditions and the heat source in a Fourier series. The sines are also called eigenfunctions of the differential operator.

If the heat conductivity varies across the medium, then we can no longer explicitly compute solutions, however, the smoothing out of temperature distributions still holds: if we separate time and spatial variables, we still have a a Sturm-Liouville problem in the spatial coordinate.

If we assume that that the conductivity is piecewise constant, and the transport between layers satisfies a radiation condition, then we define an inner product to match the differential operator. Using this inner product, we can even approximate the eigenfunctions and eigenvalues using a straightforward modal FEM analysis of the problem. By decomposing the initial condition and the heat sources into the eigenfunctions, we can come up with a solution, whose only approximative component is the FEM modal analysis. The time integration can be carried out analytically.

This was all researched by Arjan Schuurmans, who obtained a engineering degree in 1996 with this project. My traineeship was about coming up with an efficient implementation. The report of this traineeship is available here. (gzipped postscript 230k). Unfortunately, it is completely in Dutch.

I rewrote the above from its previous text, which will be quite unintellegible if you are not a mathematician. That text is here: The temperature in a layered medium satisfies a parabolic differential equation. It can be proved that the differentation operator D that corresponds to the differential equation is positive, self-adjoint and has a compact inverse (when a suitable inner product is chosen). Hence, D has a sequence of eigenvalues (leading to infinity), with orthogonal eigenfunctions. This system of eigenfunctions is complete. We can therefore solve the diffusion problem analogously to the unlayered problem (ie. using Fourier sequences, and the eigenfunctions take the place of the sines).

However, there is a practical problem: the eigenvalues and eigenfunctions can not be computed exactly, since that computation involves transcedental equations. Therefore the eigenfunctions are approximated by a finite-element procedure: the interval is subdivided into set of intervals, and using this subdivision, a (finite) basis B_n taken from the domain of D is constructed. The space spanned by B `approximates' the domain of D, and by projecting D onto span(B) we can approximate D.

The differentiation operator D is projected onto this finite-dimensional space, yielding a positive definite matrix. The eigenvectors of this matrix are the coefficients (in basis B) of the eigenfunctions that we are interested in.