# Introduction to finite elements/Linear heat equation

Finite element methods are used to solve boundary value problems (BVP) or initial boundary value problems (IBVP) in engineering. BVPs are mathematical models of real-life situations. Such situations can be physical, biological, economic, and so on.

We will explore the problem of heat conduction and see how we build a finite element model and solve this problem. The first step will be to build a model. The model is the IBVP in the form of a partial differential equation or a variational problem. We will discuss whether the problem is well-posed. Then we will try to construct an approximate solution to the problem. Finally, we will discuss how good the solution is.[1]

## Construction of a Model

Figure 1 shows a region ${\displaystyle \Omega \in \mathbb {R} ^{3}}$ through which heat is flowing. Points in the medium are represented by ${\displaystyle \mathbf {x} }$ with components ${\displaystyle (x_{1},x_{2},x_{3})}$ with respect to the basis ${\displaystyle (\mathbf {e} _{1},\mathbf {e} _{2},\mathbf {e} _{3})}$ and the origin. The heat flux into (or out of) the medium is ${\displaystyle \mathbf {q} (\mathbf {x} ,t)}$ where ${\displaystyle t}$ is the time. The heat flux takes place across part of the boundary of the medium (${\displaystyle \Gamma _{q}}$). The temperature on the remainder of the boundary (${\displaystyle \Gamma _{T}}$) has a known value ${\displaystyle {\overline {T}}}$. Heat sources inside the medium (for example, chemical reactions and plastic deformation) are given as a function ${\displaystyle s(\mathbf {x} ,t)}$.

 Figure 1. The problem of heat conduction.

The goal is to find the temperature distribution in the medium (${\displaystyle T(\mathbf {x} ,t)}$) as time evolves.

A realistic model for this problems needs two components:

1. A balance (or conservation) equation for energy.
2. A constitutive equation for the medium.

### Balance of energy

Let us assume that there are no sources of energy other than thermal energy. Then, the balance of energy states that:

Rate of change of thermal energy = heat generated by internal sources + heat flow into the body.

We have to convert this statement into mathematical form. To do this, consider an arbitrary part of the body (${\displaystyle \Omega _{1}}$) with boundary ${\displaystyle \Gamma _{1}}$. Let ${\displaystyle \mathbf {n} }$ be the outward unit normal to the boundary.

The total thermal energy in ${\displaystyle \Omega _{1}}$ at a particular time is given by the heat capacity of the body. The heat capacity ${\displaystyle C_{v}}$ is the amount of energy needed to raise the temperature (of a unit mass of the body) by one unit.

Let the mass density be ${\displaystyle \rho (\mathbf {x} )}$. Then, the total thermal energy in ${\displaystyle \Omega _{1}}$ at time ${\displaystyle t}$ is

${\displaystyle {\text{(2)}}\qquad \int _{\Omega _{1}}C_{v}(\mathbf {x} )~\rho (\mathbf {x} )~T(\mathbf {x} ,t)~dV~.}$

The total heat generated by internal sources is

${\displaystyle {\text{(3)}}\qquad \int _{\Omega _{1}}s(\mathbf {x} ,t)~dV~.}$

And the total heat flux into ${\displaystyle \Omega _{1}}$ is

${\displaystyle {\text{(4)}}\qquad -\int _{\Gamma _{1}}\mathbf {q} (\mathbf {x} ,t)\bullet \mathbf {n} ~dA~.}$

Apply the divergence theorem to (4) to get

${\displaystyle {\text{(5)}}\qquad -\int _{\Gamma _{1}}\mathbf {q} (\mathbf {x} ,t)\bullet \mathbf {n} ~dA=-\int _{\Omega _{1}}{\boldsymbol {\nabla }}\bullet \mathbf {q} (\mathbf {x} ,t)~dV~.}$

Put (2), (3), and (5) together, and apply the balance of energy to get

${\displaystyle {\text{(6)}}\qquad {\cfrac {d}{dt}}\int _{\Omega _{1}}C_{v}(\mathbf {x} )~\rho (\mathbf {x} )~T(\mathbf {x} ,t)~dV=\int _{\Omega _{1}}s(\mathbf {x} ,t)~dV-\int _{\Omega _{1}}{\boldsymbol {\nabla }}\bullet \mathbf {q} (\mathbf {x} ,t)~dV~.}$

The limits of integration are fixed. So we can write (6) as

${\displaystyle {\text{(7)}}\qquad \int _{\Omega _{1}}C_{v}(\mathbf {x} )~\rho (\mathbf {x} )~{\frac {\partial T(\mathbf {x} ,t)}{\partial t}}~dV=\int _{\Omega _{1}}s(\mathbf {x} ,t)~dV-\int _{\Omega _{1}}{\boldsymbol {\nabla }}\bullet \mathbf {q} (\mathbf {x} ,t)~dV~.}$

Since ${\displaystyle \Omega _{1}}$ is arbitrary, if the functions that appear in (7) are smooth enough, the equation is equivalent to

${\displaystyle {\text{(8)}}\qquad {C_{v}(\mathbf {x} )~\rho (\mathbf {x} )~{\frac {\partial T(\mathbf {x} ,t)}{\partial t}}+{\boldsymbol {\nabla }}\bullet \mathbf {q} (\mathbf {x} ,t)=s(\mathbf {x} ,t)~.}}$

We can write equation (8) in index notation as

${\displaystyle {C_{v}~\rho ~{\dot {T}}+q_{i,i}=s}}$

where we have assumed that the components are with respect to the Cartesian basis (${\displaystyle \mathbf {e} _{1},\mathbf {e} _{2},\mathbf {e} _{3}}$). In the following, whenever we use index notation, the components are assumed to be with respect to that Cartesian basis.

### Constitutive equation

There are two unknowns in equation (8). These are ${\displaystyle T(\mathbf {x} ,t)}$ and ${\displaystyle \mathbf {q} (\mathbf {x} ,t)}$ and only one equation. Therefore, we need another equation that characterizes the material.

One possibility is the Fourier law which states that:

the heat flux is linearly related to the temperature gradient.

In mathematical form,

${\displaystyle {\text{(9)}}\qquad {\mathbf {q} (\mathbf {x} ,t)=-{\boldsymbol {\kappa }}(\mathbf {x} )\bullet {\boldsymbol {\nabla }}T(\mathbf {x} ,t)~;~~~{\boldsymbol {\kappa }}={\boldsymbol {\kappa }}^{T}~.}}$

The quantity ${\displaystyle {\boldsymbol {\kappa }}}$ is a second-order tensor called the thermal conductivity tensor. The minus sign shows that heat flows from hot to cold. Recall that a second-order tensor takes a vector to another vector (in this case a temperature gradient to a heat flux).

In index notation, we can write equation (9) as

${\displaystyle {q_{i}=-\kappa _{ij}~T_{,j}~;~~~\kappa _{ij}=\kappa _{ji}~.}}$

If the region ${\displaystyle \Omega }$ is homogeneous then ${\displaystyle {\boldsymbol {\kappa }}}$ is constant. If the region ${\displaystyle \Omega }$ is isotropic, the thermal conductivity tensor takes the form

${\displaystyle \kappa _{ij}(\mathbf {x} )=\kappa (\mathbf {x} )\delta _{ij}}$

where ${\displaystyle \kappa }$ is a scalar.

### Heat equation

To get the heat equation, combine (9) with (8) to get

${\displaystyle {\text{(10)}}\qquad {C_{v}(\mathbf {x} )~\rho (\mathbf {x} )~{\frac {\partial T(\mathbf {x} ,t)}{\partial t}}-{\boldsymbol {\nabla }}\bullet [{\boldsymbol {\kappa }}(\mathbf {x} )\bullet {\boldsymbol {\nabla }}T(\mathbf {x} ,t)]=s(\mathbf {x} ,t)~.}}$

Rearrange to get

${\displaystyle {\text{(11)}}\qquad {{\frac {\partial T(\mathbf {x} ,t)}{\partial t}}-{\frac {1}{C_{v}(\mathbf {x} )~\rho (\mathbf {x} )}}{\boldsymbol {\nabla }}\bullet [{\boldsymbol {\kappa }}(\mathbf {x} )\bullet {\boldsymbol {\nabla }}T(\mathbf {x} ,t)]=Q(\mathbf {x} ,t)}}$

where ${\displaystyle Q(\mathbf {x} ,t):=s(\mathbf {x} ,t)/[C_{v}(\mathbf {x} )~\rho (\mathbf {x} )]}$.

In index notation, equation (11) is

${\displaystyle {{\dot {T}}-{\cfrac {1}{C_{v}~\rho }}~(\kappa _{ij}~T_{,j})_{,i}=Q~.}}$

In expanded form, equation (11) reads

{\displaystyle {\text{(12)}}\qquad {\begin{aligned}{\frac {\partial T}{\partial t}}-{\frac {1}{C_{v}~\rho }}&\left[{\frac {\partial }{\partial x_{1}}}\left(\kappa _{11}{\frac {\partial T}{\partial x_{1}}}+\kappa _{12}{\frac {\partial T}{\partial x_{2}}}+\kappa _{13}{\frac {\partial T}{\partial x_{3}}}\right)+\right.\\&~{\frac {\partial }{\partial x_{2}}}\left(\kappa _{12}{\frac {\partial T}{\partial x_{1}}}+\kappa _{22}{\frac {\partial T}{\partial x_{2}}}+\kappa _{23}{\frac {\partial T}{\partial x_{3}}}\right)+\\&\left.~{\frac {\partial }{\partial x_{3}}}\left(\kappa _{13}{\frac {\partial T}{\partial x_{1}}}+\kappa _{23}{\frac {\partial T}{\partial x_{2}}}+\kappa _{33}{\frac {\partial T}{\partial x_{3}}}\right)\right]=Q~.\end{aligned}}}

Equation (12) is the transient, inhomogeneous, heat equation.

### Boundary conditions

Boundary conditions (BCs) are needed to make sure that we get a unique solution to equation (12).

The temperature is prescribed on ${\displaystyle \Gamma _{T}}$. Prescribed boundary conditions are also called Dirichlet BCs or essential BCs. In this case

${\displaystyle {\text{(13)}}\qquad {T={\overline {T}}(\mathbf {x} ,t)~~~~{\text{on}}~~\Gamma _{T}~.}}$

The heat flux is given on the remainder of the boundary (${\displaystyle \Gamma _{q}}$). Such flux boundary conditions are also known as Neumann BCs or natural BCs. The flux boundary condition is

${\displaystyle {\text{(14)}}\qquad {\mathbf {q} (\mathbf {x} ,t)\bullet \mathbf {n} (\mathbf {x} )={\overline {q}}(\mathbf {x} ,t)~~~~{\text{on}}~~\Gamma _{q}~.}}$

In index notation, the essential boundary condition is

${\displaystyle {q_{i}~n_{i}={\overline {q}}\qquad {\text{on}}~\Gamma _{q}~.}}$

Plug equation (9) into (14). We get

${\displaystyle {\text{(15)}}\qquad -({\boldsymbol {\kappa }}\bullet {\boldsymbol {\nabla }}T)\mathbf {n} ={\overline {q}}~.}$

If the region ${\displaystyle \Omega }$ is isotropic with thermal conductivity ${\displaystyle \kappa }$, we can define ${\displaystyle g:=-{\overline {q}}/\kappa }$ and ${\displaystyle \partial T/\partial n:={\boldsymbol {\nabla }}T\bullet \mathbf {n} }$ (also called the normal derivative). Then the flux BC simplifies to

${\displaystyle {{\frac {\partial T}{\partial n}}=g~~~~{\text{on}}~~\Gamma _{q}~.}}$

If the boundary ${\displaystyle \Gamma _{q}}$ is insulated, then ${\displaystyle {\overline {q}}}$ = 0 = ${\displaystyle g}$.

### Initial conditions

If the problem depends on time, we also need an initial condition for the temperature in the body,

${\displaystyle {T(\mathbf {x} ,0)=T_{0}(\mathbf {x} )~.}}$

### The complete model

The model initial boundary value problem (IBVP) for heat conduction is

{\displaystyle {\text{(16)}}\qquad {\begin{aligned}&&{\mathsf {The~initial~boundary~value~problem~for~heat~conduction}}\\&&\\&{\text{PDE:}}~~~&~~~{\frac {\partial T}{\partial t}}-{\frac {1}{C_{v}~\rho }}{\boldsymbol {\nabla }}\bullet ({\boldsymbol {\kappa }}\bullet {\boldsymbol {\nabla }}T)=Q~~{\text{in}}~~\Omega ,t>0\quad \\&{\text{BCs:}}~~~&~~~T={\overline {T}}(\mathbf {x} ,t)~~{\text{on}}~~\Gamma _{T}~~{\text{and}}~~\mathbf {q} \bullet \mathbf {n} ={\overline {q}}~~{\text{on}}~~\Gamma _{q}\quad \\&{\text{IC:}}~~~~&~~~T(\mathbf {x} ,0)=T_{0}(\mathbf {x} )~~{\text{in}}~~\Omega \quad \end{aligned}}}


## References

1. This writeup is based on the introductory chapter in Introductory Functional Analysis with Applications to Boundary Value Problems and Finite Elements by B. Daya Reddy, Springer, 1998, and the chapter on parabolic and hyperbolic problems in The Finite Element Method: Linear Static and Dynamic Finite Element Analysis by T.J.R Hughes.