# 1-Dimensional Transient Heat Conduction

One famous example on solving Parabolic Ordinary Differential Equations (ODE) using finite differences is 1-D transient heat conduction which would be described in this section.

# Diffusion problem - Implementation

An insulated rod of length 1 meter initially has a temperature of $T(x,0)=20 ° C$. At $t=0$ the rod is placed between two hot walls of the same fixed temperature $T=100 ° C$ (Fig. b.3.1). The problem is to numerically find the subsequent temperature $T(x,t)$ of any point in the rod.

Fig. b.3.1: Unsteady heat conduction in a rod

The transient heat conduction in one dimension is governed by:
$$\frac{\partial T}{\partial t} - \sigma \frac{\partial^2 T}{\partial x^2} = 0,$$ where $T$ is the temperature and $\sigma$ is the thermal diffusivity. By means of finite difference expressions, it is possible to discretize the equation to the forward time, centered space (FTCS) scheme:
$$\frac{T_j^{n+1} - T_j^n}{\Delta t} - \sigma \frac{\left( T_{j-1}^n - 2T_j^n + T_{j+1}^n \right)}{\Delta x^2} = 0,$$ with step sizes $\Delta t$ ($n = 0, 1, ..., N$) and $\Delta x$ ($j = 1, 2, ..., J$).

# Initial and boundary conditions

For the solution of the diffusion equation, the initial conditions:
$$T(x,0) = T_0(x) = 20 ° C \ \ , \ \ 0 \leq x \leq 1,$$ serve to initialize the dependent variable $T_j^0$, and the boundary values $T_1^{n+1}$ and $T_J^{n+1}$ are given by the boundary conditions: $$T(0,t) = T_A = 100 ° C \ \ , \ \ T(1,t) = T_B = 100 ° C.$$

# Numerical and exact solutions

Rewriting the Eq. (2) using $s = \sigma \Delta t/\Delta x^2$, keeping the term $T_j^{n+1}$ in the left-hand side, and transferring the rest to the right-hand side (RHS), we have:
$$T_j^{n+1} = sT_{j-1}^n + (1-2s)T_j^n + sT_{j+1}^n.$$ Now considering the initial and boundary conditions, a loop to solve the temperature value at each node can be employed. For simplicity, these assumptions are considered to solve the diffusion ODE at the first step, Case (1-1): $\Delta x = 0.1$, $\Delta t = 500$, $\sigma = 10^{-5}$, and $t_{max} = 3000$; therefore, $s = 0.5$. The solutions for each time step is shown in Fig. b.3.2.

Fig. b.3.2: Temperature distribution, Case (1-1)

The values of $T(0,0)$ and $T(1,0)$ can be averaged between their values implied by the initial conditions ($20 ° C$) and those implied by the boundary conditions ($100 ° C$). Implementing these changes at initial conditions results in Case (2-1) which is illustrated in Fig. b.3.3.

Fig. b.3.3: Temperature distribution, Case (2-1)

In order to check the accuracy of the results, the exact solution is needed which can be computed (obtained by the separation of variables approach) from: $$\text{TE}(x_j,t_n) = T_A - \sum\limits_{m=1}^{\text{MAX}} \left[\frac{4 (T_A - T_0)}{(2m-1)\pi}\sin \left[(2m-1)\pi x_j\right]e^{[-\alpha (2m-1)^2\pi ^2 t_n]}\right],$$ where MAX (=10) is the number of terms in the exact equation. Fig. b.3.4 compares the results of cases (1-1) and (2-1) at the last time-step with the obtained exact solution.

Fig. b.3.4: Temperature distribution at $t = 3000$, numerical vs. exact solutions

It is clear that the second case is more accurate even with such a course grid.

Moreover, the rms error could be obtained from:
$$\text{RMS} = \sqrt{\left.\sum\limits_{j=1}^{J_{max}} (T_j - \text{TE}_j)^2 \right/ J_{max}.}$$ Computing the rms error helps to compare two mentioned cases more accurately. This comparison will be done after computing the problem for different values of $s$.

# More accurate numerical solutions

Decreasing the value of $s$ results in a more accurate solution. But also the results depend on $\Delta x$ and $\Delta t$ variations. It means that the magnitude of error would change for a fixed value of $s$, but different combinations of $\Delta x$ and $\Delta t$. To show this, the second case (which seems to be more accurate than the first case) is modified as it is shown in Table b.3.1. A comparison of all cases is illustrated in this table.

Table b.3.1: RMS error of numerical solutions

As an example, comparing the rms error of cases (2-2) and (2-4) shows that although the value of $s$ is decreased, the error is larger. Because the grid is not fine enough. Refining the grid in case (2-5) leads to a much better results. The following figures show the temperature distribution for some cases together with the exact solution.

Fig. b.3.5: Case (2-2), temperature distribution

Fig. b.3.6: Case (2-3), temperature distribution

# Stability

Now it is easily possible to check the stability of the scheme by running the case (2-7) in which the value of $s>0.5$ makes it unstable. This is plotted in Fig. b.3.7 to show the unstable solution.

Fig. b.3.7: Case (2-7), unstable solution