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:

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 ($\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

${\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

← Model finite element approximation

Finite elements

Steady state heat conduction →