Finite elements/Solution of heat equation

Finite element solution for the Heat equation edit

Let us now try to create a finite element approximation for the variational initial boundary value problem for the heat equation . This problem can be stated as


Let   be a finite dimensional approximation to   and let   be a finite dimensional approximation to  . Let the weighting functions   be such that   on  . Similarly, any other function   also goes to zero on the boundary  .

Assume that the trial solutions   can be represented as


The additional quantity   results in the satisfaction of the boundary condition   on  .

If we substitute these quantities into the variational initial boundary value problem, we get the Galerkin formulation. The steps in this process are shown below.

Approximate IBVP edit

The approximate form of the variational IBVP is


Replacing   with   we get


After expanding and rearranging terms, we get


In compact notation


Similarly, the initial condition can be written as


Substituting for   and expanding out terms, we have


In compact form,


Therefore the approximate form of the variational BVP can be written as


Note that the derivatives with respect to time are still continuous in this approximate variational BVP.

Finite element approximation edit

We now discretize our domain into elements   where   and   is the total number of elements. Nodes may exist anywhere within the element domains. But the most common locations are at the element vertices and along the interelement boundaries.

Let the set of global node numbers be


Let the set nodes at which a temperature (essential BC) is prescribed be


Then the set of nodes at which a temperature is to be determined is the complement


Let   be the number of nodes in  . Then the number of equations to be solved is also  . See Figure 1 to get an idea about what these sets of nodes refer to.

Figure 1. FEM discretization for the heat conduction problem.

As we have seen before, a typical weighting function   is assumed to have the form


Note that   only if   for every  . On the other hand   for every node in  .

Similarly, a typical trial solution   can be written as


where   is the unknown temperature at node   at time  .

We also often represent the approximate form of the essential BCs in a similar way, i.e.,


Substitute equations (42), (43), and (44) into the first of equations (40). You will get


Assuming that   does not change with time, we can take the constants and time-dependent parameters outside the integrals to get


Using the usual argument about the arbitrary nature of  , we get


In compact notation,


where  , and  .

Let us rename parts of the above equation as follows:


Then we get


In matrix form,


As before, we can break up the global matrices into sums over elements. Thus,




and   is the number of nodes in an element.

We can do the same for the initial conditions. However, a more common approach is to specify the values of   directly at the nodes instead of going through an assembly process.

Computing M, K, f edit

The finite element system of equations at time  :


Isoparametric Map edit

Isoparametric map.

Coordinate Transformation edit


Integrating Stiffness Matrix edit

Stiffness matrix terms:


Index notation:


Expanded out (in 2D) (assume  ):


Transformation edit

Chain Rule:


Matrix notation:


Inverse Transformation edit




Returning to the integration of the stiffness matrix terms:


In parent element coordinates:


Gaussian Quadrature edit


Gaussian Quadrature Example edit

Find the constants Co, C1, and X so that the quadrature formula


This has the highest possible degree of precision.


Since there are three unknowns, Co, C1 and X1, we will expect the formula to be exact for




Equation 2 and 3 will yield.








Thus the degree of the precision is 2

Robin boundary conditions edit

A similar approach can be used when w:Robin boundary conditions are needed to be applied. Let us assume that the Robin boundary conditions take the form


Recall from the previous section that the boundary term in the weak form of the heat equation is




Let us assume, for simplicity, that the material is isotropic. In that case   becomes a scalar quantity   and we can write


Now, we can write the Robin boundary condition as


Using this expression we can get of the flux term in   to get


Then the boundary term takes the form


If we ignore the contributions to the trial function from the Dirichlet boundary conditions we can write


Plugging these expressions into the boundary term leads to


Now recall that spatially discretized equation for transient heat conduction can be expressed as (in simplified form; after getting rid of contributions due to Dirichlet boundary conditions and with a scalar  )


Substituting the last term on the right hand side with the new expression for the boundary term, we have


Invoking the arbitrariness of   and rearranging, we get




Then the finite element system of equations is


or, in matrix form


Note that the   matrix and the   vector are different from the case where there are no Robin boundary conditions. The boundary terms in the discretized system of equations can be computed in the usual manner and will add terms to the diagonal of the   matrix.