Heat Diffusion Transient

This section discretizes the heat diffusion equation in three dimensions for any arbitrary shape mesh cell using the finite volume method (FVM) for transient conditions. The discretization is based on the general form of the heat diffusion equation, which includes both dependent and independent source terms.

Generic Form of the Heat Diffusion Equation

The heat diffusion equation in three dimensions can be written as:

\[\rho c_p \frac{\partial T}{\partial t} = \frac{\partial}{\partial x_i}(k \frac{\partial T}{\partial x_i}) + S_T\]

where:

  • \(T\) is the temperature (dependent variable).

  • \(\alpha\) is the thermal diffusivity (\(\alpha = \frac{k}{\rho c_p}\)) where

    • \(k\) is the thermal conductivity.

    • \(\rho\) is the density of the material.

    • \(c_p\) is the specific heat capacity.

  • \(S_T\) is a source term.

Discretization in 3D

Continuing from the steady state case

Example SVG

The equation now have the time derivative term which can be integrated over a time step \(\Delta t\) to yield the following integral form:

\begin{align} \int_{t}^{t+\Delta t} \int_{\Delta V} \rho c \frac{\partial T}{\partial t} dV dt & = \int_{t}^{t+\Delta t} \int_{\Delta V} \left[\frac{\partial}{\partial x}\left(k \frac{\partial T}{\partial x}\right) + S_T\right] dV dt \notag \\ \rho c (T_{i}-T_{i}^{0}) \Delta V & = \int_{t}^{t+\Delta t} \left( \left[kA \frac{T_{i+1} - T_i}{\Delta x} \right]_{i+1} \right. \notag \\ & \left. + \left[kA \frac{T_{i-1} - T_i}{\Delta x} \right]_{i-1} + S_u + S_i T_i \right)\ dt \end{align}

On the temperature integral of the right side of the equation,we can generalize the equation by means of a weighting parameter \(\theta\) which can be used to approximate the temperature at the next time step as follows:

\begin{align*} T_{i}^{n+1} = T_{i}^{n} + \theta \left( T_{i}^{n+1} - T_{i}^{n} \right) \end{align*}

where \(\theta\) is a weighting factor that can be set to 0 for implicit, 0.5 for the Crank-Nicolson method or 1 for the explicit method.

Dropping the superscript for the future time step and using \(n = 0\) for the current time step, we can express the equation as:

\begin{align*} T_{i} & = T_{i}^{0} + \theta \left( T_{i} - T_{i}^{0} \right) \\ T_{i} & = \theta T_{i} + (1-\theta) T_{i}^{0} \\ \int_{t}^{t+\Delta t} T_{i} dt & = \left[ \theta T_{i} + (1-\theta) T_{i}^{0} \right] \Delta t \end{align*}

Thus the equation (1) can be re-written as:

\begin{align} \rho c (T_{i}-T_{i}^{0}) \Delta V & = \theta\left(\left[kA \frac{T_{i+1} - T_i}{\Delta x} \right]_{i+1} + \left[kA \frac{T_{i-1} - T_i}{\Delta x} \right]_{i-1} + S_i T_{i}\right) \Delta t \notag \\ + (1-\theta) & \left(\left[kA \frac{T_{i+1}^{0} - T_i^{0}}{\Delta x} \right]_{i+1} + \left[kA \frac{T_{i-1}^{0} - T_i^{0}}{\Delta x} \right]_{i-1} + S_i T_{i}^{0} \right) \Delta t \notag \\ & + S_u \Delta V \Delta t \notag \\ \Rightarrow \rho c (T_{i}-T_{i}^{0}) \frac{\Delta V}{\Delta t} & = \theta\left(\left[kA \frac{T_{i+1} - T_i}{\Delta x} \right]_{i+1} + \left[kA \frac{T_{i-1} - T_i}{\Delta x} \right]_{i-1} + S_i T_{i}\right) \notag \\ + (1-\theta) & \left(\left[kA \frac{T_{i+1}^{0} - T_i^{0}}{\Delta x} \right]_{i+1} + \left[kA \frac{T_{i-1}^{0} - T_i^{0}}{\Delta x} \right]_{i-1} + S_i T_{i}^{0} \right) \notag \\ & + S_u \Delta V \end{align}

Rearrenging them to organize all the unknowns on the left side and knowns on the right side, we can express the equation as:

\begin{align} & - \theta \left[ \frac{kA}{\Delta x} \right]_{i-1} T_{i-1} \notag \\ & + \left[ \frac {\rho c \Delta V}{\Delta t} + \theta \left\{ \left( \frac{kA}{\Delta x} \right)_{i+1} + \left( \frac{kA}{\Delta x} \right)_{i-1} -S_i \right\} \right] T_{i} \notag \\ & -\theta \left[ \frac{kA}{\Delta x} \right]_{i+1} T_{i+1} \notag \\ & = (1-\theta) \left[ \frac{kA}{\Delta x} \right]_{i-1} T^{0}_{i-1}\notag \\ & + \left[ \frac {\rho c \Delta V}{\Delta t} - (1-\theta) \left\{ \left(\frac{kA}{\Delta x} \right)_{i+1} + \left( \frac{kA}{\Delta x} \right)_{i-1} -S_i \right\}\right] T_{i}^{0} \notag \\ & +(1-\theta) \left[ \frac{kA}{\Delta x} \right]_{i+1} T^{0}_{i+1} + S_u \Delta V \end{align}

This is the generic discretization equation for the transient heat diffusion equation using the finite volume method. The equation can be applied to any arbitrary shape mesh cell by appropriately defining the coefficients based on the geometry, cell connectivity, boundary faces, and material properties of the mesh cell.

If \(\theta = 0\), the equation becomes explicit, meaning that the temperature at the next time step is calculated directly from the current temperature and source terms. If \(\theta = 1\), it becomes implicit, requiring a system of equations to be solved at each time step. For \(\theta = 0.5\), it represents the Crank-Nicolson method, which is a time-centered scheme providing a balance between stability and accuracy. The equation is organized such that all the unknowns (temperatures at the next time step) are on the left side, while all known values (temperatures at the current time step and source terms) are on the right side.

For \(\theta = 0\), the equation simplifies to an explicit form as follows

\begin{align} & \left[ \frac {\rho c \Delta V}{\Delta t} \right] T_{i} \notag \\ & = \left[ \frac{kA}{\Delta x} \right]_{i-1} T^{0}_{i-1}\notag \\ & + \left[ \frac {\rho c \Delta V}{\Delta t} - \left\{ \left(\frac{kA}{\Delta x} \right)_{i+1} + \left( \frac{kA}{\Delta x} \right)_{i-1} -S_i\right\} \right] T_{i}^{0} \notag \\ & +\left[ \frac{kA}{\Delta x} \right]_{i+1} T^{0}_{i+1} + S_u \Delta V \end{align}

For \(\theta = 1\), the equation simplifies to an implicit form as follows

\begin{align} & - \left[ \frac{kA}{\Delta x} \right] T_{i-1} \notag \\ & + \left[ \frac {\rho c \Delta V}{\Delta t} + \left\{ \left( \frac{kA}{\Delta x} \right)_{i+1} + \left( \frac{kA}{\Delta x} \right)_{i-1} -S_i \right\} \right] T_{i} \notag \\ & -\left[ \frac{kA}{\Delta x} \right] T_{i+1} \notag \\ & = \left[ \frac {\rho c \Delta V}{\Delta t} \right] T_{i}^{0} \notag \\ & + S_u \Delta V \end{align}

For \(\theta = \frac{1}{2}\), the equation translates into a Crank-Nicolson form.

Most interestingly, the equation (3) can be used to derive the steady state heat diffusion equation by setting \(\Delta t \to \infty\) and \(\theta = 1\) , which leads to the steady state form of the heat diffusion equation.

\begin{align} & - \left[ \frac{kA}{\Delta x} \right]_{i-1} T_{i-1} \notag \\ & + \left[ \left( \frac{kA}{\Delta x} \right)_{i+1} + \left( \frac{kA}{\Delta x} \right)_{i-1} -S_i \right] T_{i} \notag \\ & -\left[ \frac{kA}{\Delta x} \right]_{i+1} T_{i+1} \notag \\ & = S_u \Delta V \end{align}

Boundary conditions

Considering the equation (2), we can estimate the discretization for the different types of boundary conditions.

1. Dirichlet Boundary Condition

For a Dirichlet boundary condition, the temperature at the boundary is specified. For example, if the temperature at the boundary is fixed at \(T_B\), then (2) becomes:

\begin{align*} \rho c (T_{i}-T_{i}^{0}) \frac{\Delta V}{\Delta t} & = \theta\left(\left[kA \frac{T_{i+1} - T_i}{\Delta x} \right]_{i+1} + \left[kA \frac{T_{B} - T_i}{\Delta x} \right]_{B} + S_i T_{i}\right) \notag \\ + (1-\theta) & \left(\left[kA \frac{T_{i+1}^{0} - T_i^{0}}{\Delta x} \right]_{i+1} + \left[kA \frac{T_{B}^{0} - T_i^{0}}{\Delta x} \right]_{B} + S_i T_{i}^{0} \right) \notag \\ & + S_u \Delta V \end{align*}

Rearrenging the equation to organize all the unknowns on the left side and knowns on the right side similar to (3), we can express the equation as:

\begin{align} & - \theta \left[ 0 \right]_B T_{i-1} \notag \\ & + \left[ \frac {\rho c \Delta V}{\Delta t} + \theta \left\{ \left( \frac{kA}{\Delta x} \right)_{i+1} + \left( \frac{kA}{\Delta x} \right)_B -S_i \right\} \right] T_{i} \notag \\ & -\theta \left[ \frac{kA}{\Delta x} \right]_{i+1} T_{i+1} \notag \\ & = (1-\theta) \left[ 0 \right]_B T^{0}_{i-1}\notag \\ & + \left[ \frac {\rho c \Delta V}{\Delta t} - (1-\theta) \left\{ \left(\frac{kA}{\Delta x} \right)_{i+1} + \left( \frac{kA}{\Delta x} \right)_{B} -S_i \right\}\right] T_{i}^{0} \notag \\ & +(1-\theta) \left[ \frac{kA}{\Delta x} \right]_{i+1} T^{0}_{i+1} \notag \\ & + S_u \Delta V + \theta \left[ \frac{k A T_B }{\Delta x} \right]_B + (1-\theta) \left[ \frac{k A T_B^0 }{\Delta x} \right]_B \notag \end{align}

Since the temperature at the boundary is fixed and known, meaning \(T_B = T_B^0\), we can express the equation as:

\begin{align} & - \theta \left[ 0 \right]_B T_{i-1} \notag \\ & + \left[ \frac {\rho c \Delta V}{\Delta t} + \theta \left\{ \left( \frac{kA}{\Delta x} \right)_{i+1} + \left( \frac{kA}{\Delta x} \right)_B -S_i \right\} \right] T_{i} \notag \\ & -\theta \left[ \frac{kA}{\Delta x} \right]_{i+1} T_{i+1} \notag \\ & = (1-\theta) \left[ 0 \right]_B T^{0}_{i-1}\notag \\ & + \left[ \frac {\rho c \Delta V}{\Delta t} - (1-\theta) \left\{ \left(\frac{kA}{\Delta x} \right)_{i+1} + \left( \frac{kA}{\Delta x} \right)_{B} -S_i \right\}\right] T_{i}^{0} \notag \\ & +(1-\theta) \left[ \frac{kA}{\Delta x} \right]_{i+1} T^{0}_{i+1} \notag \\ & + S_u \Delta V + \theta \left[ \frac{k A T_B^0 }{\Delta x} \right]_B + (1-\theta) \left[ \frac{k A T_B^0 }{\Delta x} \right]_B \notag \end{align}

Therefore the coefficients with the term \(\theta\) cancels out, and the equation simplifies to:

\begin{align} & - 0 \notag \\ & + \left[ \frac {\rho c \Delta V}{\Delta t} + \theta \left\{ \left( \frac{kA}{\Delta x} \right)_{i+1} + \left( \frac{kA}{\Delta x} \right)_B -S_i \right\} \right] T_{i} \notag \\ & -\theta \left[ \frac{kA}{\Delta x} \right]_{i+1} T_{i+1} \notag \\ & = 0\notag \\ & + \left[ \frac {\rho c \Delta V}{\Delta t} - (1-\theta) \left\{ \left(\frac{kA}{\Delta x} \right)_{i+1} + \left( \frac{kA}{\Delta x} \right)_{B} -S_i \right\}\right] T_{i}^{0} \notag \\ & +(1-\theta) \left[ \frac{kA}{\Delta x} \right]_{i+1} T^{0}_{i+1} \notag \\ & + S_i^u \Delta V + \left[ \frac{k A T_B^0 }{\Delta x} \right]_B \end{align}

2. Neumann Boundary Condition

The flux at the boundary is specified for the Neumann boundary condition, which can be expressed as:

  1. Thermal Insulation \(\frac {\partial T}{\partial x} = 0\)

  2. a convective heat transfer with coefficient \(h\) and ambient temperature \(T_\infty\) as \(\frac {\partial T}{\partial x} = hA(T_\infty - T_B)\),

  3. a radiative heat loss with emissivity \(\epsilon\), Stefan Boltzmann constant \(\epsilon\), and ambient temperature \(T_\infty\) as \(\frac {\partial T}{\partial x} = \sigma \epsilon A(T_\infty^4 - T_B^4)\).

We can consider the temperature of the boundary node and the cell almost equal as a mean to simplify where \(T_B = T_i\). In an attempt to develop a generic boundary condition that covers all three of the cases, we can add all of them together and replace the flux at the boundary. Then the equation (2) becomes:

\begin{align*} \rho c (T_{i}-T_{i}^{0}) \frac{\Delta V}{\Delta t} & = \theta\left(\left[kA \frac{T_{i+1} - T_i}{\Delta x} \right]_{i+1} \notag \right. \\ & \left. + \left[0 + hA(T_B-T_\infty) + \sigma \epsilon A(T_B^4 - T_\infty^4) \right]_{B} + S_i T_{i} \vphantom{\frac{\Delta V}{\Delta t}}\right) \notag \\ + (1-\theta) & \left(\left[kA \frac{T^{0}_{i+1} - T^{0}_i}{\Delta x} \right]_{i+1} \notag \right. \\ & \left. + \left[0 + hA(T_\infty - T_B) + \sigma \epsilon A(T_\infty^4 - T_B^4) \right]_{B} + S_i T_{i} \vphantom{\frac{\Delta V}{\Delta t}}\right) \notag \\ & + S_u \Delta V \end{align*}

Rearrenging the equation to organize all the unknowns on the left side and knowns on the right side similar to (3) where the known boundary values having the coefficient \(\theta\) cancels out and thus, we can express the equation as:

\begin{align} & -0 + \left[ \frac {\rho c \Delta V}{\Delta t} + \theta \left\{ \left( \frac{kA}{\Delta x} \right)_{i+1} -S_i \right\} \right] T_{i} \notag \\ & -\theta \left[ \frac{kA}{\Delta x} \right]_{i+1} T_{i+1} \notag \\ & = 0 + \left[ \frac {\rho c \Delta V}{\Delta t} - (1-\theta) \left\{ \left(\frac{kA}{\Delta x} \right)_{i+1} -S_i \right\}\right] T_{i}^{0} \notag \\ & +(1-\theta) \left[ \frac{kA}{\Delta x} \right]_{i+1} T^{0}_{i+1} \notag \\ & + S_i^u \Delta V + \left[ hA(T_\infty - T_B) + \sigma \epsilon A(T_\infty^4 - T_B^4) \right]_{B} \end{align}

Generic Discretization for any dimensions

In scenarios of higher dimensions, the cell connectivity is more than two thus the generic equation for cells without any boundary faces are a modified form of the equation (6) as follows:

\begin{align} & - \theta \sum_{j=1}^n \left[ \frac{kA}{\Delta x} \right]_{i, j} T_{j} \notag \\ & + \left[ \frac {\rho c \Delta V}{\Delta t} + \theta \left\{ \sum_{j=1}^n \left( \frac{kA}{\Delta x} \right)_{i, j} - S_i \right\} \right] T_{i} \notag \\ & = (1-\theta) \sum_{j=1}^n \left[ \frac{kA}{\Delta x} \right]_{i-1} T^{0}_{i-1}\notag \\ & + \left[ \frac {\rho c \Delta V}{\Delta t} - (1-\theta) \left\{ \sum_{j=1}^n \left(\frac{kA}{\Delta x} \right)_{i, j} - S_i \right\}\right] T_{i}^{0} \notag \\ & + S_i^u \Delta V \end{align}

The above equation works for any arbitrary shape mesh cell in any dimensions, where \(n\) is the number of shared cells (connected cells) with the cell indexed with \(i\). The coefficients are defined as follows:

\begin{align} a_{i,j} & = - \theta \sum_{j=1}^n \left[ \frac{kA}{\Delta x} \right]_{i, j} \notag \\ a_{i,i} & = \left[ \frac {\rho c \Delta V}{\Delta t} + \theta \left\{ \sum_{j=1}^n \left( \frac{kA}{\Delta x} \right)_{i, j} - S_i \right\} \right] \notag \\ b_i & = (1-\theta) \sum_{j=1}^n \left[ \frac{kA}{\Delta x} \right]_{i-1} T^{0}_{i-1}\notag \\ & + \left[ \frac {\rho c \Delta V}{\Delta t} - (1-\theta) \left\{ \sum_{j=1}^n \left(\frac{kA}{\Delta x} \right)_{i, j} - S_i \right\}\right] \notag \\ & + S_i^u \Delta V \end{align}

for a generic system of linear algebraic equation of the form:

\begin{align*} A_{i,j} X_{i} = b_{i} \end{align*}

where \(A\) is the coefficient matrix, \(X\) is the vector of unknown temperatures, and \(b\) is the vector of known values.

For a given cell if there are boundary faces, only then special equation is to be used. If it is a Dirichlet boundary conditions

\begin{align*} & \left[ \frac {\rho c \Delta V}{\Delta t} + \theta \left\{ \sum_{j=1}^n \left( \frac{kA}{\Delta x} \right)_{i, B_j} - S_i \right\} \right] T_{i} \notag \\ & = \left[ \frac {\rho c \Delta V}{\Delta t} - (1-\theta) \left\{ \sum_{j=1}^n \left(\frac{kA}{\Delta x} \right)_{i, B_j} - S_i \right\}\right] T_{i}^{0} \notag \\ & + S_i^u \Delta V + \sum_{j=1}^n \left[ \frac{k A T_{j,B}^0 }{\Delta x} \right]_B \end{align*}

The \(\frac{\rho c \Delta V}{\Delta t}\) term should be already added to the coefficient matrix \(A\) for the cell indexed with \(i\) and thus, it is unnecessary/redundent to include this term in the equation above. Therefore, the co-efficient update can be simplified to:

\begin{align} a_{i,i} & += \theta \left\{ \sum_{j=1}^n \left( \frac{kA}{\Delta x} \right)_{i, B_j} - S_i \right\} \notag \\ b_i & += - (1-\theta) \left\{ \sum_{j=1}^n \left(\frac{kA}{\Delta x} \right)_{i, B_j} - S_i \right\} T_i^0 \notag \\ & + \sum_{j=1}^n \left[ \frac{k A T_{j,B}^0 }{\Delta x} \right]_B \end{align}

For Neumann boundary conditions

\begin{align*} & \left[ \frac {\rho c \Delta V}{\Delta t} + \theta \left\{ - S_i \right\} \right] T_{i} \notag \\ & = \left[ \frac {\rho c \Delta V}{\Delta t} - (1-\theta) \left\{- S_i \right\}\right] T_{i}^{0} \notag \\ & + S_i^u \Delta V + \sum_{j=1}^n \left[ hA(T_\infty - T_{j,B}) + \sigma \epsilon A(T_\infty^4 - T_{j, B}^4) \right]_{B} \end{align*}

where

\begin{align} a_{i,i} & += -\theta S_i \notag \\ b_i & += (1-\theta) S_i T_i^0 \notag \notag \\ & + \sum_{j=1}^n \left[ hA(T_\infty - T_{j,B}) + \sigma \epsilon A(T_\infty^4 - T_{j, B}^4) \right]_{B} \end{align}

Equation (10), (11), and (12) are all we need for the construction of matrix \(A\) and the known vector \(b\) for FVM solution for any n-dimensional heat diffusion transient/steady state problem for any boundary conditions. We simply set \(\theta = 0\) for explicit, \(\theta = 1\) for implicit, and \(\theta = 0.5\) for Crank-Nicolson method and \(\Delta t = \infty\) with np.inf for steady state cases.