# Nonlinear 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)$ .

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}} 