One-dimensional heat-flow simulations using simple digital and analogue models

The basic physical principles of heat flow in one dimension have long been established, and at a school or college level are first encountered in classic experiments such as the determination of thermal conductivity by the famous Searle's bar method. The conduction of heat down a perfectly insulated metal rod is an example of thermal diffusion .

As is well known, this process is described by the partial differential equation below:

$$\begin{equation}\frac{{\partial \theta }}{{\partial t\,}} = \,\alpha .\frac{{{\partial ^2}\theta }}{{\partial {x^2}}}\end{equation} \tag{ 1 }$$

where θ is temperature, t is time, and x is distance along the conductor. The quantity α is called the thermal diffusivity and is a property of the material of the conductor.

We define α by

$$\begin{equation}\alpha \equiv \frac{\lambda }{{\rho S}}\end{equation} \tag{ 2 }$$

Here $\lambda \,$ is the thermal conductivity, $\rho$ is the density, and S is the specific heat capacity of the material.

Exact solutions of equation (1) are, of course, known but the methods involved in their derivation are non-trivial and require a high degree of mathematical sophistication. At an introductory level, students of physics need to understand heat flow in a simple way which avoids unnecessary mathematical complexities. This can be achieved by using a basic step-wise approach which, with the aid of a computer, can provide approximate but useful solutions.

Figure 1(a) shows the one-dimensional heat flow situation to be modelled. On the left is a heat source of fixed temperature. At the other end of the thermally insulated conductor is a heat sink of fixed lower temperature. Students measuring thermal conductivity typically use a heat source at 100 °C and a heat sink at 0 °C. They are traditionally told to allow the temperature gradient down the bar to reach a steady state. Only then can the temperatures of two inserted thermometers be read and used in their calculations.

Zoom In Zoom Out Reset image size

The question which arises is 'how long is needed to reach the required steady state with a linear temperature profile down the bar?'

A step-wise analysis can answer this question.

The basis of such an analysis is shown in figure 1(b) where the conducting bar, of cross-sectional area A, is split into sections each of length h. Here we show three neighbouring sections with corresponding temperatures θ _{j −1}, θ _j and θ _{j +1} where j is an integer subscript. The problem is to find the values of θ as a function of distance down the bar.

At time t there will be a temperature gradient G_in between the left and central sections and a temperature gradient G_out between the central and right-hand section.

So, we can write the approximations:

$$\begin{equation}{G_{{\text{in }}}} \approx { }\frac{{{\theta _{j - 1}} - { }{\theta _j}}}{h}\end{equation} \tag{ 3 }$$

and

$$\begin{equation}{G_{{\text{out }}}} \approx { }\frac{{{\theta _j} - { }{\theta _{j + 1}}}}{h}.\end{equation} \tag{ 4 }$$

Also, the heat flux into the central element $= {{\lambda }}A\,{G_{{\text{in }}}}$ and the heat flux out of the central element ${{ = }}\,{{\lambda }}A\,\,{G_{{\text{out }}}}$

So, the energy absorbed by the central element in interval Δt is given by

$$\begin{equation}\Delta q\,{{ = }}\,{{\lambda }}A\left( {{G_{{\text{in }}}} - {G_{{\text{out}}}}} \right) \cdot \,\Delta t.\end{equation} \tag{ 5 }$$

But

$$\begin{equation*}\Delta q\, = m.S.\,\,\Delta \theta \end{equation*}$$

where m is the mass of each section and S is the specific heat capacity of the material of the bar. This can be written

$$\begin{equation}\Delta q = \rho Ah.S.\,\Delta \theta .\end{equation} \tag{ 6 }$$

Combining (3), (4), (5) and (6) leads to

$$\begin{equation}\Delta \theta \approx k\left( {{\theta _{j - {\mathbf{1}}{ }}} + {\theta _{j + {\mathbf{1}}}} - 2{\theta _j}} \right)\Delta t\end{equation} \tag{ 7 }$$

where $k \equiv \,\frac{\lambda }{{\rho S{h^2}}}$ = $\frac{\alpha }{{\,{h^2}}}$ (using equation (2)).

We thus have a means of performing a step-wise update of the temperature θ_j :

$$\begin{equation}{\theta _{j{ }}} \leftarrow {\theta _{j{ }}} + k\left( {{\theta _{j - {\mathbf{1}}{ }}} + {\theta _{j + {\mathbf{1}}}} - 2\,{\theta _j}} \right)\Delta t.\end{equation} \tag{ 8 }$$

The freely downloadable package SMath Studio [1, 2] provides a very simple and pedagogically transparent means of realising the above model on a computer.

For example, let us imagine the system in figure 1(a) which consists of seven thermal elements. This is depicted in figure 2.

Zoom In Zoom Out Reset image size

The values θ1−θ7 are the temperatures at the centres of each element, the end elements being the heat source and heat sink respectively. Elements 2–6 are the five sections into which the bar of length L is divided. Figure 3(a) depicts a SMath worksheet showing how the algorithm represented by expression (8) can be implemented. The length of the copper bar is 125 mm. This value has been chosen for reasons which will become clear later. As the bar is divided into five sections, the value of h is 25 mm.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The initial temperatures are placed in the row vector θ while the corresponding positions are placed in the row vector x. These vectors are changed into column vectors using the transpose functions θ: = θ^T and x: = x^T.

Expression (8) is implemented in the two nested for-loops and the values of x along with the updated values of θ are combined into a single matrix M using the augment function.

The theoretical linear steady-state function y(X) is also evaluated using the values of θ₁, θ₇, L and h.

The versatile SMath X–Y plot is supplied with the matrix M and the function y(X). Left clicking on the graph area brings up a formatting window where plotting types, axes, and colours can be selected. The result shown is the calculated temperature profile after a run time t of 10 s. Students and/or lecturers can easily experiment with different values of t and time step Δt. It is also possible to experiment with the behaviour of materials other than copper by changing the values of λ, ρ, and S.

The dashed red line represents the steady-state temperature profile.

Once confidence in using SMath is gained, it is then instructive to change other parameters in the system such as the bar length and number of bar sections. Because no cumbersome programming is needed, SMath enables the user to concentrate almost exclusively on the physics.

An interesting extension is to modify the worksheet to produce a plot showing temperature profiles for different run times. This is very straightforward and is achieved by means of a third outer for-loop which changes the values of t.

It is left as an exercise for interested readers to produce the curves in figure 3(b) with run times of 4 s, 10 s, and 16 s, and twice as many sections.

It can be seen that the longer run times produce curves which move closer and closer to the linearity of the steady state (red dashed line).

A useful analogue of one-dimensional heat flow is the RC ladder.

Figure 4 shows part of such a resistor-capacitor chain. It should be compared with figure 1(b). Considering the currents into and out of the central section, we have

$$\begin{equation}{I_{{\text{in}}}} = \frac{{{V_{j - 1}} - {V_j}}}{R}\end{equation} \tag{ 9 }$$

Zoom In Zoom Out Reset image size

and

$$\begin{equation}{I_{{\text{out}}}} = \frac{{{V_j} - {V_{j + 1}}}}{R}\end{equation} \tag{ 10 }$$

where V is electric potential and R is resistance.

The current into the capacitor is just

$$\begin{equation}I \approx \Delta Q/\Delta t\end{equation} \tag{ 11 }$$

and so, from (9), (10) and (11), we get

$$\begin{equation}\Delta Q \approx \left( {\frac{1}{{R\,}}} \right).\left( {{V_{j - 1}} + {V_{j + 1}} - 2{V_j}} \right).\,\Delta t\end{equation} \tag{ 12 }$$

where $\Delta Q$ is the charge flowing into the capacitor in time $\Delta t.$

But the change in potential

$$\begin{equation}\Delta V \approx \frac{{\Delta Q}}{C}\end{equation} \tag{ 13 }$$

Hence, from (12) and (13), we obtain

$$\begin{equation}\Delta V \approx \left( {\frac{1}{{RC\,}}} \right).\left( {{V_{j - 1}} + {V_{j + 1}} - 2{V_j}} \right).\,\Delta t\end{equation} \tag{ 14 }$$

It can be seen that equation (14) is the exact electrical analogue of the thermal equation (7). Instead of the constant k, we have the constant (1/RC).

It should be obvious from the dimensionality of equation (7) that the constant k has the unit s⁻¹. Similarly equation (14) shows that the unit of (1/RC) is also s⁻¹.

Hence, we could define the quantities τ_T and τ_E, the thermal and electrical time constants :

$$\begin{equation*}{\tau _T} \equiv \frac{1}{k}\,\,\,\,{\text{and}}\,\,\,\,{\tau _E} \equiv {{RC}}.\end{equation*}$$

These parameters determine the time taken for a given thermal or voltage profile to be reached.

If we use C = 470 μF and R = 10 kΩ, then ${\tau _E}$ is approximately 5 s. It will now be seen why the choice of h = 25 mm was made as this makes ${\tau _T}$ approximately the same i.e. 5 s. If ${\tau _T}$ ≈ ${\tau _E}$, we can use the RC ladder to directly mimic the behaviour of heat flow down our insulated bar. We can also arrange for a supply voltage of 5 V to represent a source temperature of 100 °C as used in our SMath model.

The above procedure was adopted because a limited number of 470 μF capacitors and 10 kΩ resistors were to hand and so it was easier to adjust h in order to match the two methods. A more professional method would be to (a) select a conductor of length L, (b) decide on the number of temperature-measurement points N, (c) obtain h from L/N, (d) calculate ${\tau _{T,}}$ and then select precision components with C and R values to give the same time constant.

Figure 5 shows the system I have used as an electric analogue of heat flow.

Zoom In Zoom Out Reset image size

Figure 6 shows the circuit board onto which the RC ladder is soldered. The voltages appearing on the capacitors are communicated to the analog ports A1–A5 on the Arduino which is connected to the PC.

Zoom In Zoom Out Reset image size

The PC is loaded with the Arduino sketch/program (see appendix) along with MakerPlot ¹ software to graph the voltages on five channels. The MakerPlot interface selected is the bars interface . This displays the voltage values (scaled to represent temperature as a bar graph. The capacitors are initially discharged and the input to the ladder connected to a 5 V power supply at the same time as the MakerPlot software is activated (red button below the display).

The bar graph shows the temperature profile evolving in real time over a period of maybe half a minute.

Figures 7 and 8 show typical results for 10 s and the steady state respectively.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The bar graphs and corresponding temperature values show pleasing agreement with the SMath results in figure 3(a).

I have shown how the physics of one-dimensional heat flow can be modelled in both a simple step-wise fashion and by using an equally simple electrical analogue. The use of SMath Studio makes the analysis easily comprehensible to students without the need for laborious computer programming. Equally, the electrical analogue is readily constructed in the teaching laboratory and demonstrates yet another pedagogical use for the versatile Arduino microcontroller. A comparison of the two approaches shows the power of analogues in physics and should form an important part of learning to do physics.

Both approaches can be easily modified and extended to deal with different circumstances. For example, an interesting exercise would be to arrange for the starting temperature distribution to be other than the one chosen above, e.g. starting with the bar at room temperature.

Once again, I should like to thank Andrey Ivashov (the creator of SMath) for giving his permission to publish SMath screen-shots.

My appreciations also go to Simon Gray for his expert, critical, and constructive comments.

//Heat-flow RC ladder analog with 5 channels Mk2 Nov. 2020

int N2,N3,N4,N5,N6; // declare N to be an integer (in range 0-1023)

float V2,V3,V4,V5,V6; // declare V to be a floating-point number (in range 0-5)

float TH2,TH3,TH4,TH5,TH6; // temperatures

float TH1=100; // source temperature

float G=0.0048876; // conversion factor 5/1023

float f=20;// conversion factor volts to degrees C. 5 volt = 100 deg C

float c=f*G; // composite conversion factor

void setup()

{

Serial.begin(9600); // set baud rate to 9600

}

void loop()

{

N2=analogRead(A1); // compute temperatures

TH2=c*N2;

N3=analogRead(A2);

TH3=c*N3;

N4=analogRead(A3);

TH4=c*N4;

N5=analogRead(A4);

TH5=c*N5;

N6=analogRead(A5);

TH6=c*N6;

Serial.print(TH1); // output temperatures

Serial.print(",");

Serial.print(TH2);

Serial.print(",");

Serial.print(TH3);

Serial.print(",");

Serial.print(TH4);

Serial.print(",");

Serial.print(TH5);

Serial.print(",");

Serial.println(TH6);

}