Heat Diffusion in a Thin Metal Rod

The solution to the heat diffusion equation meets the Fourier series

What would happen if you heated a small section of an insulated metal rod and left it alone for a while? Our daily experience of heat diffusion allows us to predict that the temperature will smooth out until it becomes uniform. In a scenario of perfect insulation, the heat will remain in the metal forever.

That is a correct qualitative description of the phenomenon, but how to describe it quantitatively?

We consider the one-dimensional problem of a thin metal rod wrapped in an insulating material. The insulation prevents the heat from escaping the rod from the side, but the heat can flow along the rod axis.

You can find the code used in this story here.

The heat diffusion equation

The heat diffusion equation is a simple second-order differential equation in two variables:

x ∈ [0, L] is the position along the rod, t is the time, u(x, t) is the temperature, and α is the thermal diffusivity of the material.

What intuition can we obtain about the temperature evolution by examining the heat diffusion equation?

Equation (1) states that the local rate of temperature change is proportional to the curvature, i.e., the second derivative with respect to x, of the temperature profile.

Figure 1 shows a temperature profile with three sections. The first section is linear; the second section has a negative second derivative, and the third section has a positive second derivative. The red arrows show the rate of change in temperature along the rod.

If ever a steady state where ∂u/∂t = 0 is reached, the temperature profile will have to smooth out up to the point where the temperature profile is linear.

The solution to the heat diffusion equation

The solution¹ to the heat diffusion equation (1) is:

You can verify by differentiating (2) that it does satisfy the differential equation (1). For those interested in the derivation, see Annex I.

The coefficients {Aₙ}, {Bₙ}, {λₙ}, C, D, and E are constants that must be fit from the initial and boundary conditions of the case. The work we did studying the Fourier series will play!

Boundary conditions

The boundary conditions are the constraints imposed at x=0 and x=L. We encounter two types of constraints in practical scenarios:

• Insulation, which translates into ∂u/∂x=0 at the rod extremity. This constraint prevents the heat from flowing in or out of the rod;
• Fixed temperature at the rod extremity: for example, the rod tip could be heated or cooled by a thermoelectric cooler, keeping it at a desired temperature.

The combination of constraint types will dictate the appropriate flavor of the Fourier series to represent the initial temperature profile.

Both ends insulated

When both rod ends are insulated, the gradient of the temperature profile gets set to zero at x=0 and x=L:

The initial condition is the temperature profile along the rod at t=0. Assume that for some obscure reason — perhaps the rod was possessed by an evil force — the temperature profile looks like this:

To run our simulation of the temperature evolution, we need to match equation (2) evaluated at t=0 with this function. We know the initial temperature profile through sample points but not its analytical expression. That is a task for a Fourier series expansion.

From our work on the Fourier series, we observed that an even half-range expansion yields a function whose derivative is zero at both extremities. That is what we need in this case.

Figure 3 shows the even half-range expansion of the function from Figure 2:

Although the finite number of terms used in the reconstruction creates some wiggling at the discontinuities, the derivative is zero at the extremities.

Equating equations (4), (5), (6), and (7) with equation (2) evaluated at t=0:

We can solve the constants:

Take a closer look at (14). This expression states that λₙ is proportional to the square of n, which is the number of half-periods that a particular cosine term goes through in the range [0, L]. In other words, n is proportional to the spatial frequency. Equation (2) includes an exponential factor exp(λₙt), forcing each frequency component to dampen over time. Since λₙ grows like the square of the frequency, we predict that the high-frequency components of the initial temperature profile will get damped much faster than the low-frequency components.

Figure 4 shows a plot of u(x, t) over the first second. We observe that the higher frequency component of the right-hand side disappears within 0.1 s. The moderate frequency component in the central section considerably fades but is still visible after 1 s.

When the simulation is run for 100 seconds, we get an almost uniform temperature:

Both ends at a fixed temperature

With both ends kept at a constant temperature, we have boundary conditions of the form:

The set of Fourier series that we studied in the previous post didn’t include the case of boundary temperatures fixed at non-zero values. We need to reformulate the initial temperature profile u₀(x) to develop a function that evaluates 0 at x=0 and x=L. Let us define a shifted initial temperature profile û₀(x):

The newly defined function û₀(x) linearly shifts the initial temperature profile u₀(x) such that û₀(0) = û₀(L) = 0.

As an illustration, Figure 6 shows an arbitrary initial temperature profile u₀, with set temperatures of 30 at x=0 and 70 at x=0.3. The green line (Cx + D) goes from (0, 30) to (0.3, 70). The orange curve represents û₀(x) = u₀(x) — Cx — D:

The shifted initial temperature profile û₀(x), going through zero at both ends, can be expanded with odd half-range expansion:

Equating equation (2) with (17), (18), (19), (20), and (21):

We can solve the constants:

The simulation of the temperature profile over time u(x, t) can now run, from equation (2):

In a permanent regime, the temperature profile is linear between the two set points, and constant heat flows through the rod.

Insulation at the left end, fixed temperature at the right end

We have these boundary conditions:

We follow essentially the same procedure as before. This time, we model the initial temperature profile with an even quarter-range expansion to get a zero derivative at the left end and a fixed value at the right end:

Which leads to the following constants:

The simulation over 1000 seconds shows the expected behavior. The left-hand extremity has a null temperature gradient, and the right-hand extremity stays at constant temperature. The permanent regime is a rod at a uniform temperature:

Conclusion

We reviewed the problem of the temperature profile dynamics in a thin metal rod. Starting from the governing differential equation, we derived the general solution.

We considered various boundary configurations. The boundary scenarios led us to express the initial temperature profile according to one of the Fourier series flavors we derived in the previous post. The Fourier series expression of the initial temperature profile allowed us to solve the integration constants and run the simulation of u(x, t).

Thank you for your time. You can experiment with the code in this repository. Let me know what you think!

¹ If a term is missing, please let me know in the comments.

Annex I

We want to demonstrate that the solution to the heat diffusion equation

is:

Let’s first acknowledge that, if a function u*(x, t) satisfies (A1.1), then the function u*(x, t) + Cx + D + Et + E/(2α) x² also satisfies (A1.1)

Proof:

Therefore, the general solution must include terms in Cx + D + Et + E/(2α) x².

Now, the leap of faith: the hypothesis of separability.

Let us assume that a solution u(x, t) has the following form:

Why make such an assumption?

Because it will make the solution easier to find. The justification will appear a posteriori if we can find a valid solution. In this case, we don’t risk running into erroneous conclusions based on a false assumption since we can always differentiate the found solution and check whether it satisfies (A1.1).

Inserting (A1.4) and (A1.5) into (A1.1):

In equation (A1.9), we observed that we equated a function of t with a function of x. The only way to satisfy this equation is to make both functions constant. Hence, we introduced the constant λₙ that must match both expressions.

(a) Solving for Gₙ(t), from (A1.9):

Inserting (A1.11) and (A1.12) into (A1.10):

Inserting (A1.15) into (A1.11):

(b) Solving for Fₙ(x), from (A1.9):

Inserting (A1.18) and (A1.20) into (A1.17):

Inserting (A1.24) into (A1.18):

We assume that Fₙ(x) is unique and the +- sign inside the sine function gets absorbed by the constant Bₙ.

Inserting (A1.16) and (A1.26) into the general solution:

Without loss of generality, we can set G₀=1 and let it get absorbed by the constants Aₙ and Bₙ:

□

Heat Diffusion in a Thin Metal Rod was originally published in Towards Data Science on Medium, where people are continuing the conversation by highlighting and responding to this story.