#### Numerical Algorithms for Simulating Three Modes of Heat Transfer

By Charles Xie ✉

Energy2D is an open-source interactive multiphysics simulation tool that I have been developing for research and education purposes. It can be used to simulate all the three modes of heat transfer: conduction, convection, and radiation. The three modes are also coupled to some extent (a completely accurate coupling like what Mother Nature does may be too hard to achieve). This article opens the black box and describes the numerical algorithms under the hood of Energy2D. It introduces only the algorithms for heat transfer and fluid dynamics. For the coupling with particle dynamics, see this article.

##### Simulating thermal conduction

Figure 1: Comparing conduction from sources at different temperatures

The thermal conduction for heterogeneous media is modeled using the Heat Equation as follows:

Equation 1

where *T*(*x*, *y*, *z*, *t*) is the temperature distribution, *k* is the
thermal conductivity, *c* is the
specific heat capacity,
*ρ* is the density, **v** is the velocity field, and *q*(*x*, *y*, *z*, *t*) is the internal heat generation
(as opposed to external heating implemented as a
boundary condition). The Heat Equation has two parts: the
diffusion part characterized by the conductivity field function
*k*(*x*, *y*, *z*) and the
advection part characterized by the velocity field function **v**(*x*, *y*, *z*, *t*).
The internal heat generation function *q*(*x*, *y*, *z*, *t*) can be thought of as a driving force for diffusion and advection.
Equation 1 is discretized in both space and time domains through a
finite-difference time-domain (FDTD) method.
An implicit FDTD method is used to achieve better numerical stability:

Equation 2

where *M*>*i*≥0, *N*>*j*≥0, and *L*>*k*≥0, with *M*, *N*, and *L* representing the numbers of
cubic voxels in *x*, *y*, and *z* directions,
and Δ*x*, Δ*y*, and Δ*z* are the lengths of the cubic voxels in the three directions, respectively. Δ*t* is the time step.
In theory, an implicit method is *unconditionally stable*,
allowing any time step to be taken, and the simulation will not “blow” up.
In practice, however, the advection part of the Heat Equation can cause numeric instability, especially when the velocity is high.
What is more, a numerically stable solver does not necessarily produce scientifically accurate results.
Note that in Equation 2 we also avoid directly calculating the spatial derivatives of the conductivity field function,
which could be very large at the interfaces between two distinct types of media and therefore could result in undesirable numeric instability.

Let:

Equation 3

Equation 2 can be rewritten as:

Equation 4

This is a very sparse matrix that can be solved efficiently using a relaxation method.

Several types of boundary conditions can be used. The first kind is the Dirichilet boundary condition, which maintains the temperature at the boundary. The second kind is the Neumann boundary condition, which maintains the flux at the boundary. Mixed boundary conditions may also be applied.

##### Simulation thermal convection

Figure 2: Comparing convection from sources at different temperatures

Fluids can transfer heat through forced and natural convection. Convection differs from conduction in that heat moves along with molecules in convection but molecules exchange only momenta — not positions — in conduction. In Energy2D, we only consider incompressible Newtonian fluids, which can be modeled by the following form of the Navier-Stokes equation:

Equation 5

where **v** is the velocity vector, *a* is the kinematic viscosity, *p* is the pressure, and **f** is the body force such as gravity or buoyancy.
The first term on the right-hand side of Equation 5 describes the diffusion of momentum. The second term describes advection.
The continuity condition
ensures the conservation of mass (as our fluids are incompressible, the mass in every voxel must not change):

Equation 6

For the temperature distribution to affect the flow of the fluid, we can add a thermal buoyancy term, known as the
Boussinesq approximation, to the body force **f**(*x*, *y*, *z*, *t*):

Equation 7

where *T*_{0}(*x*, *y*, *t*) is the average temperature of the column of fluid at the location (*x*, *y*) and at time *t* and *g* is the
thermal expansion coefficient.
The Navier-Stokes equation (Equation 5) can be decomposed into two steps (assuming a uniform distribution of pressure or incorporating the pressure gradient term into the body force term):

**1. The diffusion step**

Equation 8

**2. The advection step**

Equation 9

The diffusion step Equation 8 can be solved efficiently using a relaxation method like in the case of solving the Heat Equation, whereas the advection step Equation 9 can be solved using the MacCormack method, which consists of the following two steps:

**1. The predictor step**

Equation 10

**2. The corrector step**

Equation 11

The function *ξ* can be any component of the velocity vector (*u*, *v*, *w*).

The solution of Equation 5 does not automatically satisfy the conservation of mass. The mass continuity condition expressed in Equation 6 must be imposed at each step. This condition is also important in simulating convective loops as it forces the velocity fields to have vortices. The discretized form of the continuity equation requires that:

Equation 12

Using the Helmholtz decomposition (the fundamental theorem of vector calculus), a vector field can be decomposed into the following form:

Equation 13

where **V**’ is a divergence-free vector field

Equation 14

and *ϕ* is a scalar field. Applying the gradient operator to both sides of Equation 13, we obtain:

Equation 15

This is a Poisson equation for the scalar field *ϕ*, given that **V** is known.
If *ϕ* can be solved, then we can obtain the divergence-free solution by subtracting ∇*ϕ* from **V**
(see Modeling the Motion of a Hot, Turbulent Gas by Foster and Metaxas):

Equation 16

Equation 16 can be discretized as follows:

Equation 17

Again, the Poisson equation can be solved by using a relaxation method. The method uses the following iteration until convergence is reached:

Equation 18

A multigrid method may be used to accelerate the convergence. If accuracy is less important, such is in the case of producing certain visual effects for animations or games, 5–20 iteration steps can be used without worrying about whether or not the solution has truly converged.

##### Simulation thermal radiation

Figure 3: Comparing radiation from sources at different temperatures

Heat can transfer at the speed of light through thermal radiation.
The Stefan-Boltzmann Law describes how the power *E*
emitted from an object depends on its temperature:

Equation 19

where *σ* is called the Stefan constant, *ε* is the emissivity
(1 if the object is a perfect black body and 0 if it is a perfect reflector), and *T* is the temperature in Kelvin.

For a system that has many objects, we have to discretize every one of them into many small patches {d*A*_{i}}
in order to model the radiative heat transfer among them. The power density given away by each patch is called
radiosity (denoted by *B*),
which is governed by the following integral equation
(this is in fact an inhomogeneous Fredholm equation of the second kind):

Equation 20

where *B*_{i} is the outgoing power per unit area from patch *i*, *E*_{i} is the emission power per unit area of the patch,
*ρ*_{i} is the reflectivity of the patch, d*A*_{i} is the area of the patch, *F*_{ji} is the
view factor (dimensionless),
and *V*_{ji} is the visibility function (0 if patches *i* and *j* cannot “see” each other, otherwise 1).
The visibility function can be calculated using collision detection
between the rays emitted by the source patch and the target patch, which means that the radiation solver can handle both the
convex and
concave shapes.

The view factor from a surface *A*_{1} to a surface *A*_{2} defines the proportion of the radiation leaving *A*_{1} that ends up striking
*A*_{2}. In 3D, it can be calculated through the following formula based on the etendue of a ray:

Equation 21

where *r*_{12} is the distance between the two surfaces, and *θ*_{1} and *θ*_{2} are the acute angles between the surface normal and the distance vector.
In 2D, it can be calculated by replacing the areas with lengths and *r*^{2}_{12} with *r*_{12}:

Equation 22

The net energy that converts into thermal energy at a patch comes from the difference between the portion of energy reflected from all other patches that becomes absorbed by the patch and the energy that it emits:

Equation 23

where *α*_{i} is the absorptivity (the ability to absorb electromagnetic radiation) of the patch’s material.
This equation can be coupled with the Heat Equation and the Navier-Stokes Equation to allow radiation to interact with the flow of heat and mass.
Equations 21 and 23 can be discretized into the following matrix forms:

Equation 24

Equation 25

Equation 24 can be solved using a relaxation method, such as the
Gauss-Seidel method, using the emission power at each patch as the initial guess.
In a transient simulation, only a few relaxation steps are needed as the time-dependent solution is already somehow iterative.
The thermal power generated at each patch *Q*_{i} can then be derived using Equation 25.