Computing Heat Flow

This practice was originally called difference equations emphasizing the connection to the “underlying” differential equations. More recently this practice is called finite element analysis and sometimes described as an approximation to “integral equations” that themselves are nearly equivalent to differential equations.

The first computing plan usually tried in such attempts is called explicit. In an explicit calculation field values are somehow known for all of the zones for some particular time and the corresponding values for the next time are computed from those of the current time in the current and immediately adjacent zones. This means that cause and effect can, within the simulation, travel thru the zones no faster than mesh velocity—one zone per time step.

I will talk here first about the one dimensional heat diffusion equation for it illustrates an interesting problem and several interesting general techniques for doing such simulations. Higher dimensional meshes are discussed too.

Imagine a rod where we are interested in the temperatures of equally spaced points along the rod and how they change over time. Suppose that heat flowing along the rod is the significant factor in knowing what the temperatures will be. Here is an interactive presentation of the heat equation, and another. Hear this one.

Let ϕ^k_j denote the temperature in zone j at time t^k. We can then compute directly from the equations:
[ϕ^k+1_j − ϕ^k_j]/Δt = [e_j(ϕ^k_j−1 − ϕ^k_j) + f_j(ϕ^k_j+1 − ϕ^k_j)]/Δx²
(This comes directly from ∂ϕ/∂t = e(x)∂²ϕ/∂x²)
The expressions in the parenthesis are temperature difference between adjacent zones and those divided by Δx are the temperature gradients which determine the rate at which heat is flowing to point k. The values e_j and f_j are conductivity of the bar and may in some situations depend on ϕ. These equations determine how fast each point is heating up, which is just the term on the left of the equation. I have assumed heat capacity to be one to simplify equations. (Alternatively varying heat capacity can be coded in the e’s and f’s.)

These equations can easily be put into the form
ϕ^k+1_j = a_jϕ^k_j−1 + b_jϕ^k_j + c_jϕ^k_j+1

The a’s, b’s and c’s depend on Δt, Δx, the e’s and f’s. Choosing a large Δt would seem to get our answer sooner but our equations, as they stand, say that the heat flows during all of a time interval at the rate determined by the gradients at the beginning of the interval whereas in fact the temperature gradients may significantly decrease during the interval. For Δt too large the temperature changes will overshoot and oscillate wildly and soon lose all fidelity to the original problem and differential equation. The heat thus tries to flow too fast and this is a truncation error which we must deal with.

The following implicit equations are much more accurate.
[ϕ^k+1_j − ϕ^k_j]/Δt = [e_j(ϕ^k+½_j−1 − ϕ^k+½_j) + f_j(ϕ^k+½_j+1 − ϕ^k+½_j)]/Δx²
where ϕ^k+½_j = (ϕ^k_j + ϕ^k+1_j)/2.
This says that the flow rate is determined by the temperature differences averaged over the time interval. Actually it is not quite that good; it uses the average of the beginning and ending gradients. These latter equations are called second order for the error is proportional to Δt². The explicit equations are first order.

But it is not now immediately clear how to compute the k+1 values from the known k values. These equations can be viewed as a set of linear equations, one per zone. With only local zone calculations the equations can be put in the form
A_jϕ^k+1_j−1 + B_jϕ^k+1_j + C_jϕ^k+1_j+1 = D_j for 0 < j < 100.
and now we solve for the ϕs. We establish boundary conditions by saying that:
ϕ^k₀ = ϕ^k₁₀₀ = 0 for k > 0
and that ϕ⁰_j for 0 < j < 100 are provided as initial data somehow. At each time step we now have 99 unknowns.

This may seem daunting but the attentive first year algebra student who can program and has understood the first lecture on solving simultaneous equations is likely to come up with the best solution. Those who have also studied linear transformations and matrices may do much worse.

We are taught in first year algebra to eliminate unknowns from a set of equations, one at a time. The natural and strategic first candidate to eliminate is ϕ¹₁. We take equation A₁ϕ¹₀ + B₁ϕ¹₁ + C₁ϕ¹₂ = D₁ and solve for ϕ¹₁ and get:
ϕ¹₁ = (D₁−A₁ϕ¹₀−C₁ϕ¹₂)/B₁.

Note that ϕ¹₀ is known from the boundary conditions. Of the remaining equations only the next, with j=2, involves ϕ¹₁. We solve that equation for ϕ¹₂ and each time there is just one equation into which to substitute the new expression. We have luckily avoided the quadratic explosion that accompanies this elimination procedure in the general case.

When we have solved the last equation for ϕ¹₉₉ there are no unknowns in the resulting expression (ϕ¹₁₀₀ is known from boundary conditions) from which we have a concrete value for ϕ¹₉₉. The previous equation which expressed ϕ¹₉₈ in terms of ϕ¹₉₉ now immediately provides ϕ¹₉₈. This pattern continues in a backwards sweep thru the mesh until we have ϕ¹_j for all j.

This method was old hat in 1955 when I learned of it by word of mouth. It is numerically well behaved for physics problems such as these. It effectively moves the temperature signal several zones in one time step. Here is another more complete description of this calculation.

If our original problem had been for a metal loop instead of a metal rod we would not had a natural initial unknown to eliminate. There would have remained a lingering unknown just to the left of our initially eliminated unknown. This would have been resolved when we came around to it from the other side. It turns out that topology rules in these matters. We were driven to simplicity towards a forwards sweep followed by a backwards sweep. The Naïve solution was indeed strategic! There is a way to understand these problems abstracted from the equations. Here is a way to understand the accuracy and stability of this scheme.

Parallel