# 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.

## Construction of a Model

Figure 1 shows a region $\Omega \in \mathbb {R} ^{3}$ through which heat is flowing. Points in the medium are represented by $\mathbf {x}$ with components $(x_{1},x_{2},x_{3})$ with respect to the basis $(\mathbf {e} _{1},\mathbf {e} _{2},\mathbf {e} _{3})$ and the origin. The heat flux into (or out of) the medium is $\mathbf {q} (\mathbf {x} ,t)$ where $t$ is the time. The heat flux takes place across part of the boundary of the medium ($\Gamma _{q}$). The temperature on the remainder of the boundary ($\Gamma _{T}$) has a known value ${\overline {T}}$. Heat sources inside the medium (for example, chemical reactions and plastic deformation) are given as a function $s(\mathbf {x} ,t)$.
* Figure 1. The problem of heat conduction.
The goal is to find the temperature distribution in the medium ($T(\mathbf {x} ,t)$) as time evolves.
A realistic model for this problems needs two components:
* A balance (or conservation) equation for energy.
* 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 ($\Omega _{1}$) with boundary $\Gamma _{1}$. Let $\mathbf {n}$ be the outward unit normal to the boundary.
The total thermal energy in $\Omega _{1}$ at a particular time is given by the heat capacity of the body. The heat capacity $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 $\rho (\mathbf {x} )$. Then, the total thermal energy in $\Omega _{1}$ at time $t$ is
    ${\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
    ${\text{(3)}}\qquad \int _{\Omega _{1}}s(\mathbf {x} ,t)~dV~.$
And the total heat flux into $\Omega _{1}$ is
    ${\text{(4)}}\qquad -\int _{\Gamma _{1}}\mathbf {q} (\mathbf {x} ,t)\bullet \mathbf {n} ~dA~.$
Apply the divergence theorem to (4) to get
    ${\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
    ${\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
    ${\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 $\Omega _{1}$ is arbitrary, if the functions that appear in (7) are smooth enough, the equation is equivalent to
    ${\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
    ${C_{v}~\rho ~{\dot {T}}+q_{i,i}=s}$
where we have assumed that the components are with respect to the Cartesian basis ($\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 $T(\mathbf {x} ,t)$ and $\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,
    ${\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 ${\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
    ${q_{i}=-\kappa _{ij}~T_{,j}~;~~~\kappa _{ij}=\kappa _{ji}~.}$
If the region $\Omega$ is homogeneous then ${\boldsymbol {\kappa }}$ is constant. If the region $\Omega$ is isotropic, the thermal conductivity tensor takes the form
    $\kappa _{ij}(\mathbf {x} )=\kappa (\mathbf {x} )\delta _{ij}$
where $\kappa$ is a scalar.

### Heat equation

To get the heat equation, combine (9) with (8) to get
    ${\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
    ${\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 $Q(\mathbf {x} ,t):=s(\mathbf {x} ,t)/[C_{v}(\mathbf {x} )~\rho (\mathbf {x} )]$.
In index notation, equation (11) is
    ${{\dot {T}}-{\cfrac {1}{C_{v}~\rho }}~(\kappa _{ij}~T_{,j})_{,i}=Q~.}$
In expanded form, equation (11) reads
    ${\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 $\Gamma _{T}$. Prescribed boundary conditions are also called Dirichlet BCs or essential BCs. In this case
    ${\text{(13)}}\qquad {T={\overline {T}}(\mathbf {x} ,t)~~~~{\text{on}}~~\Gamma _{T}~.}$
The heat flux is given on the remainder of the boundary ($\Gamma _{q}$). Such flux boundary conditions are also known as Neumann BCs or natural BCs. The flux boundary condition is
    ${\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
    ${q_{i}~n_{i}={\overline {q}}\qquad {\text{on}}~\Gamma _{q}~.}$
Plug equation (9) into (14). We get
    ${\text{(15)}}\qquad -({\boldsymbol {\kappa }}\bullet {\boldsymbol {\nabla }}T)\mathbf {n} ={\overline {q}}~.$
If the region $\Omega$ is isotropic with thermal conductivity $\kappa$, we can define $g:=-{\overline {q}}/\kappa$ and $\partial T/\partial n:={\boldsymbol {\nabla }}T\bullet \mathbf {n}$ (also called the normal derivative). Then the flux BC simplifies to
    ${{\frac {\partial T}{\partial n}}=g~~~~{\text{on}}~~\Gamma _{q}~.}$
If the boundary $\Gamma _{q}$ is insulated, then ${\overline {q}}$ = 0 = $g$.

### Initial conditions

If the problem depends on time, we also need an initial condition for the temperature in the body,
    ${T(\mathbf {x} ,0)=T_{0}(\mathbf {x} )~.}$

### The complete model

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