# Numerical Solution of 1-D Unsteady Conduction

## Discretization of Governing Equations

For 1-D transient conduction problems, the energy equation can be expressed as:

$\frac{1}{{A(x)}}\frac{d}{{dx}}\left( {kA(x)\frac{{dT}}{{dx}}} \right) + S = \rho c\frac{{\partial T}}{{\partial t}} \qquad \qquad(1)$

Similar to the case of steady-state conduction, it is assumed that the source term is a linear function of temperature, S = SC + SPT. Substituting this linearized source term into eq. (1) and multiplying the resulting equation by A(x) yields

Control volume for one-dimensional heat conduction
$\frac{d}{{dx}}\left( {kA(x)\frac{{dT}}{{dx}}} \right) + ({S_C} + {S_P}T)A(x) = \rho cA(x)\frac{{\partial T}}{{\partial t}}$

Integrating the above equation with respect to t in the interval of (t, t+Δt) and with respect to x over the control volume P (shaded area in the figure), we have

$\begin{array}{l} \int_t^{t + \Delta t} {\int_e^w {\frac{d}{{dx}}\left( {kA(x)\frac{{dT}}{{dx}}} \right)} dxdt} + \int_t^{t + \Delta t} {\int_e^w {({S_C} + {S_P}T)A(x)} dxdt} \\ = \int_e^w {\int_t^{t + \Delta t} {\rho cA(x)\frac{{\partial T}}{{\partial t}}dt} dx} \\ \end{array} \qquad \qquad(2)$

The first integral on the left-hand side becomes

$\begin{array}{l} \int_t^{t + \Delta t} {\int_e^w {\frac{d}{{dx}}\left( {kA(x)\frac{{dT}}{{dx}}} \right)} dxdt} = \int_t^{t + \Delta t} {\left[ {{{\left( {kA\frac{{dT}}{{dx}}} \right)}_e} - {{\left( {kA\frac{{dT}}{{dx}}} \right)}_w}} \right]} dt \\ = \int_t^{t + \Delta t} {\left[ {{k_e}{A_e}\frac{{{T_E} - {T_P}}}{{{{(\delta x)}_e}}} - {k_w}{A_w}\frac{{{T_P} - {T_W}}}{{{{(\delta x)}_w}}}} \right]} dt \\ \end{array}$
Variation of temperature within one time step

which requires knowledge of how temperature varies in the time interval of (t,t + t). There are commonly three kinds of choices (see figure on the right): (1) ${T_P}(t) \equiv T_P^t$ where $T_P^t$ is the temperature of grid point P at time t, (2) ${T_P}(t) \equiv T_P^{t + \Delta t}$ where $T_P^{t + \Delta t}$ is the temperature of grid point P at time t + t, or (3) TP(t) linearly varies from $T_P^t$ at timetto $T_P^{t + \Delta t}$ at time t + t. The following result is valid for all three variations:

$\int_t^{t + \Delta t} {{T_P}(t)} dt = \left[ {fT_P^{t + \Delta t} + (1 - f)T_P^t} \right]\Delta t$

where the above three cases correspond to f = 0, f = 1, and f = 1 / 2. Similar expressions for TE and TW can be obtained and we have

$\begin{array}{l} \int_t^{t + \Delta t} {\int_e^w {\frac{d}{{dx}}\left( {kA(x)\frac{{dT}}{{dx}}} \right)} dxdt = \left\{ {f\left[ {{k_e}{A_e}\frac{{{T_E} - {T_P}}}{{{{(\delta x)}_e}}} - {k_w}{A_w}\frac{{{T_P} - {T_W}}}{{{{(\delta x)}_w}}}} \right]} \right.} \\ {\rm{ }}\left. { + (1 - f)\left[ {{k_e}{A_e}\frac{{T_E^0 - T_P^0}}{{{{(\delta x)}_e}}} - {k_w}{A_w}\frac{{T_P^0 - T_W^0}}{{{{(\delta x)}_w}}}} \right]} \right\}\Delta t \\ \end{array}$

where the superscript t is replaced by 0 and the superscript t+Δt has been dropped, i.e., temperatures for grid point P at time t and t + t are represented by $T_P^0$ and at time t + t is represented by TP, respectively.

The second integral on the left-hand side of eq. (2) can be expressed as

$\int_t^{t + \Delta t} {\int_e^w {({S_C} + {S_P}T)A(x)} dxdt} = \left\{ {{S_C} + {S_P}[fT_P^0 + (1 - f){T_P}]} \right\}{A_P}{(\Delta x)_P}\Delta t$

The integral on the right hand of eq. (2) can be evaluated as

$\int_e^w {\int_t^{t + \Delta t} {\rho cA(x)\frac{{\partial T}}{{\partial t}}dt} dx} = {(\rho c)_P}{A_P}({T_P} - T_P^0){(\Delta x)_P}$

Substituting the above equations into eq. (2) the discretized equation for one-dimensional transient conduction becomes

$\left\{ {f\left[ {{k_e}{A_e}\frac{{{T_E} - {T_P}}}{{{{(\delta x)}_e}}} - {k_w}{A_w}\frac{{{T_P} - {T_W}}}{{{{(\delta x)}_w}}}} \right]} \right.\left. { + (1 - f)\left[ {{k_e}{A_e}\frac{{T_E^0 - T_P^0}}{{{{(\delta x)}_e}}} - {k_w}{A_w}\frac{{T_P^0 - T_W^0}}{{{{(\delta x)}_w}}}} \right]} \right\}\Delta t$
$+ \left\{ {{S_C} + {S_P}[fT_P^0 + (1 - f){T_P}]} \right\}{A_P}{(\Delta x)_P}\Delta t = {(\rho c)_P}{A_P}({T_P} - T_P^0){(\Delta x)_P} \qquad \qquad(3)$

which can be rearranged to obtain (Tao, 2001)

${a_P}{T_P} = {a_E}[f{T_E} + (1 - f)T_E^0] + {a_W}[f{T_W} + (1 - f)T_W^0] + b \qquad \qquad(4)$

where

${a_E} = \frac{{{k_e}{A_e}}}{{{{(\delta x)}_e}}},{\rm{ }}{a_W} = \frac{{{k_w}{A_w}}}{{{{(\delta x)}_w}}} \qquad \qquad(5)$
${a_P} = f{a_E} + f{a_W} + a_P^0 - f{S_P}{A_P}{(\Delta x)_P} \qquad \qquad(6)$
$b = \left[ {a_P^0 - (1 - f){a_E} - (1 - f){a_W} + (1 - f){S_P}{A_P}{{(\Delta x)}_P}} \right]T_P^0 + {S_C}{A_P}{(\Delta x)_P} \qquad \qquad(7)$
$a_P^0 = \frac{{\rho c{{(\Delta x)}_P}}}{{\Delta t}} \qquad \qquad(8)$

## Explicit and Implicit Schemes

Different schemes can be obtained when f in eqs. (4) – (8) takes different values. When f = 0, i.e., the temperature at grid point P in the entire interval of (t,t + t) is equal to the temperature at time t, eq. (4) becomes

${a_P}{T_P} = {a_E}T_E^0 + {a_W}T_W^0 + b \qquad \qquad(9)$

which indicates that the temperature for grid point P at time t + Δt can be explicitly obtained from temperatures at the previous time step, because it is not related to temperatures of grid points E and W at time t + t. The scheme for the case of f = 0 is referred to as the explicit scheme. While an explicit scheme is computationally very convenient, it is not unconditionally stable. For the case of constant thermophysical properties, without an internal heat source, and with uniform grid size, eq. (9) becomes

${T_P} = \frac{{\alpha \Delta t}}{{{{(\Delta x)}^2}}}(T_E^0 + T_W^0) + \left[ {1 - \frac{{2\alpha \Delta t}}{{{{(\Delta x)}^2}}}} \right]T_P^0 \qquad \qquad(10)$

Although the criterion of the stability for eq. (10) can be rigorously derived from Von Neumann analysis, we will use a simple reasoning to obtain this criterion. If the temperatures at E and W remain unchanged, we would expect that TP is higher if $T_P^0$ is higher and TP is lower if $T_P^0$ is lower. This can be assured only if the coefficient of $T_P^0$ in eq. (10) is positive, i.e.,

${\rm{F}}{{\rm{o}}_\Delta } = \frac{{\alpha \Delta t}}{{{{(\Delta x)}^2}}} < \frac{1}{2} \qquad \qquad(11)$

where FoΔ is grid Fourier number. Equation (11) indicates that for a system with uniform grid size of Δx, the time step Δt must be small enough to ensure the grid Fourier number is less than 0.5. Equation (10) can be rewritten in terms of the grid Fourier’s number as

${T_P} = {\rm{F}}{{\rm{o}}_\Delta }(T_E^0 + T_W^0) + (1 - 2{\rm{F}}{{\rm{o}}_\Delta })T_P^0 \qquad \qquad(12)$

If f = 0.5, i.e., the temperature at grid point P varies linearly in the interval of (t,t + t). Equation (4) becomes

${a_P}{T_P} = \frac{{{a_E}}}{2}[{T_E} + T_E^0] + \frac{{{a_W}}}{2}[{T_W} + T_W^0] + b \qquad \qquad(13)$

which is referred to as the Crank-Nicolson (C-N) scheme. Although the C-N scheme has often been described as unconditionally stable, its result is physically meaningful only if the grid Fourier number is less than 1. For the case that FoΔ > 1, C-N scheme is mathematically stable but the solution physically does not make sense. When f = 1, eqs. (4), (6) and (7) become

${a_P}{T_P} = {a_E}{T_E} + {a_W}{T_W} + b \qquad \qquad(14)$
${a_P} = {a_E} + {a_W} + a_P^0 - {S_P}{A_P}{(\Delta x)_P} \qquad \qquad(15)$
$b = \left[ {a_P^0 + {S_P}{A_P}{{(\Delta x)}_P}} \right] + {S_C}{A_P}{(\Delta x)_P} \qquad \qquad(16)$

which is referred to as a fully-implicit scheme. It is truly unconditionally stable and physically meaningful under any grid Fourier number. This advantage makes the fully-implicit scheme top choice for numerical solution of conduction problems. Since fully-implicit schemes are stable under any time step, one can set $\Delta t \to \infty$ and the discretized equations for transient conduction become identical to those for the steady-state conduction problem. Thus, we can run the program for transient heat conduction with implicit schemes for very large time steps, say 1020, to get a solution for a steady state conduction problem. Since the temperature of grid point P at time t + t is related to the temperatures of its neighbor points E and W, the TDMA discussed in the previous subsection can be employed.

## References

Tao, W.Q., 2001, Numerical Heat Transfer, 2nd Ed., Xi’an Jiaotong University Press, Xi’an, China (in Chinese).

Faghri, A., Zhang, Y., and Howell, J. R., 2010, Advanced Heat and Mass Transfer, Global Digital Press, Columbia, MO.