How to Solve the Heat Equation?

Transient problems in heat transfer are restricted to simple geometries and simple analytical solutions. Nevertheless, in many cases the geometry and boundary conditions make difficult the use of analytical techniques, therefore we can make use of the finite-difference methods.

Numerical method for the heat equation

For the two-dimensional system of Figure 1 under transient conditions with constant properties, the heat equation is:



Figure 1. Discretized computational domain with finite-differences.

We may use the central-difference approximations to the spatial derivates. The subscripts (i, j) are used to designate the locations of the  nodal points,in x and y respectively. The heat equation equation must be discretized in space and in time to obtain a numerical solution which depends on x, y and t. To achieve this purpose, the integer p is used to calculate the time:


Finite-difference approximation to the time derivate is represented as:


p denotes the time dependence of T,  the time derivative is expressed in terms of the difference in temperatures associated with the new temporal step (p+1)  and previous temporal step (p). Finite-difference approximation to the spatial derivate is represented as:

figura 2

The explicit form of the finite-difference equation for the interior node i, j is:


The last expresion is found just subtituting the explicit form of the spatial derivative in the heat equation. Then, solving for the nodal temperature T(i,j) at the new temporal step (p+1) and assuming that Δx= Δy, it follows that:


The last equation is explicit because unknown nodal temperatures for the new time are determined by known nodal temperatures at the previous time, and F0  is the Fourier number defined as:


The temperature of each interior node is known at t=0 (p=0) from prescribed initial conditions, the calculations begin at t= Δt , (p=1). With temperatures known for t=Δt, the appropriate finite-difference equation is the applied at each node to determinate its temperature at t=2 Δt (p=2). The transient temperature distribution is obtained by marching out in time, using intervals of Δt.

In order to prevent errors due to instability of the method, the prescribed value of Δt must be maintained below a certain limit, which depends on Δx and other parameters of the system. The stability criterion is determined by requiring that the coefficient associated with the node of interest at the previous time is greater than or equal to zero. The stability criterion depends on Fourier number, and for a two dimensional interior node is:



An easy example is a square plate. At first we are going to consider a material with Fourier number of 0.2, and boundary conditions with constant temperature. Temperature on the top is 20°C, and the others sides are keep constant at 100°C. The temperature at the beginning of the time is also 100°C.

The computational domain is shown below and, it represents those points in which the temperature can be calculated.


According to the domain, we are going to implement the heat equation in each node, and results are shown in the next table. We must notice that nodes (1,1) (2,1) (3,1) (4,1) (1,2) (1,3) (1,4) (4,1) (4,2) (4,3) and (4,4) can not be calculated because they are constant boundary conditions.


The result for three dimensional  heat equation with p=1500 steps is shown below.

Download EDNA

Download f90

EDNA is a powerful program which simulates the field temperature of a 3D solid body by mean of the heat diffusion equation solution.


“Fundamentals of Heat and Mass Transfer”, Incropera, sixth edition.



Introduce tus datos o haz clic en un icono para iniciar sesión:

Logo de

Estás comentando usando tu cuenta de Cerrar sesión / Cambiar )

Imagen de Twitter

Estás comentando usando tu cuenta de Twitter. Cerrar sesión / Cambiar )

Foto de Facebook

Estás comentando usando tu cuenta de Facebook. Cerrar sesión / Cambiar )

Google+ photo

Estás comentando usando tu cuenta de Google+. Cerrar sesión / Cambiar )

Conectando a %s