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{align} \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{align} } $$ 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{align} & &\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{align} } $$