Heat transfer Modeling and Theory

1. What is Heat Transfer?

1.1. Heat transfer modes

Heat transfer is the energy transfer process as a result of temperature difference.

A thermal analysis is undertaken to predict temperatures and heat transfer within and around bodies. This information may then be used to model temperature dependent phenomenon such as thermally induced stresses or the effect on fluid flow in the case of a solidifying metal. Heat flow has been categorised into three different modes.


In a solid body, the energy is transferred from a high temperature region to a low temperature region. The rate of heat transfer per unit area is proportional to the temperature gradient in the normal direction. This mode of heat transfer is referred to as conduction.

Figure 1. Conduction through a solid or a stationary fluid

Convection heat transfer determines the heat flow between a solid body and the surrounding atmosphere of either gases or liquid. In the heat transfer literature, the term 'fluid' generally encompasses both gases and liquid. The convection heat transfer is either forced convection where the fluid is forced to flow around the solid body or natural convection where the fluid flow is due to density variations arising from the heat transfer process.

Figure 2. Convection from a surface to a moving fluid

In this mode of heat transfer, the heat flow can occur between a solid body and the surrounding atmosphere with or without the presence of gases or liquid. The heat flows by electromagnetic radiation. This means that heat flow can occur even if the solid body is kept in a vacuum.

Figure 3. Net radiation heat exchange between two surfaces

The heat conducted through a solid is convected as well as radiated from its surface to the atmosphere. The Figure below illustrates a particular case of this phenomenon. It is assumed that the temperature at the surface 1 is maintained at \(T_{1} K\) and is greater than the ambient temperature \(T_{\infty} K.\) As a result, the heat will conduct through the solid and dissipate into the atmosphere by the means of radiation and/or convection across the boundary surface 2.

Figure 4. Three different modes of heat transfer
\(q, q_{c}, q_{r}\)

Rate of heat transfer due to conduction, convection and radiation respectively \((W)\)


Cross sectional area normal to heat flow (at surface 1 and } 2 \(\left(m^{2}\right)\)

\(k, h, \varepsilon, \sigma\)

Conductivity \((W / m K),\) convection heat transfer coefficient \(\left(W / m^{2} K\right),\) emissivity and the Stefan Boltzmann constant \(\left(W / m^{2} K^{4}\right)\) respectively.

\(T_{1}, T_{2}, T_{\infty}\)

Temperature at surfaces 1 and \(2,\) and the ambient temperature \((K)\) respectively \(\left(T_{1}>T_{2}>T_{\infty}\right)\)

From the modelling viewpoint, the conduction equation is numerically solved within a solid domain to determine the temperature variation.

Heat transfer due to radiation and convection occur only at the boundary of the solid and are, hence, included as boundary conditions.

The boundary of the solid may also be maintained at a fixed temperature or it may be insulated. Such identification of boundary conditions is necessary in order to solve the conduction equation.

1.2. Why thermal analysis?

There are many reasons for undertaking thermal analyses. In the context of engineering and in its simplest form, thermal analysis may be used to

Estimate the temperature distribution within a component or assembly

This may be important where temperature is a key factor with regard to integrity or where differential thermal expansion causes internal stresses. Contrasting examples include microelectronic components that are subjected to a working load cycle and large structures that are subject to a fire risk.

Calculate the temperature distribution in a component

it may be possible to establish heat flows that are present where energy transfer is a primary consideration (e.g. heat exchanger and condenser).

Thermal calculations are frequently coupled with other types of analysis:

  • fluid flow and

  • structural analysis.

1.2.1. Heat transfert and Structural Snalysis

The coupling of heat transfer and structural analysis comprises a common application: A structure or a material either contracts or expands when subjected to temperature changes.

The thermal contraction or expansion is directly proportional to the temperature changes and is characterised by a coefficient of linear expansion

\[\Delta \varepsilon=\alpha \Delta T\]

where \(\Delta \varepsilon\) is the change in strain taking place in the body, \(\alpha\) is the coefficient of thermal expansion and \(\Delta T\) is the temperature change. The value of the coefficient of linear expansion is determined experimentally on very small samples and therefore represents an unrestrained system. However, in a structure, it represents initial strains in the body and can be converted into stresses and forces to be used as the load vector in the finite element structural analysis. In most practical applications, the contraction or expansion of the structure is partially or fully constrained thus inducing internal stresses. This can be important where thermally induced stresses may cause failure, where thermal distortion control is a requirement, or where residual stresses need to be evaluated in the case of components that are cooled (e.g. castings, heat treatment etc.).

1.2.2. Heat transfer and Other Physics

More complex systems also require thermal analysis:

Chemical reaction

it must incorporate a model of thermal generation that arises as a consequence of the reaction.

Phase change

Heat transfer with phase change has also become an important field because a large number of industrial applications (e.g. metal casting, plastics, injection moulding and other forming processes.) address this as part of a manufacturing process simulation. For example, the transient temperatures in a casting (or an injection moulded component) are important since they point to potential areas in the casting that may be defective, or they provide an indication of the duration of the cooling cycle.

1.3. Purpose of Thermal Analysis

Before an analysis is undertaken, it is important to establish the purpose of the investigation, since this motivates the study and defines some of the inputs and the output requirements from the analysis. * if temperature limits are a concern, then the temperature field throughout the system is the important output. * if the analysis goal is to establish cooling, then heat flux is likely to be more important and this will be displayed as an output.

A good understanding of the underlying physical principles of the thermal analysis is a basic requirement for its correct application in any modelling task, in that the model must be 'fit for purpose'.

The purpose of the analysis will also influence the input detail for the model. This is a key issue since this actually defines the thermal model. These details include the geometry, which may need to be abstracted from a design drawing or even from a sketch at the design feasibility stage, material properties and thermal boundary conditions. In the case of a transient thermal analysis this also requires the input of an initial temperature field together with time stepping information. This, similar to a structural analysis, the model input comprises abstracted geometry, material properties, constraints and loads.

The thermal analysis is classified into four categories. The classification is based on whether the heat transfer problem is time dependent and material properties or boundary conditions vary with respect to the temperature.

The simple flow diagram summarises the classification process.

Guideline to select the type of Thermal Analysis
Figure 5. Guideline to select the type of Thermal Analysis
Table 1. Summary of heat transfer processes



Rate Equation

Transport Property or Coefficient


Diffusion of energy due to random molecular motion

\(q_{x}^{\prime \prime}\left(\mathrm{W} / \mathrm{m}^{2}\right)=-k \nabla T\)

\(k(\mathrm{W} / \mathrm{m} \cdot \mathrm{K})\)


Diffusion of energy due to random molecular motion plus energy transfer due to bulk motion (advection)

\(q^{\prime \prime}\left(\mathrm{W} / \mathrm{m}^{2}\right)=h\left(T_{s}-T_{\infty}\right)\)

\(k(\mathrm{W} / \mathrm{m}^2 \cdot \mathrm{K})\)


Energy transfer by electromagnetic waves

\(\begin{aligned} q^{\prime \prime}\left(\mathrm{W} / \mathrm{m}^{2}\right) &=\varepsilon \sigma\left(T_{s}^{4}-T_{\mathrm{sur}}^{4}\right) \\ \text { or } q(\mathrm{W}) &=h_{r} A\left(T_{s}-T_{\mathrm{sur}}\right) \end{aligned}\)

\(\begin{aligned}\varepsilon\\ h_{r}\left(\mathrm{W} / \mathrm{m}^{2} \cdot \mathrm{K}\right)\end{aligned}\)

2. Nomenclature and Units

Table 2. Notations and units
Notation Quantity Unit



\(m^{2}=10^{6} \mathrm{mm}^{2}=10^{4} \mathrm{cm}^{2}\)



\(k g / m^{3}=10^{3} g / c m^{3}\)



\(S\) (second)



\(C\) (degree Celsius), \(K\) (degree Kelvin)


Specific heat

\(J / k g C=J / k g K\)



\(W / m C=W / m K\)


Convection heat transfer coefficient

\(W / m^{2} C=W / m^{2} K\)


Rate of heat transfer



Stefan Boltzmann constant

\(W / m^{2} K^{4}\)



Table 3. Typical values of the convection heat transfer coefficient


\(h\) \(\left(\mathrm{W} / \mathrm{m}^{2} \cdot \mathrm{K}\right)\)

Free convection





Forced convection





Convection with phase change

Boiling or condensation


3. Heat Transfer Theory


3.1. Equations

heat equation
\[\rho C_p \frac{\partial T}{\partial t} - \nabla \cdot \left( k \nabla T \right) = Q, \quad \text{ in } \Omega\]

which is completed with boundary conditions and initial value

\[\text{at } t=0, \quad T(x,0) = T_0(x)\]

3.1.1. Convective heat transfer

convective heat equation
\[\rho C_p \left( \frac{\partial T}{\partial t} + \boldsymbol{u} \cdot \nabla T \right) - \nabla \cdot \left( k \nabla T \right) = Q, \quad \text{ in } \Omega\]

3.1.2. Steady case

steady heat equation
\[ - \nabla \cdot \left( k \nabla T \right) = Q, \quad \text{ in } \Omega\]
steady convective heat equation
\[\rho C_p \boldsymbol{u} \cdot \nabla T - \nabla \cdot \left( k \nabla T \right) = Q, \quad \text{ in } \Omega\]

3.1.3. Multi-materials

Given a domain \(\Omega \subset \mathbb{R}^d, d=1,2,3\), \(\Omega\) is partitioned into \(N_r\) regions \(\Omega_i,i=1,\ldots,N_r\) corresponding to different materials (solid or fluid). We consider \(\rho_i\), \(C_{p,i}\) and \(k_i\) the material properties defined in each regions \(\Omega_i\). We define also \(\boldsymbol{n}_i\) the outward unit normal vector associated to the boundary \(\partial \Omega_i\).

\[\begin{eqnarray} \rho_i C_{p,i} \frac{\partial T}{\partial t} - \nabla \cdot \left( k_i \nabla T \right) &=& Q, \quad &\text{ in }& \Omega_i,i=1,\ldots,N_r \\ T_{|_{\Omega_i}} &=& T_{|_{\Omega_j}}, \quad &\text{ on }& \partial \Omega_i \cap \Omega_j = \Gamma_{ ij}, \forall i \neq j \\ -k_i \nabla T \cdot \boldsymbol{n}_i &=& k_j \nabla T \cdot \boldsymbol{n}_j, \quad &\text{ on }& \partial \Omega_i \cap \Omega_j = \Gamma_{ ij}, \forall i \neq j \\ \end{eqnarray}\]

We assume the operator \(\mathcal{L}\) tel que \(\mathcal{L} T = \rho_i C_{p,i} \frac{\partial T}{\partial t} - \nabla \cdot \left( k_i \nabla T \right)\) is elliptical.

We multiply \(\mathcal{L} u = Q\) by a function test \(v \in \mathbf{V}\) and integrates by part on \(\Omega_i\). Which give:

\[\rho C_{p,i} \displaystyle \int_{\Omega_i} \frac{\partial T}{\partial t} v - \int_{\Omega_i} \nabla \cdot \left[ k_i \nabla T \right] v = \int_{\Omega_i} Qv, \quad \forall v \in H^1_0(\Omega) \quad for i = 1, \cdots , N_r\]

By the formula of Green, we get

\[\rho C_{p,i} \displaystyle \int_{\Omega_i} \frac{\partial T}{\partial t} v + \int_{\Omega_i} k_i(y) \nabla T \cdot \nabla v- \int_{\partial \Omega} k_i \nabla T \cdot \boldsymbol{n}_i v = \int_{\Omega_i} Qv \quad \forall v \in \mathbf{V}\]

Additivity of the integral, we have

\[\sum_{ i=1}^{N_r} \left( \rho C_{p,i} \displaystyle \int_{\Omega_i} \frac{\partial T}{\partial t} v + \int_{\Omega_i} k_i \nabla T \cdot \nabla v- \int_{\partial \Omega_i} k_i \nabla T \cdot \boldsymbol{n}_i v \right) = \sum_{ i=1}^{N_r} \left( \int_{\Omega_i} Qv \right) \forall v \in \mathbf{V}\]

Note that

\[\bigcup_{ i=1}^{ N } \partial \Omega_i = \bigcup_{ i,j} \Gamma_{ ij} \cup \partial \Omega\]

Use the conditions in the interfaces, we get

\[\sum_{ i=1}^{N_r} \left( \rho C_{p,i} \displaystyle \int_{\Omega_i} \frac{\partial T}{\partial t} v + \int_{\Omega_i} k_i \nabla T \cdot \nabla v- \int_{\partial \Omega} k_i \nabla T \cdot \boldsymbol{n} v \right) = \sum_{ i=1}^{N_r} \left( \int_{\Omega_i} Qv \right) \forall v \in \mathbf{V}\]

Using the implicit Euler method for the time term:

\[\frac{\partial T}{\partial t} (t^{ k+1}) \approx \frac{ T (t^{ k+1}) - T(t^k)}{ dt} \quad \forall t^k \in \mathbb{ R^+} \text{ et } k \in \mathbb{N}\]

Denoting \(T^k = T(t^k)\), we write the formula in \(t^{ k+1}\), we obtain:

\[\sum_{ i=1}^{N_r} \left( \rho C_{p,i} \displaystyle \int_{\Omega_i} \frac{ T^{k+1}}{dt} v + \int_{\Omega_i} k_i \nabla T^{k+1} \cdot \nabla v - \int_{\partial \Omega} k_i \nabla T^{k+1} \cdot \boldsymbol{n} v \right) = \sum_{ i=1}^{N_r} \left( \int_{\Omega_i} \frac{T^{k}}{dt} v + \int_{\Omega_i} Qv \right) \quad \forall v \in \mathbf{V}\]

So, the weak wording becomes:

The weak formulation
\[\text{ On cherche } T \in \mathbf{H} \text{ telle que:} \\ a(T^{k+1}, v) = l(v) \quad \forall v \in \mathbf{V} \\ \text{ and} \quad a(T^{k+1}, v) = \sum_{ i=1}^{N_r} \left( \rho C_{p,i} \displaystyle \int_{\Omega_i} \frac{ T^{k+1}}{dt} v + \int_{\Omega_i} k_i \nabla T^{k+1} \cdot \nabla v - \int_{\partial \Omega} k_i \nabla T^{k+1} \cdot \boldsymbol{n} v \right) \\ l(v) = \sum_{ i=1}^{N_r} \left( \int_{\Omega_i} \frac{T^{k}}{dt} v + \int_{\Omega_i} Qv \right)\]

So we have \(a(u_{k+1},v)\) a continuous bilinear form coercive in \(v \in \mathbf{V}\) and \(l(\phi)\) a continuous linear form . We are in a Hilbert space, so we have all the conditions for the application of the Lax-Milgram theorem. So this problem is well posed.

Correct approximation:

We use the Galerkin approximation method:

Let \(\{ \mathcal{T}_h \}\) a family of meshes of \(:\Omega\).

Let \(\{ \mathcal{K}, P, \sum \}\) a finite element of Lagrange of reference of the degree \(k \geq 1\).

Let \(P^k_{c,h}\) the conforming approximation space defined by

\[P^k_{ c,h} = \{ v \in C^0(\Omega), \forall \mathcal{K} \in \mathcal{T}_h, v|_{\mathcal{K}} \in \mathbb{P}_k(\mathcal{K}) \}\]

To obtain a conformal approximation in V, we add the boundary conditions

\[V_h = P^k_{c,h} \cap V\]

Discrete problem is written:

Problème discrète
\[\text{ Find } T_h \in V_h \text{ such that} \\ a(T_h, v_h) = l(v_h) \quad \forall v_h \in V_h\]

Let \(\{ \varphi_1, \varphi_2, ..., \varphi_N \}\) the base of \(V_h\). An element \(T_h \in V_h\) is written as

\[T_h = \sum^{N}_{l=1} T_l \varphi_l\]

Using \(v\) as a basic function of \(V_h\), our problem becomes

\[\sum_{ i=1}^{N_r} \left( \rho C_{p,i} \displaystyle \int_{\Omega_{i}} \sum_{ l=1}^N T^{k+1}_l \frac{ \varphi_l }{dt} \varphi_j + \int_{\Omega_i} k_i \sum_{ l=1}^N T^{k+1}_l \nabla \varphi_l \cdot \nabla \varphi_j - \int_{\partial \Omega} k_i \sum_{ l=1}^N T^{k+1}_l \nabla \varphi_l \cdot \boldsymbol{n} \varphi_j \right) = \sum_{ i=1}^{N_r} \left( \int_{\Omega_i} \sum_{ l=1}^N T^{k}_l \frac{ \varphi_l }{dt} \varphi_j + \int_{\Omega_i} Q \varphi_j \right)\]

The variational problem of approximation is then equivalent to a linear system

Algebraic problem
\[\text{Determine } T_l \text{ satisfying} \\ \sum_{ l=1}^N a(\varphi_l, \varphi_j) T^{k+1}_l = l(\varphi_j) \forall j = 1, \cdots , N\]


\[A = (a(\varphi_i , \varphi_j)), \quad 1 \leq i,j \leq N , \\ U^{k+1} = (T_1^{k+1}, T_2^{k+1}, ..., T_N^{k+1}) \in \mathbb{R}^{N}, \\ F = (l(\varphi_1), l(\varphi_2), ..., l(\varphi_N)) \in \mathbb{R}^{N}\]

We write the system in matrix form

\[AU = F\]

3.1.4. Variational formulation and discretization of the heat equation with radiative boundary conditions on several surfaces

Radiative heat transfer is not yet available in the toolboxes. An application implementing radiative heat is currently available in feelpp/doc/manual/heat.

Let \(\partial\Omega_D\) and \(\partial \Omega_N\) be the portions of the boundary where Dirichlet and Neumann boundary conditions are applied, respectively. Let us write the variational formulation of the heat equation: find \(T \in H^1((0,T);H^1(\Omega))\) such that, for all \(\phi \in H^1_{0,\partial \Omega_D}(\Omega)\)

\[\int_\Omega \rho c_p \partial_t T \phi + \int_\Omega k \nabla T \cdot \nabla \phi + \int_{\partial \Omega_R} k \nabla T \cdot \vec{n} \phi = \int_\Omega S \phi - \int_{\partial \Omega_N} k \nabla T \cdot \vec{n} \phi.\]

When the radiative boundary \(\partial \Omega_R\) is composed of several subsurfaces that can exchange heat through radiation, the associated radiative boundary condition is complex. In fact, each surface receives heat contributions from the other ones, proportionally to the values of the corresponding view factors. From Modest’s book [Radiative_heat_transfer], equation (5.28), the radiative heat flux at point \(x\) of \(\partial \Omega_R\) is

\[\frac{q(x)}{\epsilon(x)} - \int_{\partial \Omega_R} (\frac{1}{\epsilon(x')}-1)q(x')dF_{x-x'} = E_b(x) - \int_{\partial \Omega_R} E_b(x') dF_{x-x'}.\]

In this equation, \(q=\nabla T \cdot \vec{n}, E_b(x)=\sigma T^4\), \(\epsilon(x)\) is the emittance and \(dF_{x-x'}\) is the view factor between the infinitesimal areas surrounding points \(x,x'\).

Let us now propose a variational formulation for this equation. Let \(\psi \in \mathbb{P}^0_{d}(\partial \Omega_R)\) be discontinuous, piecewise constant basis functions. Functions \(\psi\) are elementwise discontinuous; however, one could choose to work, as a first approximation, with \(\Psi \in P\), where \(P\) a space of discontinuous, piecewise constant basis functions, where discontinuities are not elementwise, but for example discontinuous in correspondence of different radiating surfaces. In the following, we will use \(\psi\) to denote test functions and \(N_h\) to denote the cardinality of the test space. We have

\[\begin{multline} \int_{\partial \Omega_R} \frac{q(x)}{\epsilon(x)} \psi_j - \int_{\partial \Omega_R}\int_{\partial \Omega_R} (\frac{1}{\epsilon(x')}-1)q(x')dF_{x-x'} \psi_j \\ = \int_{\partial \Omega_R} E_b(x) \psi_j- \int_{\partial \Omega_R} \int_{\partial \Omega_R} E_b(x') dF_{x-x'} \psi_j. \end{multline}\]

By decomposing \(q(x)=\sum_{i=0}^{N_h} q_i \psi_i(x)\), the first term of the left-hand side gives rise to a mass matrix \(M_{ij} = \int_{\partial \Omega_R} \epsilon_i^{-1} \psi_i\psi_j\). The non-linearity in \(q\) is handled iteratively: the second term on the left-hand side is treated explicitly and moved to the right-hand side. It can be decomposed as a function of coefficients \(q_i\) as \(N_{j}(q) = \int_{\partial \Omega_R} \Big( \sum_{A_k} \frac{1}{A_k} \int_{A_k} (\frac{1}{\epsilon_k}-1) q_k F_{ik} \, dA \Big)\psi_j\). The right-hand side is of the form \(D_j(T) = \int_{\partial \Omega_R} ( \sigma T^4 + \sum_{A_k} \frac{1}{A_k} \int_{A_k} T^4_k F_{jk} \, dA) \psi_j\).

Due to the presence of the fourth power of the temperature and the non-linearity of the equation with respect to temperature and heat flux, two nested iterative loops are proposed. In the following algorithm, \(n\) denotes time indices, \(k\) denotes the indices of the temperature loop and \(l\) denotes the indices of the flux loop.

Double iterative loop
    \caption{Solution of the heat equation with radiative BC on the timestep $\Delta t$.}
    \STATE $T^{n}=T^{n,0}$, $q^{n}=q^{n,0,0}$
    \STATE $T^{n+1,0} = Heat_{\Delta t}(T^n,q^{n})$;
    \STATE $q^{n+1,0,0} = q^{n}$
    \WHILE{$||T^{n+1,k+1} -T^{n+1,k}||/||T^{n+1,0}|| > \tau_T$}
        \STATE $T^{n+1,k} \leftarrow T^{n+1,k+1}$
        \STATE $q^{n+1,k,0} = M_{ij}^{-1} (D_j(T^{n+1,k})-N_{ij}(q^{n+1,k+1,0}))$
        \WHILE{$||q^{n+1,k,l} -q^{n+1,k,l-1}||/||q^{n+1,k,0}|| > \tau_q$}
            \STATE $q^{n+1,k,l-1} \leftarrow q^{n+1,k,l}$
            \STATE $q^{n+1,k,l} = M_{ij}^{-1} (D_j(T^{n+1,k})-N_{j}(q^{n+1,k,l-1}))$
        \STATE $q^{n+1,k+1,0} \leftarrow q^{n+1,k,l}$
        \STATE $T^{n+1,k+1} = Heat_{\Delta t}(T^{n+1,k},q^{n+1,k+1,0})$
    \STATE $T^{n+1} \leftarrow T^{n+1,k+1}$
    \STATE $q^{n+1} \leftarrow q^{n+1,k+1,0}$

3.2. Bibliography

[Radiative_heat_transfer] Modest, M.F., Radiative Heat Transfer, Elsevier Science (2013) doi.org/10.1016/C2010-0-65874-3



Conductivity of a solid is the heat flow per unit area in the normal direction within the solid body as a result of temperature difference.

Heat flux(conduction)

The heat flux by conduction \(q''(W/m2)\) is the heat transfer rate per unit area perpendicular to the direction of transfer, and it is proportional to the temperature gradient in this direction. We have the Fourier law in the multi-dimensional case:

\[q'' = -k \nabla T\]
Heat Rate(conduction)

The heat rate by conduction, \(q (W)\), through a plane wall of area \(A\) is then the product of the flux and the area, \(q = q'' \cdot A\).

Specific heat

It is the energy required to raise the temperature of unit mass through one degree temperature rise.

Adiabatic condition

An adiabatic condition means that heat is not allowed to transfer across the boundary. In other words it means a perfectly insulated boundary. Symmetry boundaries are adiabatic.


This terminology is used in the context of radiation. It is the ratio of the total energy emitted by a surface to the total energy emitted by a black surface at the same temperature.

View factor

In the radiation heat transfer context, it is the ratio of the energy incident on the second body to the energy emitted by the first body. Sometimes, it is also referred to as a geometric factor or form factor.

Thermal stress

Stress induced in a structure or body as a result of thermal loading.

Phase change

Substance subjected to temperature changes may transfer from solid state to liquid and/or gaseous state. During change of phase the latent heat is either released or absorbed.