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,

 

where

 

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

 

with

 

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.

Solution.

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

 

Thus

 

Equation 2 and 3 will yield.

 

Hence

 

Now,

 

And

 

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

 

where

 .

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

 

Define

 

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.